Venomic Analysis of the Poorly Studied Desert Coral Snake, Micrurus tschudii tschudii, Supports the 3FTx/PLA2 Dichotomy across Micrurus Venoms

The venom proteome of the poorly studied desert coral snake Micrurus tschudii tschudii was unveiled using a venomic approach, which identified ≥38 proteins belonging to only four snake venom protein families. The three-finger toxins (3FTxs) constitute, both in number of isoforms (~30) and total abundance (93.6% of the venom proteome), the major protein family of the desert coral snake venom. Phospholipases A2 (PLA2s; seven isoforms, 4.1% of the venom proteome), 1–3 Kunitz-type proteins (1.6%), and 1–2 l-amino acid oxidases (LAO, 0.7%) complete the toxin arsenal of M. t. tschudii. Our results add to the growing evidence that the occurrence of two divergent venom phenotypes, i.e., 3FTx- and PLA2-predominant venom proteomes, may constitute a general trend across the cladogenesis of Micrurus. The occurrence of a similar pattern of venom phenotypic variability among true sea snake (Hydrophiinae) venoms suggests that the 3FTx/PLA2 dichotomy may be widely distributed among Elapidae venoms.


Introduction
New World genus Micrurus (Elapidae) (Wagler, 1824) [1] represents a monophyletic clade of some 80 currently recognized species of venomous coral snakes [2][3][4][5][6][7][8], although the topology of the tree is still unresolved. Coral snakes are widely distributed in tropical and subtropical regions from the southern United States to northeastern Argentina, including several continental islands inhabited by endemic forms [2,3]. Coral snakes are considered by herpetologists to be among the most beautiful snakes of the planet, as they are adorned with unique combinations of red-, black-, and yellow-colored banding. Based on their body ring color pattern characteristics, tail proportion with respect to body length, and hemipenial morphology, coral snakes are traditionally arranged in four species groups: a tricolored monadal group, a bicolored group, a Central American tricolored triadal group, and a South American triadal group [3,[9][10][11][12]. The monadal group comprises long-tailed species with a single black ring separating the white and red rings; the bicolored group contains short-tailed species with white or red rings separated by black rings; the Central American triadal group taxa are long-tailed species with three black rings separated by white rings between red rings; the South American triadal group is represented by short-tailed species with the same color pattern as the Central American triadal group species but is restricted geographically to South America, extending from Panama to southern Argentina. In the last decade, the increasing application of omics techniques to the study of snake venoms has greatly enhanced our knowledge on their composition, evolution, biological activities, and clinical effects [13][14][15][16][17]. However, the venoms of only a handful of the vast number (~130) of species and subspecies that constitute the genus Micrurus have been the subject of proteomic studies (consult Table 2 of [18]). To understand this fact it should be taken into consideration that, despite producing among the most potent neurotoxic venoms of any New World snake, bites and fatalities by coral snakes are very rare. This is facilitated by the fact that coral snakes are not aggressive; when confronted by humans, coral snakes will almost always attempt to flee, and bite only as a last resort. On the other hand, coral snakes generally inhabit sparsely populated areas and are thus infrequently encountered, and they have short fangs that cannot penetrate thick leather clothing. Consequently, envenomings by coral snakes are less frequent than those produced by sympatric pitvipers, by far the most dangerous snakes of South America, accounting for less than 3% of the snakebites recorded in the Americas [19][20][21].
The desert coral snake M. tschudii tschudii (January 1858) [22], one of the lesser-studied species of the genus Micrurus, belongs to the South American triad coral snakes [11,12,23]. Named after the Swiss naturalist, explorer and diplomat Johann Jakob von Tschudi (1818-1889), who traveled extensively in South America and published important works in herpetology, M. t. tschudii is a small (adults average 45 to 55 cm in length), tri-colored (Figure 1), mainly diurnal and terrestrial coral snake, found in tropical deciduous forest, dry tropical forest, and thorn scrub, mainly along watercourses, from near sea level to 1450 m elevation, in the western slopes of the Andes on the semi-arid Pacific coast of South America, from southern Ecuador to southwestern Peru [2,3,24]. The desert coral snake feeds on geckkonid lizards (Phyllodactilus spp.), Amphisbaenids (A. occidentalis), and colubrids of the species Mastigodryas haathi [3]. However, available literature about its reproduction and activity is scarce, and in particular virtually nothing (not even a single entry in PubMed) is known about its venom. Here, we have applied a venomics approach [13,14] to gain an insight into the spectrum of toxins that make up the venom proteome of M. t. tschudii.

The Venom Proteome of the Desert Coral Snake
The venom of M. t. tschudii was fractionated into 26 RP-HPLC fractions ( Figure 1A). Each chromatographic fraction was analyzed by SDS-PAGE ( Figure 1B), and the protein bands were excised and submitted to mass spectrometric analysis. The MS/MS data, listed in Supplementary Table S1, resulted in the identification of ě38 proteins belonging to only four snake venom protein families: three-finger toxin (3FTx), phospholipase A 2 (PLA 2 ), Kunitz-type, and L-amino acid oxidase (LAO) ( Table 1). The relative abundances of M. t. tschudii venom proteins and toxin families are displayed, respectively, in Table S1 and Figure 1C. The 3FTxs constitute, by far, both in number of isoforms (~30) and total abundance (93.6% of the total venom proteins), the major protein family of the desert coral snake venom. Seven PLA 2 isoforms, 1-3 Kunitz-type proteins, and 1-2 LAO molecules (accounting, respectively, for 4.1%, 1.6%, and 0.7% of the venom proteome) complete the toxin arsenal of the desert coral snake.

The Venom Proteome of the Desert Coral Snake
The venom of M. t. tschudii was fractionated into 26 RP-HPLC fractions ( Figure 1A). Each chromatographic fraction was analyzed by SDS-PAGE ( Figure 1B), and the protein bands were excised and submitted to mass spectrometric analysis. The MS/MS data, listed in Supplementary Table S1, resulted in the identification of ≥38 proteins belonging to only four snake venom protein families: three-finger toxin (3FTx), phospholipase A2 (PLA2), Kunitz-type, and L-amino acid oxidase (LAO) ( Table 1). The relative abundances of M. t. tschudii venom proteins and toxin families are displayed, respectively, in Table S1 and Figure 1C. The 3FTxs constitute, by far, both in number of isoforms (~30) and total abundance (93.6% of the total venom proteins), the major protein family of the desert coral snake venom. Seven PLA2 isoforms, 1-3 Kunitz-type proteins, and 1-2 LAO molecules (accounting, respectively, for 4.1%, 1.6%, and 0.7% of the venom proteome) complete the toxin arsenal of the desert coral snake.
The 3FTx and PLA2 molecules are hallmark components of the venoms of Elapidae. Catalytically active PLA2 molecules typically exhibit presynaptic neurotoxic activity, myotoxic activity, or both, although some forms exhibit antiplatelet activity [25,26]. Despite their pronounced structural similarity, members of the 3FTx family exhibit a wide variety of pharmacological effects including postsynaptic neurotoxicity, cytotoxicity, cardiotoxicity, and anticoagulant, and antiplatelet activities [27,28]. In addition, L-type Ca 2+ channel antagonists of the 3FTx family may act synergistically with muscarinic three-finger toxins to promote hypotension [29]. Kunitz-type serine protease inhibitor isolated from elapid venoms blocked the activity of a range of serine proteases, producing an antihemorrhagic effect [30]. Non-covalent PLA2-KUN complexes have been characterized from The 3FTx and PLA 2 molecules are hallmark components of the venoms of Elapidae. Catalytically active PLA 2 molecules typically exhibit presynaptic neurotoxic activity, myotoxic activity, or both, although some forms exhibit antiplatelet activity [25,26]. Despite their pronounced structural similarity, members of the 3FTx family exhibit a wide variety of pharmacological effects including postsynaptic neurotoxicity, cytotoxicity, cardiotoxicity, and anticoagulant, and antiplatelet activities [27,28]. In addition, L-type Ca 2+ channel antagonists of the 3FTx family may act synergistically with muscarinic three-finger toxins to promote hypotension [29]. Kunitz-type serine protease inhibitor isolated from elapid venoms blocked the activity of a range of serine proteases, producing an antihemorrhagic effect [30]. Non-covalent PLA 2 -KUN complexes have been characterized from venom of some Micrurus species [18 and references cited]. These complexes, originally discovered in the venom of the Texas coral snake M. tener, are target acid-sensing receptors ASIC1a/2 evoking pain [31,32]. There are many examples of venom-derived toxins that elicit notoriously intense pain [33]. Pain response may cause paralysis and serve as a warning signal for discouraging potentially threatening predators by triggering lasting acute physiological distress. The L-amino acid oxidases are flavoenzymes that catalyze oxidative deamination of L-amino acids to form corresponding α-keto acids, hydrogen peroxide and ammonia. LAOs are thought to contribute to the toxicity of the venom due to the production of hydrogen peroxide during the oxidation reaction [34,35], although their contribution to the envenoming process remains elusive.
The low yield of venom and the challenge of long-term maintenance of the desert coral snake in captivity precluded a detailed toxicovenomics analysis of the venom components of this poorly studied Micrurus species.
Margres and coworkers [44] have detected strong evidence of positive selection for the 3FTx and PLA 2 toxin families of M. fulvius. This accelerated evolution is most likely due to their direct involvement in fitness. Figure 2 shows a clear geographical distribution of PLA 2 -and 3FTx-predominant Micrurus venoms along their north-south dispersal axis, further supporting the adaptive nature of this phenotypic dichotomy.  (Table 1). Distribution ranges were adapted from [3] and The Reptile Database (http://www.reptile-database.org), and are color-coded: green, PLA2rich venom phenotype; red, 3FTx-predominant venom composition. The arrow highlights the trend towards diverging venom phenotypes along the Micrurus north-south dispersal, suggesting the epicenter of the divergence in Mesoamerica.
The evolutionary origin and adaptive relevance of the puzzling 3FTx/PLA2 dichotomy remains elusive since 3FTx-and PLA2-predominant venoms are scattered through the phylogenetic tree of Micrurus. The dominant protein families in the venom proteome of C. bivirgata flaviceps are PLA2 (41%), 3FTx (22.6%) and SVMP (18.7%) [65], suggesting that the ancestor venom phenotype of coral snakes might have been of the PLA2-predominant type. Alternatively, 3FTx-and PLA2-rich venom proteomes may represent pedomorphic and ontogenetic traits, as has been documented in Crotalus. Rattlesnake venoms belong to one of two distinct phenotypes, which broadly correspond to type I (high levels of SVMPs and low toxicity, LD50 > 1 mg/g mouse body weight) and type II (low metalloproteinase activity and high toxicity, LD50 < 1 mg/g mouse body weight) [66]. In Neotropical rattlesnakes, the adaptive pressure for type I to type II transition was the gain of neurotoxicity and  (Table 1). Distribution ranges were adapted from [3] and The Reptile Database (http://www.reptile-database.org), and are color-coded: green, PLA 2 -rich venom phenotype; red, 3FTx-predominant venom composition. The arrow highlights the trend towards diverging venom phenotypes along the Micrurus north-south dispersal, suggesting the epicenter of the divergence in Mesoamerica.

Tracing the Evolutionary Origin of the 3FTx/PLA 2 Dichotomy across Elapidae
Elapid snakes are a relatively young group. The crown elapid radiation is approximately 38 million years (My) old, yet elapid snakes exhibit some of the highest diversification rates in reptiles [57]. In particular, two clades, Hydrophis and Micrurus, show anomalously high rates of diversification within Elapidae [36]. The Asian coral snake genus Calliophis emerges as monophyletic and the sister group to all other American and Asian coral snakes [47], suggesting an Asiatic origin for the common ancestor of these elapids [58]. Coral snakes diversified around 30-25 My ago, and sea snakes (Hydrophiini) are approximately 16 My old [36]. Estimated divergence times suggest that Hydrophiini is a young and rapidly speciating clade that had a common ancestor~8 My ago, although the majority of extant lineages diversified more recently, over the last~1.5-3.5 My [36,59].
The evolutionary origin and adaptive relevance of the puzzling 3FTx/PLA 2 dichotomy remains elusive since 3FTx-and PLA 2 -predominant venoms are scattered through the phylogenetic tree of Micrurus. The dominant protein families in the venom proteome of C. bivirgata flaviceps are PLA 2 (41%), 3FTx (22.6%) and SVMP (18.7%) [65], suggesting that the ancestor venom phenotype of coral snakes might have been of the PLA 2 -predominant type. Alternatively, 3FTx-and PLA 2 -rich venom proteomes may represent pedomorphic and ontogenetic traits, as has been documented in Crotalus. Rattlesnake venoms belong to one of two distinct phenotypes, which broadly correspond to type I (high levels of SVMPs and low toxicity, LD 50 > 1 mg/g mouse body weight) and type II (low metalloproteinase activity and high toxicity, LD 50 < 1 mg/g mouse body weight) [66]. In Neotropical rattlesnakes, the adaptive pressure for type I to type II transition was the gain of neurotoxicity and lethality to rodents, and this transition represents a miRNA-modulated pedomorphic trait that correlates with the increased concentration of crotoxin along the axis of Crotalus radiation in South America [67,68]. In Nearctic species, such as C. s. scutulatus (Css), the venom phenotype changes geographically from SVMP-rich to Mojave toxin-rich (type-II) as one moves from south central to southeastern Arizona, with a transitional zone between the SVMP and Mojave toxin phenotypes [69]. Understanding the phylogenetic origin of the 3FTx-rich and PLA 2 -predominant venom phenotypes across Micrurus, and whether there is parallelism between the 3FTx/PLA 2 (Micrurus) and the type I/type II (Crotalus) venom dichotomies, requires genus-wide profiling of the venom proteomes and the venom gland transcriptomes of adult and juvenile Crotalus and Micrurus congeneric specimens.

Concluding Remarks and Perspectives
While genomic studies on model organisms have been widely applied to identify traits involved in maintaining functional genetic variation, few studies have drawn the links between genotype, phenotype, and fitness, and the environmental pressures that act to maintain variation that affects organismal phenotypes [70,71]. Venom represents a key adaptive trophic trait that has played an important role in the radiation of advanced snakes, and thus could be a particularly powerful model system to investigate players and mechanisms of adaptive variation at the phenotype level [17,70,72]. In this regard, structural and functional studies of Micrurus venoms point to balancing selection as the mechanism acting to maintain a 3FTx/PLA 2 venom dichotomy. The realization that this 3FTx/PLA 2 dichotomy may have deep roots in the evolution of coral snakes may also have important translational implications. Thus, both types of venoms should be part of immunization mixtures aimed at generating broad-spectrum antivenoms.

Venom
Adult M. t. tschudii specimens were caught in the Peruvian Pacific coast regions of La Libertad, Ancash, and Lima, and kept in captivity at the serpentarium of the Instituto Nacional de Salud, Lima, Perú. Venom from nine adult M. t. tschudii specimens was collected and pooled during the first three months of captivity. Venom was lyophilized and stored at´20˝C until used.

Isolation and Characterization of Venom Proteins
First 0.5 milligrams of crude, lyophilized venom were dissolved in 200 µL of 5% acetonitrile in water containing 0.1% trifluoroacetic acid (TFA), centrifuged to remove debris, and separated by reverse-phase HPLC using a Teknokroma Europa Protein 300 C18 (0.4 cmˆ25 cm, 5 µm particle size, 300 Å pore size) column and an LC 1100 High Pressure Gradient System (Agilent Technologies, Santa Clara, CA, USA) equipped with DAD detector and micro-Auto-sampler [73]. The flow rate was set to 1 mL/min and the column was developed with a linear gradient of 0.1% TFA in water (solution A) and acetonitrile (solution B) using the following column elution conditions: isocratically (5% B) for 5 min, followed by 5%-25% B for 10 min, 25%-45% B for 60 min, and 45%-70% for 10 min. Protein detection was carried out at 215 nm with a reference wavelength of 400 nm. Fractions were collected manually, dried in a vacuum centrifuge (Savant), and redissolved in water, and submitted to molecular mass determination using a SYNAPT ® G2 High Definition Mass Spectrometry System (Waters Corp., Milford, MA, USA), and SDS-PAGE analysis in 15% polyacrylamide gels, under reducing and non-reducing conditions. Gels were stained with Coomassie Brilliant Blue R-250 (Sigma-Aldrich, St. Louis, MO, USA).

Characterization of the Venom Peptidome and Proteome
Electrophoretic protein bands were excised from a Coomassie Brilliant Blue-stained SDS-PAGE gel and subjected to in-gel reduction (10 mM dithiothreitol) and alkylation (50 mM iodoacetamide), followed by overnight sequencing-grade trypsin digestion (66 ng/µL in 25 mM ammonium bicarbonate, 10% acetonitrile; 0.25 µg/sample) in an automated processor (ProGest Protein Digestion Workstation, Genomic Solution Ltd., Cambridgeshire, UK) following the manufacturer's instructions. Tryptic digests were dried in a SpeedVac (Savant™, Thermo Scientific Inc., West Palm Beach, FL, USA), redissolved in 15 µL of 0.1% formic acid in water, and submitted to LC-MS/MS. To this end, tryptic peptides were separated by nano-Acquity UltraPerformance LC ® (UPLC ® , Waters Corporation, Milford, MA, USA) using BEH130 C18 (100 µmˆ100 mm, 1.7 µm particle size) column in-line with a SYNAPT ® G2 High Definition Mass Spectrometry System (Waters Corp., Milford, MA, USA). The flow rate was set to 0.6 µL/min and the column was developed with a linear gradient of 0.1% formic acid in water (solution A) and 0.1% formic acid in acetonitrile (solution B), isocratically 1% B for 1 min, followed by 1%-12% B for 1 min, 12%-40% B for 15 min, 40%-85% B for 2 min. Doubly and triply charged ions were selected for collision-induced dissociation (CID) MS/MS. Fragmentation spectra were interpreted (a) manually (de novo sequencing); (b) using the on-line form of the MASCOT program at http://www.matrixscience.com against NCBInr database, a comprehensive, non-identical protein database compiled from GenBank CDS translations, PIR, SwissProt, PRF, and PDB; and (c) processed in Waters Corporation's (Milford, MA, USA) ProteinLynx Global SERVER 2013 version 2.5.2. (with Expression version 2.0) and the generated .pkl peak list files were exported to MASCOT for protein identification against the NCBInr database. MS/MS mass tolerance was set to˘0.6 Da. Carbamidomethyl cysteine and oxidation of methionine were selected as fixed and variable modifications, respectively. Cut-off for MASCOT reporting was set to top 10 hits All MASCOT identifications were manually verified. Amino acid sequence similarity searches were performed against the NCBInr and UniProtKB databases using the BLASTP program implemented in the WU-BLAST2 search engine at http://www.bork.embl-heidelberg.de.
The relative abundances (expressed as percentage of the total venom proteins) of the different protein families were calculated as the ratio of the sum of the areas of the reverse-phase chromatographic peaks containing proteins from the same family to the total area of venom protein peaks in the reverse-phase chromatogram [14,74]. When more than one protein band was present in a reverse-phase fraction, their proportions were estimated by densitometry of Coomassie-stained SDS-polyacrylamide gels using ImageJ version 1.47 (Free Software Foundation, Boston, MA, USA) (http://rsbweb.nih.gov/ij). Conversely, the relative abundances of different proteins contained in the same SDS-PAGE band were estimated based on the relative ion intensities of the three more abundant peptide ions associated with each protein by MS/MS analysis. Finally, protein family abundances were estimated as the percentages of the total venom proteome.

Determination of LD 50 for Mice
Animal experiments were performed in accordance with protocols approved by the Institutional Committee for the Care and Use of Laboratory Animals of the University of Costa Rica (CICUA 041-15, 21 October 2015), using CD-1 mice of either sex, provided by Instituto Clodomiro Picado. To evaluate the lethal activity, various doses of M. t. tschudii venom, dissolved in 200 µL of PBS, were injected into groups of five CD-1 mice (16-18 g body weight), by the intraperitoneal (i.p.) route. Deaths were scored over a 48 h period and the median lethal dose (LD 50 ) was calculated by probits [75].