Squeezers and Leaf-cutters: Differential Diversification and Degeneration of the Venom System in Toxicoferan Reptiles*

Although it has been established that all toxicoferan squamates share a common venomous ancestor, it has remained unclear whether the maxillary and mandibular venom glands are evolving on separate gene expression trajectories or if they remain under shared genetic control. We show that identical transcripts are simultaneously expressed not only in the mandibular and maxillary glands, but also in the enigmatic snake rictal gland. Toxin molecular frameworks recovered in this study were three-finger toxin (3FTx), CRiSP, crotamine (beta-defensin), cobra venom factor, cystatin, epididymal secretory protein, kunitz, l-amino acid oxidase, lectin, renin aspartate protease, veficolin, and vespryn. We also discovered a novel low-molecular weight disulfide bridged peptide class in pythonid snake glands. In the iguanian lizards, the most highly expressed are potentially antimicrobial in nature (crotamine (beta-defensin) and cystatin), with crotamine (beta-defensin) also the most diverse. However, a number of proteins characterized from anguimorph lizards and caenophidian snakes with hemotoxic or neurotoxic activities were recruited in the common toxicoferan ancestor and remain expressed, albeit in low levels, even in the iguanian lizards. In contrast, the henophidian snakes express 3FTx and lectin toxins as the dominant transcripts. Even in the constricting pythonid and boid snakes, where the glands are predominantly mucous-secreting, low-levels of toxin transcripts can be detected. Venom thus appears to play little role in feeding behavior of most iguanian lizards or the powerful constricting snakes, and the low levels of expression argue against a defensive role. However, clearly the incipient or secondarily atrophied venom systems of these taxa may be a source of novel compounds useful in drug design and discovery.

Venom secretory proteins evolve by a process in which the duplicate of a gene encoding for a key regulatory or bioactive protein is recruited for selective, tissue-specific expression in a venom gland (1,2). Subsequently, the venom gland specific gene isoform is further duplicated, often explosively. Evidence of accelerated evolution and positive selection is characteristic of this "birth and death" model of venom protein evolution (1,(3)(4)(5)(6)(7)(8)(9)(10)(11). For example, the A-superfamily of conotoxin genes isolated from venomous cone snails is among the most rapidly evolving protein-coding gene family identified in metazoans to date (5). In particular, the bioactive residues on the surface of the toxic protein molecules undergo the highest level of accelerated evolution (4,(12)(13)(14)(15). This diversification produces neo-functionalized forms with a myriad of new activities and potencies (1,4,5,7,16).
We have previously demonstrated that a core set of toxinencoding genes was expressed in the oral glands of the common ancestor of all venomous reptiles (Toxicofera clade of squamates). This ancestral set was subsequently augmented by the recruitment of additional toxins and formed the substrate for the evolution of the complex venoms observed in modern venomous lizards and snakes. (Fig. 1A) (2,7,(17)(18)(19)(20)(21)(22)(23)(24). Although Toxicofera, which consists of three lineages (Anguimorpha lizards, Iguania lizards, and Serpentes) is well established as a monophyletic clade within the squamate reptiles (18,25); phylogenetic studies based on nuclear genes have thus far failed to resolve the relative relationships of the snake/iguanian/anguimorph trichotomy [cf (2,25,26). However, other sources of evidence such as SINEs (short interspersed nuclear elements) and morphology favor the clustering of snakes with anguimorph lizards (27,28). We therefore follow this arrangement here (Fig. 1A). The ancestral toxicoferan condition was the possession of thin, simple seromucous glands in both the mandibular and maxillary regions (18). These glands produced active substances that were the substrate for the evolution of toxins. Iguanian lizards apparently split off while this system was only in an incipient stage and little derivation has subsequently occurred within this clade. The apparent lack of specialization of the venom system within the Iguania is likely because of the fact that most species are insectivorous or vegetarian; however, we have shown that toxin genes continue to be expressed by the oral glands of these lizards (18). In the snakes and Anguimorpha lizards, by contrast, the venom system has been extensively diversified (18,20,22,24,29). In the anguimorph lizards, specialization of the mandibular gland was favored, whereas in snakes it was the maxillary gland that became most specialized; in both lineages highly derived forms exist, including encapsulated glands.
Snakes are known for their extensive development of oral glands, which may include a great variety of secretory structures mainly in the upper but also in the lower jaw, some of which develop in association with the dental lamina. Intriguing questions arise in relation to the number, volume, and diversification of these glands: why are there so many glands differing in general architecture, gross anatomy, and histochemistry, and how are these differences reflected in their function? In particular, within the complex of the circumlabial oral glands, those at the corner of the mouth (rictal glands) have received little attention except for a few detailed studies and some general investigations (30 -37). McDowell (38) carried out a detailed study of the corner of the mouth of a large number of caenophidian snakes and described the rictal glands and muscles of the area together with a reinterpretation of the homology of the involved structures. Underwood (36) investigated the circumlabial oral glands of a series of henophidian snakes including Cylindrophis ruffus and the powerful constricting snake Aspidites melanocephalus and paid special attention to the glands of the rictal region. In C. ruffus, he found well-developed maxillary and mandibular glands as well as superior and inferior rictal glands. He described the lobules of the maxillary and mandibular glands as mixed: serous glands leading to mucous tubules that open into the oral cavity via short ducts along the margin of the lips. The C. ruffus superior and inferior rictal glands consisted of branched serous tubules opening into mucous tubules that led into wide mucous ducts opening directly into the corner of the mouth. Phisalix and Caius (39,40) reported that the secretions of the rictal glands of Eryx conicus and six species of uropeltid snake were toxic to birds. These interesting observations do not appear to have been taken any further in the more than 80 years since they were reported. More recently, a test of oral secretions from a wide range of Australian snakes produced a curious cross-reactivity between oral secretions of the pythonid A. melanocephalus and tiger snake venom in the Australian Snake Venom Detection Kit (41). This result may be indicative of the fact that toxins are still being secreted from the oral glands of A. melanocephalus in detectable quantities.
The relative expression of genes in mandibular versus maxillary glands has remained uninvestigated. In addition, the molecular evolution of any associated secreted proteins in henophidian snakes has remained entirely neglected and that of the iguanian lizards has been hampered by only a single species having been investigated to-date and with only limited sequencing (18). Thus, the over-arching aim of this study was to examine the diversification of these glands and their secretory proteins in henophidian snakes and iguanian lizards to elucidate several intriguing aspects of reptile venom evolution. We also examined the genome of Anolis carolinensis for genes encoding putative toxin sequences.
In this study, iguanian lizards and henophidian snakes were both shown to express (at detectable levels) toxin types typically associated with the venoms of caenophidian snakes or anguimorph lizards. Our results revealed that the mandibular and maxillary glands are under a single system of genetic control, with identical transcripts being simultaneously expressed in both glands. Additionally, we showed that the rictal glands of snakes not only also secrete toxins and thus represent an unappreciated source of venoms but that this gland is under the same gene-expression control as the mandibular and maxillary glands. MATERIALS  Histology-Histological sections were prepared from the intact head and the excised venom gland and delivery system. Whole heads were removed and a cut made to the underside to allow fast penetration of the fixative (10% neutral buffered formalin). After a minimum of 2 days excess tissue was removed and specimens were immersed in Kristensen's decalcification solution and placed on a rotor for up to 3 weeks (depending on the size of the head). Before processing the heads were bisected longitudinally for cutting transversely, at 3 microns, in two separate blocks. The processing schedule was: 10% neutral buffered formalin 2 h; absolute ethanol 4 ϫ 1 h; histolene 3 ϫ 1 h; followed by paraffin wax 2 ϫ 90 min. The sections were taken every 100 microns and matching sections were stained with Periodic Acid Schiff's stain and Masson's Trichrome stain. In other specimens, the venom gland, venom duct, and the adjacent bony and muscular tissue were excised and placed in decalcifying solution (Cal-Ex, Fisher) for 72-168 h. Each sample was dehydrated and cleared through a progressive ethanol series and Cyto-Sol (Fisher) before embedding in Paraplast (Fisher). Serial sections were cut at 10 -12 microns. All species were sectioned in the frontal plane; when available, the contralateral venom gland and delivery system were sectioned either parasagitally or transversely. Sections were stained using a variant of Van Gieson's stain, which provides clear distinction between connective tissue, muscle, and epithelium, or with hematoxylin and eosin.
Magnetic Resonance Imaging-Magnetic resonance imaging (MRI) was used to examine the three-dimensional shape and internal anatomy of the venom glands. Formalin-ethanol fixed heads were first submersed in Fomblin (Solvay Solexis) and placed under vacuum to prevent air artifacts. Depending on head size, imaging was performed on either 9.4 T (small/medium) or 17.6 T (large) vertical 89-mm-bore systems (Bruker BioSpin, Rheinstetten, Germany) with a Bruker Mi-cro2.5 gradient system of 1 T/m and transmit/receive birdcage radiofrequency coil with diameter of 10 -30 mm. Bruker ParaVision 4.0 software was used for image acquisition. Anatomical images were acquired using a 3D gradient echo sequence. The field-of-view and matrix were varied to fit the individual samples, resulting in voxel sizes between (40) 3 mm3 and (70) 3 mm 3 . Imaging parameters were: TE ϭ 8 ms, TR ϭ 40 ms, flip angle 20 0 , 4 -8 averages, total scan time 3-9 h per sample, depending on size and resolution. Image segmentation of glands was performed manually in Amira 4.1 (Mercury Computer Systems Inc.) and 3D surface renderings generated.
Transcriptomes-Mandibular venom glands of three specimens were pooled and total RNA extracted using the standard TRIzol Plus method (Invitrogen). Extracts were enriched for mRNA using standard RNeasy mRNA mini kit (Qiagen) protocol. mRNA was reverse transcribed, fragmented and ligated to a unique 10-base multiplex identifier (MID) tag prepared using standard protocols and applied to one PicoTitrePlate (PTP) for simultaneous amplification and sequencing on a Roche 454 GS FLXϩ Titanium platform (Australian Genome Research Facility). Automated grouping and analysis of sample-specific MID reads informatically separated sequences from the other transcriptomes on the plates, which were then post-processed to remove low quality sequences before de novo assembly into contiguous sequences (contigs) using the MIRA software program. Assembled contigs were processed using CLC Main Work Bench (CLC-Bio) and Blast2GO bioinformatic suite (42)(43)(44)(45) to provide Gene Ontology, BLAST, and domain/Interpro annotation. The above analyses assisted in the rationalization of the large numbers of assembled contigs into phylogenetic "groups" for detailed phylogenetic analyses outlined below. Sequence accession numbers are given in Table I and sequences are also available in supplementary File S1, nucleotide sequences.
Phylogenetics-Phylogenetic analyses were performed to allow reconstruction of the molecular evolutionary history of each toxin type for which transcripts were bioinformatically recovered. Toxin sequences were identified by comparison of the translated DNA sequences with previously characterized toxins using a BLAST search (46) of the UniProtKB protein database. Molecular phylogenetic analyses of toxin transcripts were conducted using the translated amino acid sequences. Comparative sequences from other venomous reptiles and physiological gene homologs identified from nonvenom gland transcriptomes were included in each data set as outgroup sequences. To minimize confusion, all sequences obtained in this study are referred to by their Genbank accession numbers (http:// www.ncbi.nlm.nih.gov/sites/entrez?db ϭ Nucleotide) and sequences from previous studies are referred to by their UniProtKB accession numbers (http://www.expasy.org/cgi-bin/sprot-search-ful). Resultant sequence sets were aligned using CLC Mainbench. When presented as sequence alignments, the leader sequence is shown in lowercase and cysteines are highlighted in black. Ͼ and Ͻ indicate incomplete N/5Ј or C/3Ј ends, respectively. Data sets were analyzed using Bayesian inference implemented on MrBayes, version 3.2 (47,48). Run conditions were lset rates ϭ invgamma with prset aamodelpr ϭ mixed. The analysis was performed by running a minimum of 1 ϫ 10 7 generations in four chains, and saving every 100th tree. The loglikelihood score of each saved tree was plotted against the number of generations to establish the point at which the log likelihood scores reached their asymptote, and the posterior probabilities for clades established by constructing a majority-rule consensus tree for all trees generated after completion of the burn-in phase.
Selection Analysis-Sophisticated likelihood models of coding-sequence evolution (53,54) implemented in CODEML of the PAML (55) package were employed to evaluate the selection pressures on the major toxin transcripts recovered from henophidian snake and iguanian serous and mucus glands. We first employed the lineage-specific one-ratio model which assumes a single for the entire phylogenetic tree. This model tends to be very conservative and can only detect positive selection if the ratio averaged over all the sites along the lineage is significantly greater than one. Because such lineage-specific models assume a single for the entire tree, they often fail to identify regions in proteins that might be affected by episodic selection pressures and ultimately, underestimate the strength of selection. Hence, we employed site-specific models which estimate positive selection statistically as a nonsynonymous-to-synonymous nucleotide-substitution rate ratio () significantly greater than 1. We compared likelihood values for three pairs of models with different assumed distributions as no a priori expectation exists for the same: M0 (constant rates across all sites) versus M3 (allows the to vary across sites within "n" discrete categories, nՆ3); M1a (a model of neutral evolution) where all sites are assumed to be either under negative ( Ͻ1) or neutral selection ( ϭ 1) versus M2a (a model of positive selection), which in addition to the site classes mentioned for M1a, assumes a third category of sites; sites with Ͼ1 (positive selection) and M7 (Beta) versus M8 (Beta and ), and models that mirror the evolutionary constraints of M1 and M2 but assume that values are drawn from a beta distribution (56). Only if the alternative models (M3, M2a, and M8: allow sites with Ͼ1) show a better fit in Likelihood Ratio Test (LRT) relative to their null models (M0, M1a, and M8: do not show allow sites Ͼ1), are their results considered significant. LRT is estimated as twice the difference in maximum likelihood values between nested models and compared with the 2 distribution with the appropriate degree of freedom-the difference in the number of parameters between the two models. The Bayes empirical Bayes (BEB) approach (Yang, 2005) was used to identify amino acids under positive selection by calculating the posterior probabilities that a particular amino acid belongs to a given selection class (neutral, conserved or highly variable). Sites with greater posterior probability (PP Ն 95%) of belonging to the ' Ͼ 1 class' were inferred to be positively selected.
Further, Single Likelihood Ancestor Counting (SLAC), Fixed-Effects Likelihood (FEL), and Random Effects Likelihood models (57) implemented in HyPhy (58) were employed to further provide significant support to the aforementioned analyses and to detect sites evolving under the influence of positive and negative selection. The more advanced Mixed Effects Model Evolution (MEME) (59) was also used to detect episodic diversifying selection. MEME employs FEL along the sites and Random-effects likelihood (REL) across the branches to detect episodic diversifying selection.
For the effective comparison of iguanian and crotaline snake crotamines (beta-defensins); iguanian and homalopsid veficolins and anguimorph lizard and caenophidian snake cysteine-rich secretory proteins (CRiSPs) 1 , we employed Codeml Clade Model C analysis (60). The significance of the analysis was tested by comparison with site model M1a.
Structural Analyses-To depict the selection pressures shaping the evolution of toxins, we mapped the sites under positive selection on the homology models created using Phyre 2 webserver (61). Pymol 1.3 (62) was used to visualize and generate the images of homology models. Consurf webserver (63) was used for mapping the evolutionary selection pressures on the three-dimensional homology models whereas Selecton webserver (64,65) was used to identify hypervariable regions on the toxins. Three-dimensional models and surface charge distributions were constructed as per (19).
Proteomics-Aspidites melanocephalus were obtained from Port Headland (Western Australia), the same population used as a "nega-tive control" by Jelinek et al. (41). Specimens were anesthetized using zoletil (3 mg/kg) and then injected subcutaneously with pilocarpine (20 mg/kg) periglandular to the labial glands. Polyethylene equipment was used to collect and process samples in all cases. Samples were subsequently filtered using 20 Å syringe filters to remove large mucoidal strands and then lyophilized. Reduction and alkylation were undertaken by redissolving 3 g of sample from each of the localities in 50 l of 100 mM ammonium carbonate. After sonication and vortexing, 50 l of 2% iodoethanol/0.5% triethylphosphine in acetonitrile was added to the re-dissolved sample. The reaction mixture was incubated for 2 h at 37°C, before being dried and re-suspended in 15 l of 2.5% acetonitrile and 0.5% formic acid. Additionally, 3 l of reduced and alkylated sample was resuspended in 30 l of 40 mM ammonium bicarbonate, before being incubated overnight with 750 ng trypsin. Digestion was stopped by addition of 1 l of concentrated formic acid and 0.75 g processed by liquid chromatography-tandem MS (LC-MS/MS). LC-MS/MS was performed using a Agilent Zorbax stable bond C18 column (2.1 mm ϫ 100 mm, 1.8 m particle size, 300 Å pore size) at a flow of 400 l/min and a gradient of 1-40% solvent B (90% acetonitrile, 0.1% formic acid) in 0.1% formic acid over 15 min on a Shimadzu Nexera UHPLC coupled with an AB SCIEX 5600 Triple TOF mass spectrometer equipped with a Turbo V ion source heated to 450°C. MS 2 spectra were acquired at a rate of 20 scans/s, with a cycle time of 2.3 s, and optimized for high resolution. Precursor ions were selected between 80 -1800 m/z, charge state of 2-5, and of an intensity of at least 120 counts per second, with a precursor selection window of 1.5 Da, and excluding isotopes within 2 Da for MS 2 . MS 2 spectra were searched against all-reading frame translated A. melanocephalus mandibular and maxillary gland transcriptomes (16216 amino acid sequences) with Proteinpilot v4.5 (ABSciex) using a thorough identification search, specifying alkylation method (iodoethanol), tryptic digestion, and allowing for biological modifications and amino acid substitutions. This was performed to maximize the identification of protein sequences from the transcriptome despite the inherent variability of toxins, potential isoform mismatch with the transcriptomic data, and account for experimental artifacts leading to chemical modifications. For specific search parameters, modification probability and amino acid substitution probability settings set for the Paragon algorithm (version 4.5.0.0, 1654) see supplementary File S2. Spectra were inspected manually to eliminate false positives, excluding spectra with low signal to noise, erroneous modification assignments, and confidence values below 95% unless justifiable by presence of a-ions. The results were further analyzed using CLC Main Workbench v6.6.
1 Dimensional Gel Electrophoresis-SDS-PAGE of 20 g A. melanocephalus secretion was carried out on 12% polyacrylamide gels as per the Laemmli-method, using a Bio-Rad Mini Protean setup. All reagents were purchased from Bio-Rad. Gels were stained overnight in Coommasie brilliant blue staining solution and de-stained again with 10% acetic acid.

RESULTS
Gland Arrangements-The rictal glands of the henophidian snake C. ruffus were large enough that they indented the mandibular and maxillary glands (Fig. 1B). In this anatomically simple and presumably plesiotypic condition, entirely serous tubules merge (Fig. 1C) to join a single duct descending vertically slightly anterior to the corner of the mouth. C. ruffus maxillary and mandibular glands were mixed glands, with the mandibular gland containing a greater proportion of serous cells than did the maxillary. The mandibular and maxillary glands of the constricting boid and pythonid snakes (Aspid- 1 The abbreviations used are: CRiSPs, cysteine-rich secretory proteins; CVF, cobra venom factor; ESP, epididymal secretory protein; RAP, renin aspartate protease; ORM, one ratio model; MYA, million years ago. ites, Boa, Corallus, Eunectes, and Morelia) were also mixed but were predominantly composed of mucoid cells with only scattered serous cells (Fig. 1D). The iguanian lizards had mixed sero-mucous maxillary and mandibular glands, with small insectivorous or vegetarian species of any size having primarily mucous cells while species that feed on vertebrates, such as Anolis equestris, contained a greater proportion of serous cells at the base of the gland arrangements (Fig. 1E). In contrast, all nontoxicoferan lizards examined in this study (Ameiva ameiva, Cordylus tropidosternum, Eumeces obsoletus, Gekko gecko*, and Lepidophyma flavimaculatum) had mandibular mucous glands (e.g., Lepidophyma flavimaculatum in Fig. 1F) and had only extremely small maxillary glands or lacked them entirely.
Transcriptome Sequencing-Sequencing of transcriptomes of iguanian lizard and henophidian snake glands recovered transcripts encoding for proteins homologous to toxins previously characterized from anguimorph lizards or caenophidian snakes: 3FTx, CRiSP, crotamine (beta-defensin), CVF (cobra venom factor), cystatin, ESP (epididymal secretory protein), kunitz peptide, L-amino oxidase, lectin, RAP (renin aspartate protease), veficolin, and vespryn (Table I). Phylogenetic analyses showed all toxicoferan sequences to be monophyletic relative to endogenous body sequences, strongly suggestive of a single early recruitment, and in most cases the lizard and snake sequences were not reciprocally monophyletic, thus conclusively demonstrating a single early origin for that protein type. Intriguingly, whenever a particular toxin type was recovered, identical toxin transcripts were present in mandibular and maxillary glands within a species and in the case of C. ruffus, the identical transcripts were also recovered from the rictal gland in addition to the aforementioned glands.
3FTx recovered from the henophidian snakes form a clade with the reptile and bird specific toxins (19,66,67) from nonfront fanged snakes along with the sequences from the viperid snake Sistrurus catenatus Fig. 2. The C. ruffus and A. melanocephalus isoforms shared a unique, newly evolved cysteine immediately before the doublet formed by the eighth and ninth ancestral cysteine. This results not only in a novel cysteine pattern but also in an odd number of cysteines (eleven). These changes would have structural implications not only for disulfide-bridging but also in the potential facilitation of dimerization. One C. ruffus variant has an exon deletion removing a significant N-terminal stretch of residues including five cysteines, resulting in a novel short peptide containing six cysteines (Fig. 3).
The iguanian crotamine (beta-defensin) sequences were also the most sequentially and phylogenetically divergent of all the iguanian transcript types (Figs. 4, 5). Some iguanian lizard crotamine (beta-defensin) sequences and the anguimorph lizard sequences unique proline-rich, highly charged N-terminal extensions not seen in the snake venom forms (Fig. 4). Sequence alignment and lack of phylogenetic affinity indicates that these novel proline rich regions are convergently evolved. A derived form from V. glauerti was notable in lacking two of the ancestral cysteines (Fig. 5).
All lizard kunitz sequences had the same two-domain structure as the sequence from the venom of the caenophidian

Toxin type Species recovered from (Genbank accession numbers) 3FTx
Aspidites snake Austrelaps labialis (Fig. 6). Phylogenetic analyses showed that the two domain forms were basal to the derived mono-domain form (Fig. 7).
The henophidian snake lectin toxins had the galactosebinding EPN tripeptide functional motif whereas the C. calyptratus form is of the lactose-binding QPD motif. Phylogenetic analysis showed that the sequences containing the functional motifs EPN (glutamic acid ϩ proline ϩ asparagine) and QPD (glutamine ϩ proline ϩ aspartic acid) were not reciprocally monophyletic (Fig. 8).
Comparing the values estimated by the above models can be misleading, because the genes being compared may have a different proportion of sites under selection. Hence, we employed the Clade Model C approach for effectively comparing the evolution of anguimorph lizard, iguanian lizard, and crotaline snake crotamines; iguanian lizard and homalopsid snake veficolins and anguimorph lizard and caenophidian snake CRiSPs (Table III). Clade model analyses indicated that the crotaline snake crotamines ( ϭ 4.43) and homalopsid snake veficolins ( ϭ 3.06) have accumulated greater variations than their iguanian lizard crotamine (beta-defensin) ( ϭ 3.61) and veficolin ( ϭ 0.33) counterparts. The varanid crotamine (beta-defensin) peptides were less diverse ( ϭ 1.97) than either the snake or iguanian lizard counterparts. It is important to note that a lower omega value is only expected among closely related taxa versus more diverged taxa if for example, they share more ecological similarity to the extent that their toxins did not need to selectively diversify as much as between ecologically diverse taxa. The analysis also revealed that the caenophidian snake CRiSPs ( ϭ 3.94) have evolved more rapidly than the anguimorph lizard CRiSPs ( ϭ 2.04). All these aforementioned analyses were significant (PϽ Ͻ0.001) in comparison with site model M1a. Single Likelihood Ancestral Counting (SLAC), Fixed Effects Likelihood (FEL), Random Effects Likelihood (REL), Mixed Effects Model of Evolution (MEME) and the identification of hypervariable regions, conclusively support these results (Tables II, III,

Molecular & Cellular Proteomics 12.7
garus Candidus venom as a template (PDB-code: 2JQP). Alternative templates (1ABT, 1CDT, 1KBA, and 2H7Z) were investigated, because not all of the sequences yielded a suitable similarity score to a single template structure. The bungatoxin structure used in this study was selected for its best fit to the cysteine sequence pattern. Regions of insertions and deletions were moved to loop structures rather than modeled in the core of the protein. Target template alignments were generated using SPDBV (68). Model structures were obtained from SwissModel server using the "project mode" option (68).
Proteomics-Proteomic analysis of the A. melanocephalus secretion was able to confirm the presence of the two dominant transcript types: 3FTx and lectin. Tandem MS spectra of tryptic fragments provided high confidence matches for each of the two toxin sequences (Fig. 10) (Supplemental Table  S12). The low number of identified proteins, despite the large amount of analyzed sample (1 g), supports the histological findings that the maxillary and mandibular glands of A. mela- nocephalus are primarily mucoidal but that low levels of protein secretion continue. The 3FTx present in the secretion was present in the identical Ͻ10kDa region in 1D gel analysis under reducing and non-reducing conditions (Fig. 11), demonstrating that the free cysteine (Fig. 3) is buried and thus does not promote dimerization. DISCUSSION Toxin gene recruitment events appear to have occurred in batches on three occasions: up to ten recruitments at the base of the toxicoferan tree; six in the common ancestor of anguimorph lizards and snakes; and eight more in the common ancestor of the caenophidian snakes (Fig. 1). Thus the Toxicofera venom system appears to have a punctuated molecular evolutionary history.
Our results significantly changed our perception of the relative timing of gene recruitment events for several toxin types (Fig. 1). 3FTx recruitment, previously thought to have occurred within the caenophidian snakes before the Viperidae split-off 54 million years ago (MYA), now appears to have occurred within the ancient Afrophidian snake clade at 104 MYA, a backwards shift of 50 million years (71). RAP was previously only known from Echis venoms within the Viperidae snakes and thus was considered to be the result of a recruitment event some 20 MYA. As we recovered it from B. constrictor, the recruitment may also have occurred within the ancestral Afrophidian snakes, a backwards shift of 84 million years. Intriguingly, a homologous renin aspartate protease sequence from the A. carolinensis genome was nested within the snake sequences, strongly suggestive of recruitment at the base of the toxicoferan reptile radiation some 166 MYA, which would be a shift of 146 million years from Echis. However, sequencing of homologous iguanian lizard mandibular or maxillary gland mRNA transcripts or proteins must first be accomplished before this timing can be absolutely confirmed. Other toxin types, however, were conclusively shown to be ancestral to the toxicoferan reptiles: CRiSP, CVF, kunitz, Lamino oxidase, lectin, and veficolin. CRiSP, CVF, lectin, and veficolin were previously sequenced from both snakes and anguimorph lizards but not from iguanian lizards, and thus this represents a backwards shift of 22 million years for the recruitment of this toxin type. Kunitz peptides and L-amino oxidase were considered exclusive to caenophidian snakes and thus our perception of the age of these toxin types shifts over 112 MYA.
This ancestral molecular diversity is particularly intriguing considering that the iguanian lizards have only poorly developed "incipient" venom glands. Despite these glands being of no obvious use in the feeding ecology of insectivorous or vegetarian species, and being predominantly mucoid in these lineages, they still secrete active peptides. Intriguingly, the highest expression level and greatest evidence of duplication and diversification within the iguanian transcripts was ob- served for the crotamine (beta-defensin) and cystatin peptides. All other transcript types were recovered in lower levels and typically with only one or two isoforms. It is noteworthy that the ancestral activity of these two peptide types is antimicrobial (72, 73)-an activity which is retained by at least some isoforms expressed in the venoms of highly toxic frontfanged snakes. Another ancestrally recruited protein type, L-amino oxidase, also has well-documented anti-microbial activity (74). Thus, the trigger for the development of the toxicoferan venom system may have been the rapid co-option of a newly developed antimicrobial peptide secreting system. However, it must be highlighted that the other toxin classes ancestrally recruited and still expressed even in the iguanian venom system (at very low levels) include those that have been developed into neurotoxins (CRiSP (75)(76)(77)(78)(79), vespryn (80) or hemotoxins: kunitz, lectin, and likely veficolin (81)). Thus, in the Iguania the evidence favors an antimicrobial function for the secretions of their incipient venom system. There is no known predatory function for the system within vegetarian or insectivorous species, although the significantly enlarged protein secreting region of species that predate on vertebrates (Fig. 1E) does point to a possible role in prey-capture. Conversely, the rapid accumulation of novel secretory frameworks appears to have underpinned the rapid radiation of the other toxicoferan reptile lineages (anguimorph lizards and the snakes). In particular, the venom system likely led to the explosive diversification of the caenophidian snakes (82).
The well-developed rictal glands in the henophidian snake C. ruffus (Fig. 1B), which are entirely serous in composition (Fig. 1C), suggest an active role for these structures in the life of these serpents. Similarly, the mandibular and maxillary glands of Cylindrophis are sero-mucous whereas those of the powerful constricting highly derived boid and pythonid snakes, stain histologically as predominantly (but not exclusively) mucoidal. Of particular note is the encoding of the same toxin transcripts in maxillary, mandibular, and rictal glands. This supports not only the hypothesized shared evolutionary history of the mandibular and maxillary glands (18) 9. Molecular evolution of lectin, 3FTx, epididymal secretory protein (A, B, and C, respectively), iguanian Cystatin (D), iguanian lizard and crotaline snake crotamines (beta-defensins) (E and F, respectively), iguanian lizard and Homalopsidae snake veficolins (G  and H, respectively). The number of positively-selected sites (Codeml: Model 8, PP Ն 0.95, Bayes-Empirical Bayes approach) and hypervariable regions are highlighted in red and orange, respectively. Various omega statistics (M8 omega value, total number of positively selected sites detected by M8 (PP Ն 0.95, BEB)) and the number of positively and negatively selected sites detected by HyPhy integrative approach (SLAC and FEL: 0.05 significance; REL 50 Bayes factor) are also indicated. but also provides the first evidence that the rictal glands are derived from the same system.
The difference in gland structures and types between C. ruffus and the boid/pythonid snakes parallels the differences in their predatory ecology: with C. ruffus being representative of the presumably ancestral snake condition in being only weak constrictors and feeding upon prey that is typically not wider than their head; while boid/pythonid snakes represent a derived clade that has developed powerful constriction as a new adaptation to predate upon massive prey items often much wider than their head. Thus, alternate selection pressures operate on the oral glands of these snakes: C. ruffus would benefit from chemical assistance in prey capture; whereas the boid/pythonid snakes easily dispatch their prey using powerful constriction but have a need for appreciable amounts of mucous to lubricate the fur of mammalian or feathers of avian prey items. However, consistent with positive reaction generated by A. melanocephalus oral secretion in snake venom detection kit trials (41), our data demonstrates that the toxin genes are not only still transcribed but also continue to be translated and secreted as relics of the venom system that the boids/pythonids have down-regulated. Although it is noteworthy that A. melanocephalus is unusual among pythonids in that it feeds primarily on reptiles (slender prey items not dissimilar to those which C. ruffus preys on), it is probable that their oral secretions are of little to no use in prey capture as the snakes rely on constriction to kill the prey items. Despite this, however, they may be problematic for human bite victims as a result of their capacity to produce false-positives in the medical diagnosis of snake envenomation.
Our sequencing of transcripts that are homologous with well-characterized toxin types also provides new evidence for the continued diversification of the secreted proteins even within the glands of henophidian snakes and iguanian lizards. One example is the cysteine changes in the 3FTx of henophidian snakes that could have profound structural implications not only on the conformation of the peptide but also the facilitation of dimerization. However, our demonstration that the 3FTx in the A. melanocephalus secretion occupy the same Ͻ10kDa location under reducing and nonreducing conditions, demonstrates that the free cysteine (located as part of the cysteine-triplet near the C terminus) is buried in the molecular core and thus does not facilitate dimerization. Our results also provide evidence for ancient derivations such as the early loss of cysteine four in the lectins with replacement by glycine in one of the early duplicates. It had previously been unclear whether the EPN motif (which blocks platelets through mannose-binding sites), or the QPD motif (which blocks platelets through galactose-binding sites) was plesiotypic as the caenophidian snakes and anguimorph lizards each have both forms (20,22). In this study, the QPD motif was present in the sole sequence recovered from an iguanian lizard (Chamaeleo calyptratus). Thus the QPD motif appears to be plesiotypic whereas the EPN is an ancient derivation. However, as neither motif is monophyletic, a degree of convergence in one or the other or both cannot be ruled out. Some henophidian snake isoforms have a glycine substituted for the fourth cysteine, a derivation shared with some caenophidian snake homologs but not present in the lizard forms. This indicates a secondary loss of this cysteine early in snake evolutionary history. In the kunitz peptides, the ancestral dual domain form previously sequenced from the caenophidian snake Austrelaps labialis was also recovered from the iguanian lizard C. calyptratus as well as the anguimorph lizards Abronia graminea and Heloderma suspectum cinctum (83). However, the derived monodomain form favored in caenophidian snake venoms was also recovered from iguanian lizards (C. calyptratus and C. kingii). Thus the loss of the first domain was a derivation that occurred before the organismal diversification of the toxicoferan ancestor.
Evidence provided by various analyses of different henophidian snake and iguanian peptides/proteins and mapping of mutations on the homology models indicates that several of these have accumulated variation under the strong influence of positive selection (Tables II, III, and Fig. 9). The henophidian snake lectins in particular seem to have evolved rapidly ( ϭ 1.51; 40 positively selected (25%) sites). This estimated variation combined with the putative hemotoxic activity revealed by the conservation of functional motifs, may indicate a functional role for orally secreted lectins in the lives of these snakes, which are not known for using venom in prey capture. These analyses also revealed hypervariable regions in iguanian lizard crotamine (beta-defensin) and veficolin toxins. The observed variation in several regions could be a result of the putative antimicrobial function of crotamine (beta-defensin) peptides.
The recovery of novel protein scaffolds from the glands studied here reinforces how little is known about the protein composition of reptile venoms or the evolutionary history of the glands. This is underscored by the number and diversity of novel scaffolds recovered despite the relatively limited sam-pling employed, including an entirely new peptide class discovered in the pythonid libraries. The toxin types recovered in this study should not be considered as the full suite because the relatively limited sampling only recovered the dominant forms. More extensive sampling will no doubt recover novel isoforms of types identified to date as well as entirely new toxin classes. Of interest for follow-up investigations should be the relative secretion by the segregated mucous lobules in the independently derived mandibular systems of varanid and helodermatid lizards as well as the derived maxillary system of snakes. Both of these systems also have mucous and protein secreting cells divided into discrete glands, with the protein glands encapsulated. In the case of the anguimorph lizards, the varanids and helodermatids have independently derived the mandibular glands into complex compartmentalized glands. Further work is also needed to examine the relationship between transcription and translation, and thus the presence and representation of the proteins themselves in the secretory mixtures. These results highlight the relatively untapped potential of such complex mixtures as sources of novel investigational ligands or lead compounds for use in drug design and development. Of particular focus of follow-up work by us will be to investigate the bioactivity of iguanian and anguimorphan crotamine (beta-defensin) peptides to determine relative degree of antimicrobial activity versus more derived activities such as neurotoxicity. It is hoped that these results will stimulate further investigation of these neglected glands and their secretory proteins in an increasingly diverse range of toxicoferan reptile species.