Enter the Dragon: The Dynamic and Multifunctional Evolution of Anguimorpha Lizard Venoms

While snake venoms have been the subject of intense study, comparatively little work has been done on lizard venoms. In this study, we have examined the structural and functional diversification of anguimorph lizard venoms and associated toxins, and related these results to dentition and predatory ecology. Venom composition was shown to be highly variable across the 20 species of Heloderma, Lanthanotus, and Varanus included in our study. While kallikrein enzymes were ubiquitous, they were also a particularly multifunctional toxin type, with differential activities on enzyme substrates and also ability to degrade alpha or beta chains of fibrinogen that reflects structural variability. Examination of other toxin types also revealed similar variability in their presence and activity levels. The high level of venom chemistry variation in varanid lizards compared to that of helodermatid lizards suggests that venom may be subject to different selection pressures in these two families. These results not only contribute to our understanding of venom evolution but also reveal anguimorph lizard venoms to be rich sources of novel bioactive molecules with potential as drug design and development lead compounds.

Further, Sweet's assertion that the genetic support for Toxicofera as a clade is erroneous due to short-branch/long-branch attractions [10], is a flawed argument as long-branch attraction is not a serious issue for recent likelihood-based methods (as it was with older parsimony based methods), particularly when good model-selection procedures are used. Also, for DNA sequence data, long branch attraction is not expected to be an issue where more than 20% of the characters favour that topology (as is the case in the Toxicofera situation [11][12][13][14][15][16][17][18][19][20][21][22]). Indeed Sweet himself notes [10] that only a minority of the evidence cited in a review by Sites, Reeder & Wiens [37] did not support the grouping of anguimorphan and iguanian lizards as part of the Toxicofera clade. Thus Sweet's main argument against Toxicofera is one of inertia due to the weight of history, as demonstrated by his statement "that dozens to hundreds of morphological, behavioral and ecological synapomorphies must be reversed is strong evidence of a false signal in the genes supporting the short Toxicoferan branch" [10]. Naturally the historical morphological, behavioural and ecological observations were viewed through the filter of the understanding of the organismal relationships at the time. The fact that these observations will have to be reinterpreted due to genetic evidence overturning morphology-driven taxonomy, should not be an impediment to the acceptance of newer, more robust data regarding the taxonomical arrangements. Just as the historical ecology and predatory observations of varanid lizards were viewed through the filter that these lizards were assumed to be non-venomous due to the perceived lack of venom glands, need to be reinterpreted (and even redone) in light of the discovery that varanid lizards in fact posses complex oral glands homologous to that of helodermatid lizards and secrete many of the same proteins with similar bioactivity [23][24][25][26][27][28][29][30][31][32][33][34]. Therefore, the weight of history should not stand in the way of scientific advancement. This lack of acceptance of more lizards being venomous than previously understood, is mirrored by the reluctance of some [35,38] to accept that all advanced snakes share a common venomous ancestor [25,29,34,36,39].
Importantly, more recent phylogenetic analyses that have combined morphological and molecular data, and simultaneously attempted to resolve the topological conflict by careful analysis and treatment of its source in the data, have again recovered Toxicofera as a robustly supported clade and uncovered new morphological synapomorphies of the clade [14,21]. Controversy regarding the phylogenetic validity of Toxicofera is often intertwined with that surrounding hypotheses concerning the evolution of venom within the clade, a situation compounded by use of the poorly defined term "Toxicofera hypothesis". Although Hargreaves et al. (HEA [40]) used this term to refer to the "single early evolution of venom in reptiles", the original publication concerning this "early origin" hypothesis almost exclusively prefers the term "venom system" over "venom" [23]-thus the original "Toxicofera hypothesis (of venom evolution) concerns the single early origin of the venom system, not "venom" per se. In that all lineages possessed uniquely derived mandibular and maxillary glands distinguished by having segregated protein and mucus secreting regions, with the enlargement of the protein-secreting region relative to that of the mucus-secreting region. More broadly, however, the term "Toxicofera hypothesis" could refer either to a hypothesis concerning the phylogenetic existence and constituency of the Toxicofera clade, or one of several interpretations of the evolutionary history of venom systems within the clade. Regardless of evidence or views on venom evolution within the clade, the weight of evidence (whether morphological or molecular) strongly supports the existence of the clade Toxicofera. It should be noted that the Hargreaves study [40] challenging the evolution of venom in lizards relied upon unreliable data such as tissues expression values up to 5000-fold in variance within a single set of replicates through to averaging only n = 2 including some in which one of the values was 0 (Supplementary Tables S5-S9 of [40]). Further to this the trees presented in the paper [40] Figures 5 and 6, and Supplementary Figures S1, S2, S4a, S7, S9, S12, S14 and S16-S18 contained nodes Toxins 2017, 9, 242 4 of 37 of less than 50% support or have the inclusion of nodes without support values. Had node support been presented in the main paper with the standard protocol of collapsing nodes below 50%, their trees would be largely polytomies with very little topological information. Given that many of their conclusions were based on the topology of those trees, we must consider that those conclusions were not supported by their data or analyses. Therefore their assertion that many reptile venom proteins are simply expression of 'housekeeping genes', and therefore calling into question the single early origin of the venom system, is not phylogenetically supported.
Some of the controversy surrounding the "Toxicofera hypothesis" (sensu HEA) of venom evolution has been blamed on inconsistencies in definitions of the word "venom" and the designation of certain species as "venomous". In reality, the various definitions of "venom" differ little-the consensus is that venom is an actively delivered (e.g., via a bite or a string) secretion that functions (i.e., has been selected for the purpose of) in the subjugation of prey or the deterrence of predators/competitors [29,38,41]. A popular wording of this definition [42], sometimes considered to be more restricted than more recent formulations [43], introduces a point of confusion by including "subjugation and/or digestion and/or in defence", which suggests that digestion is a distinct function of "venom". This would greatly complicate matters, as oral secretions (e.g., human saliva) in general aid digestion. Here, we prefer the more restricted definition-whilst digestion is undoubtedly a function of the oral secretions of the lizards examined in the present study, that fact alone does not qualify the secretions as "venoms".
Confusion has thus been fueled less by differing definitions of "venom", than by different applications of the term, and particularly its conflation with the term "venom system". This confusion stems from the use of the two terms interchangeably in papers concerning Toxicofera going back at least as far as 2006 [23]. It must be stressed that taxa possessing a "venom system" are not necessarily "venomous"-in particular we have frequently referred to the venom system of iguanian lizards as "incipient" (e.g., [30] and see [41] for a discussion of the use of the term incipient in evolutionary contexts). So, the "Toxicofera hypothesis of venom evolution", more properly considered (although this discussion does not constitute a formulation of such a hypothesis), concerns the early evolution of the venom system-i.e., the synapomorphy of the Toxicofera clade, which (as above) is the presence of protein-secreting oral glands (i.e., an incipient venom system) which may be considered exapted for the subsequent development of sophisticated venom delivery systems.
In making a selectable contribution to the subjugation of prey, it is not necessary for a venom to be capable of rapidly killing or incapacitating the prey item. Venoms do not necessarily function as stand-alone prey subjugation mechanisms-they may be used in concert with physical means of subjugation. Thus, a venom system need only marginally increase a predator's chances of successfully securing a meal, e.g., by slightly weakening a potential prey animal, making it easier to subdue physically [41]. When the effects of a varanid lizard venom are contrasted with those of the venom of a highly front-fanged snake such as a taipan as an argument against such lizards having a 'true' venom [43], the differences are indeed stark but this is a spurious comparison. This does not constitute a strong argument against the lizard being venomous, instead it only serves to illustrate the fact that venom systems exist in myriad states of development within Toxicofera, as they do within the advanced snakes themselves [39].
Criticisms of (some variants of) the Toxicofera hypothesis of venom evolution notwithstanding, it is clear that a core set of toxins are present in the venoms of all anguimorph lizards studied to date including CRiSP, kallikrein, B-type natriuretic peptides and type III phospholipase A 2 (PLA 2 ) [23][24][25][27][28][29][30][31][32][33][34]. It is an important point that the presence of a particular species of protein within a secretion does not demonstrate the function of that secretion, which is why the evolution of venom must be studied by considering the venom system as a whole-functional assays should be deployed, as well as considerations of associated anatomy (i.e., delivery mechanisms) and organismal feeding ecology. Some of the previously listed toxins are responsible for the hypotensive effect of intravenous injection of crude venom in rats [22,32,33], such as aortic smooth muscle relaxation by natriuretic peptides equipotent to forms recovered from venomous snakes [44]. In addition, kinin Toxins 2017, 9, 242 5 of 37 release from kininogen by kallikrein enzymes is another source of hypotensive effects resulting from lizard venoms [45]. In addition, it has also been previously noted that injection of rodents and birds with V. griseus venoms resulted in paralysis [46]. Reports of human envenomation by monitor lizards have been noted both in the literature and anecdotally. A bite report documented effects for a V. griseus bite to a zoo keeper that included symptoms similar to that of the bird and rodent studies: dizziness, muscular weakness and soreness in the extremities, facial and eye muscle pain, respiratory distress and pain (local and systemic) [47]. Additional reports of V. griseus bites to zoo keepers recorded similar symptoms including dysphagia, chest tightness, muscle soreness in the extremities, facial pain, dizziness, and difficulty walking [48]. In all the above V. griseus cases, the bite victims were experienced zoo keepers who did not view the bites with concern upon their occurrence. Vikrant and Verma (2014) reported a lethal bite by the related V. bengalensis that induced local pain, blood loss, as well as nausea, diaphoresis, dizziness, and breathlessness in the victim and eventually led to an acute kidney injury and cardiac arrest [49]. The actual culprit responsible for this bite has been questioned by clinical toxinologists, who considered it more likely to be a venomous snakebite [50]. However, it must be noted that this dissenting opinion by people who were not involved in the case did not provide any new and contradictory evidence. In contrast, a local wildlife officer who was present at the event confirmed the identity of the monitor lizard and the attending physicians documented the bite wounds as being inconsistent with the puncture wounds characteristic of bites by venomous snakes but rather consistent with the lacerations produced by monitor lizard bites (Vikrant personal communication). Thus, while this case is extraordinary, the possibility of dangerous bites from varanid lizards to humans, under exceptional circumstances, should not be discounted. Another case report described a bite by a juvenile V. komodoensis that led to faintness, prolonged bleeding and transient hypotension [51]. The authors attributed these effects to a vasovagal reaction despite noting themselves the similarity to reported in vitro effects of V. komodoensis venom [33] and with vasovagal reaction unable to explain the prolonged bleeding.
Anecdotally, a great many varanid lizard bites to biologists, zookeepers and amateur reptile enthusiasts have resulted in little that could be attributed to the action of toxins; however, some bite victims do report burning sensations, prolonged bleeding and inflammation disproportionate to the mechanical damage inflicted [10]. In addition to published cases, BGF has been contacted by keepers of varanid species who had symptomatic cases, sometimes with taxon-specific effects: 3 cases of V. albigularis reporting extreme muscle soreness and weakness lasting days, 4 cases of V. kordensis producing significant local swelling of the bitten finger and adjacent fingers, 5 cases of V. varius with pronounced bleeding lasting 3-4 h, a V. salvadorii bite with similar symptoms to that of V. varius and a V. komodoensis bite also producing apparent anticoagulant effects. These effects are consistent with the characterised venom chemistry of anguimorph lizards [23,29,32,33].
When compared with helodermatid lizard bites, which have received far more research attention than those of their varanoid cousins, a key point becomes evident: duration of contact. Typically, a helodermatid lizard will stay attached, chewing more venom into the bite site, while a monitor lizard is less likely to hold onto something that is not a food item but in some cases may hold on for up to a half hour when defensively biting. This likely leads to differences in the amount of oral fluids inoculated to the victim, with symptomatic human bites typically being those in which the varanid lizard chewed for a prolonged period of time. Interestingly, while feeding on large prey items, varanid lizards seem to have a tendency to shake it violently and chew vigorously until it is subdued, which may facilitate venom delivery as well as potentiating mechanical damage [52].
Varanoid lizards are characterised by a refined mandibular venom gland that is homologous with that of the helodermatid lizards [23,25,27,29,32,33,53]. Most anguimorph lizards have simple-structured mandibular venom glands, however, Heloderma and Lanthanothus/Varanus have independently evolved complex glands [32]. Both lineages derived their compartmentalised venom glands from the ancestral anguimorph lizard condition of an enlarged, mixed sero-mucous gland in which the protein-and mucous-secreting regions are not segregated into distinct glandular structures. In both cases, sophisticated structures with separated protein and mucous regions, a structured lumen for storing liquid venom, and a thick membranous cover, have evolved. Further to this, morphological as well as molecular evolutionary studies have indicated that these glands are homologous with the venom glands of snakes [23,25,27,29,32,33,53]. The fact that varanid lizards possess highly developed dental glands suggests that those glands in one way or another play an important role in their life. It has been shown with venomous snakes which switch to eggs or constriction as an alternate form of prey capture rapidly lose the functionality of their venom glands, with atrophying occurring in short periods of (evolutionary) time [25,29,54,55].
The evolution of a complex venom system is likely only possible under certain contingent circumstances-i.e., when both environmental conditions select for it and a species' overall evolutionary trajectory facilitates it. For example, in Iguania the incipient venom glands never developed any significant complexity, probably due to the mainly insectivorous/herbivorous nature of these lizards. In addition, when a species evolves an alternative method of subduing prey that renders the venom system redundant, or switches to a diet with no need for subjugation (e.g., plants), the venom system often degrades [25,29,54,55]. The cost of venom production is presumably high enough to justify the presence of active secretory and delivery apparatuses only when it confers an evolutionary advantage [56]. It is notable, however, that in iguanian species which include a large quantity of vertebrates in their diet, the glands are larger and the protein-secreting region more developed (though we stress that this is not equivalent to them being 'venomous' per se but it is strongly suggestive of a functional role in predation) [29,30].
Previous work on helodermatid lizards has demonstrated a striking level of proteomic conservation within the venoms of this clade, with the same basic toxin groups present in similar quantities [27,57] despite the most recent common ancestor (MRCA) of the five Heloderma species existing~15-20 million years ago [11]. While the venom glands of Varanus species have been compared transcriptomically, the only proteomic comparisons to-date were limited to SELDI mass spectrometry [23,32,33]. However a diversity of components have been functionally characterised from helodermatid and varanid venoms (Table 1), and since studies have demonstrated the complexity and medicinal potential of Heloderma venom [27,[58][59][60][61][62][63][64] it is of interest to study the venom system of varanoid lizards in detail.
Toxins 2017, 9,242 6 of 38 morphological as well as molecular evolutionary studies have indicated that these glands are homologous with the venom glands of snakes [23,25,27,29,32,33,53]. The fact that varanid lizards possess highly developed dental glands suggests that those glands in one way or another play an important role in their life. It has been shown with venomous snakes which switch to eggs or constriction as an alternate form of prey capture rapidly lose the functionality of their venom glands, with atrophying occurring in short periods of (evolutionary) time [25,29,54,55]. The evolution of a complex venom system is likely only possible under certain contingent circumstances-i.e., when both environmental conditions select for it and a species' overall evolutionary trajectory facilitates it. For example, in Iguania the incipient venom glands never developed any significant complexity, probably due to the mainly insectivorous/herbivorous nature of these lizards. In addition, when a species evolves an alternative method of subduing prey that renders the venom system redundant, or switches to a diet with no need for subjugation (e.g., plants), the venom system often degrades [25,29,54,55]. The cost of venom production is presumably high enough to justify the presence of active secretory and delivery apparatuses only when it confers an evolutionary advantage [56]. It is notable, however, that in iguanian species which include a large quantity of vertebrates in their diet, the glands are larger and the protein-secreting region more developed (though we stress that this is not equivalent to them being 'venomous' per se but it is strongly suggestive of a functional role in predation) [29,30].
Previous work on helodermatid lizards has demonstrated a striking level of proteomic conservation within the venoms of this clade, with the same basic toxin groups present in similar quantities [27,57] despite the most recent common ancestor (MRCA) of the five Heloderma species existing ~15-20 million years ago [11]. While the venom glands of Varanus species have been compared transcriptomically, the only proteomic comparisons to-date were limited to SELDI mass spectrometry [23,32,33]. However a diversity of components have been functionally characterised from helodermatid and varanid venoms (Table 1), and since studies have demonstrated the complexity and medicinal potential of Heloderma venom [27,[58][59][60][61][62][63][64] it is of interest to study the venom system of varanoid lizards in detail. Figure 1. Organismal relationships, sizes and ecological niches occupied (blue = aquatic, brown = terrestrial, green = arboreal). Phylogeny based on [65][66][67]. Note that as it only includes species from this study, it excludes the other anguimorph lizards such as anguids that intervene between Heloderma and Lanthantus/Varanus. Red lineages are terrestrial, green are arboreal, and blue are aquatic. Therefore, the aim of this study was to undertake a proteomic and functional comparison of varanoid lizard venoms in order to uncover patterns of venom evolution across the clade and to add to the body of knowledge regarding their controversial evolution. In view of treating the varanid lizard venom apparatus as an integrated system, the teeth of each species were also examined. Varanid lizards belong to the lizard clade Anguimorpha that also contains anguid and helodermatid Figure 1. Organismal relationships, sizes and ecological niches occupied (blue = aquatic, brown = terrestrial, green = arboreal). Phylogeny based on [65][66][67]. Note that as it only includes species from this study, it excludes the other anguimorph lizards such as anguids that intervene between Heloderma and Lanthantus/Varanus. Red lineages are terrestrial, green are arboreal, and blue are aquatic. Therefore, the aim of this study was to undertake a proteomic and functional comparison of varanoid lizard venoms in order to uncover patterns of venom evolution across the clade and to add to the body of knowledge regarding their controversial evolution. In view of treating the varanid lizard venom apparatus as an integrated system, the teeth of each species were also examined. Varanid lizards belong to the lizard clade Anguimorpha that also contains anguid and helodermatid lizards [13,14,67], of which the gila monster (Heloderma suspectum) and Komodo dragon (V. komodoensis) are the most well-known species. Varanus is the sole extant genus in the family Varanidae and its closest extant relative is Lanthanotus borneensis (Borneo earless monitor-Lanthanotidae), also included in this study. Together Varanidae and Lanthanotidae (along with numerous extinct taxa) comprise the superfamily Varanoidea, of which Heloderma was previously also considered a member. Lizards in the genus Varanus are squamate reptiles with total body length ranging from 23 cm for adult V. brevicauda to over 3 m for V. komodoensis ( Figure 1).

Teeth Scanning Electron Microscopy
In the 20 species of varanoid lizard examined (Figure 1

Proteomics Studies
Proteomic analyses revealed the venoms to contain a large diversity of components, with kallikrein predominating in all species except V. griseus and V. varius, suggesting it plays a lesser role in these species (Figures 4 and 5, Table 2). Full MS/MS data are available in Supplementary Tables S1. 2D gels show that the kallikrein toxins possess significant diversity in isoelectric points and molecular weight diversity ( Figure 5), ranging from the narrow but dense spots of V. giganteus ( Figure 5B) and V. mertensi ( Figure 5E) through to the evolution of differential forms in V. prasinus ( Figure 5G) and extending to poorly staining low levels in V. griseus ( Figure 5D) and V. varius ( Figure 5L). There were significant kallikrein expression variations in closely related taxa, ranging from broad isoelectric point variation in V. scalaris ( Figure 5J) relative to V. tristis ( Figure 5K), through to downregulation in V. varius ( Figure 5L) relative to V. salvadorii ( Figure 5I).
L. borneensis venom displayed a 1D pattern more similar to the putatively ancestral type shared with Heloderma species (Figure 4) [57]. In contrast, the varanid 1D gels (Figure 4) revealed extensive variation is detectable outside regions of the gel corresponding to kallikrein (which remains relatively invariant across species). Chitinase/chitotriosidase were present in all the venoms of the dwarf species, but virtually absent from those of larger species (Figures 4 and 5, Table 2). point variation in V. scalaris ( Figure 5J) relative to V. tristis ( Figure 5K), through to downregulation in V. varius ( Figure 5L) relative to V. salvadorii ( Figure 5I).
L. borneensis venom displayed a 1D pattern more similar to the putatively ancestral type shared with Heloderma species (Figure 4) [57]. In contrast, the varanid 1D gels (Figure 4) revealed extensive variation is detectable outside regions of the gel corresponding to kallikrein (which remains relatively invariant across species). Chitinase/chitotriosidase were present in all the venoms of the dwarf species, but virtually absent from those of larger species (Figures 4 and 5, Table 2).

Kallikrein Molecular Evolution
The molecular phylogeny of the available anguimorph kallikrein sequences showed evidence of gene duplication and diversification ( Figure 6). In contrast to a previous study challenging kallikrein as a shared toxin type between snakes and lizard venoms based on the non-monophyly of toxicoferan oral gland sequences relative to non-toxicoferan oral gland and body tissue forms [40], we recovered all toxicoferan gland sequences as a monophyletic group relative to sequences from body tissues or from non-toxicoferan species. As the authors of that study did not provide their alignments or run files despite repeated requests by us, we cannot ascertain the cause of the discrepancy as due to either alignment issues with that previous study or differential phylogenetic methods. Regardless, our tree  [40].  Molecular modelling demonstrated the overall dN/dS value for anguimorph kallikrein toxins for which nucleotide data are publicly available was 0.80, indicating that the protein as a whole has been subject to neutral or weak purifying selection. Under negative or purifying selection, less "fit" nonsynonymous substitutions accumulate more slowly than synonymous substitutions, and under diversifying or positive selection, the converse is true. However, the FUBAR (Fast Unconstrained Bayesian AppRoximation) and MEME (Mixed Effects Model of Evolution) methods detected a number of individual sites on the toxin surface ( Figure 7) that are likely to have been subject to diversifying selection. This suggests that these sites may be important in the coevolutionary arms race between anguimorphs and their prey and may well play a role in the function of the toxins [87]. It should be noted that the analyses were conducted without the inclusion of a H. horridum sequence originally stated containing an insert characteristic only of snake venom forms [84] but which was subsequently shown to be erroneous [32]. However in the intervening time, this sequencing error led a study to conclude that convergent evolution had occurred between lizard and shrew toxin forms [88] and has been repeatedly cited as evidence for such convergence .
Outside of the core kallikrein enzyme, other widely expressed protein types included CRiSP, lysosomal acid lipase and PLA 2 . However, in contrast to the in general strong presence of kallikrein, these were more variable. Shotgun MS/MS revealed the presence of additional toxin types not discernible by gel-based methods such as the small natriuretic peptides ( Table 2).

Lysosomal Acid Lipase Molecular Phylogeny
Molecular phylogeny revealed the anguimorph lizard lysosomal acid lipase sequences to form a clade distinct from those sequenced from snake venom glands (Figure 8).

Bioactivity Testing
Bioactivity testing revealed dynamic variation in activities between and within each assay type with statistical significance between species for each assay ( Figure 9, Table 3

Kallikrein Enzymatic Activity upon Fluorescent Substrates
Consistent with the evidence for kallikrein molecular diversification, the ability of the venoms to cleave serine protease substrates varied significantly among the species studied but with phylogenetic patterns evident ( Figure 10). In addition, some species' venoms were potent on one substrate but not the other, with V. mitchelli the most potent upon both substrates.

Fibrinogen Cleavage Gels
As kallikrein toxins isolated from Heloderma venoms have been shown to exert non-clotting, destructive cleavage of fibrinogen [79,85], we investigated the relative time-dependent actions of venoms in this study for such actions (Figures 11-13). Intriguingly, differences in cleavage products indicated differential cleavage sites. A consistent pattern emerged in which all venoms displayed some ability to cleave the Aα chain, while the Bβ chain was cleaved more slowly and only fully destroyed by the most potent Aα chain acting venoms. High relative Aα chain activity predicted correspondingly high relative Bβ chain activity (PGLS: t = 7.3532, p = 7.969 × 10 −7 ). Aα chain cleavage was well predicted by activity upon the substrate S-2302 (PGLS: t = 2.8968, p = 0.009611, Figure 14) as was Bβ chain activity (PGLS: t = 2.8459, p = 0.01073, Figure 15). However, the γ chain was untouched even in species which fully consumed the Bβ chain.

Phospholipase A 2 Enzymatic Activity upon Fluorescent Substrate
V. varius venom had a dramatically higher level of PLA 2 enzymatic activity than other species tested (Figure 16), and thus obscured that other species displayed significant levels of PLA 2 activity consistent with the previous transcriptomic, proteomic or functional evidence for the presence of this enzyme type [23,32,86].

Rat Ileum Contraction Organ Bath Assay
V. varius venom demonstrated strong smooth muscle contracting activity ( Figure 17) consistent with the previously documented presence of AVIT and cholecystoxin peptides and the demonstrated action for snake venom homologues of this toxin type [32,33,113]. Due to lack of sufficient venom, other species were not tested.       Ancestral state reconstructions over branches comparing substrate consumption relative to alpha chain destruction. Bars indicate 95% confidence intervals for the estimate at each node. Note that due to the high dynamicity of venom evolution these quickly become broad as one moves down the tree. Phylogeny follows [65][66][67].

Figure 14.
Ancestral state reconstructions over branches comparing substrate consumption relative to alpha chain destruction. Bars indicate 95% confidence intervals for the estimate at each node. Note that due to the high dynamicity of venom evolution these quickly become broad as one moves down the tree. Phylogeny follows [65][66][67]. Toxins 2017, 9, 242 22 of 38 Figure 15. Ancestral state reconstructions over branches comparing substrate consumption relative to beta chain destruction. Bars indicate 95% confidence intervals for the estimate at each node. Note that due to the high dynamicity of venom evolution these quickly become broad as one moves down the tree. Phylogeny follows [65][66][67]. Figure 15. Ancestral state reconstructions over branches comparing substrate consumption relative to beta chain destruction. Bars indicate 95% confidence intervals for the estimate at each node. Note that due to the high dynamicity of venom evolution these quickly become broad as one moves down the tree. Phylogeny follows [65][66][67].

Figure 16.
Ancestral state reconstructions over branches for PLA2 enzyme substrate consumption where warmer colours represent more fluorescence production. Bars indicate 95% confidence intervals for the estimate at each node. Note that due to the high dynamicity of venom evolution these quickly become broad as one moves down the tree. Phylogeny follows [65][66][67].  where warmer colours represent more fluorescence production. Bars indicate 95% confidence intervals for the estimate at each node. Note that due to the high dynamicity of venom evolution these quickly become broad as one moves down the tree. Phylogeny follows [65][66][67].
Toxins 2017, 9, 242 23 of 38 Figure 16. Ancestral state reconstructions over branches for PLA2 enzyme substrate consumption where warmer colours represent more fluorescence production. Bars indicate 95% confidence intervals for the estimate at each node. Note that due to the high dynamicity of venom evolution these quickly become broad as one moves down the tree. Phylogeny follows [65][66][67].

Discussion
The venoms of varanoid lizards remain understudied in evolutionary toxinology; however, multiple sources of evidence point to the adaptive evolution of venom in varanoid lizards. The present study revealed the differential complexity of anguimorph lizard venoms across a wide taxonomical range ( Figure 1) and considered its evolutionary context as part of a combined predatory arsenal.
Prior to this study, tooth form and function have been reported for very few species of varanid lizards [1,114,115] and no study has compared tooth morphology and function across many species. The presence of serrations can be important indicators of tooth function in reptiles. Serrated teeth cut compliant material by presenting less contact area with the material, and thus the applied force at each point of contact is relatively greater. Serrations also act to trap and cut sections of material with a 'bite and slice' mechanism allowing a sliding force to be transferred into a cutting force [116]. Serrations have been reported for the teeth of V. komodoensis, V. salvator, and V. varius [1,114,115]. Juvenile V. niloticus also appear to bear serrations on the anterior and posterior edges, but these are lost in the adult form, to become smooth and blunt, reflecting a change in diet from insectivorous to molluscivorous in adults [117].
Slightly serrated teeth appear to be the plesiomorphic condition of varanid lizards (Figures 2  and 3). Significant diversification from the inferred plesiotypic tooth state was evident in two groups: the V. varius clade and members of the Odatria subgenus of dwarf monitors (V. acanthurus, V. baritji, V. gilleni, V. mitchelli, V. scalaris, and V. tristis) (Figure 3). Each of these groups also possesses morphological features and predatory strategies characteristic of the particular clade. V. komodoensis, V. salvadorii and V. varius had the most pronounced serrations, consistent with these species typically feeding on thick-skinned large mammalian prey that they dismember. Thus while there were not any evident direct links to venom chemical composition and serrations, there was for prey dismemberment strategies. The aquatic lineage V. mertensi thus represents a secondary loss of serrations. The dwarf species (Odatria) have also secondarily lost the serrations, with the sole exception within the clade being the arboreal lineage V. scalaris which has subsequently re-evolved serrations. The reasons for the unique (amongst Odatria) evolution of serrated teeth in V. scalaris may be linked to the inclusion of large, chitinous insects in its diet, which may require dismemberment before ingestion, or may be linked to active defence of territory from conspecifics [118]. In addition, this species predates upon the large treefrog Litoria splendida in its range, which it dismembers, eating only the legs and avoiding the large toxin-secreting parotid glands (Fry personal observations).
The extreme proteomic variability of varanid venoms, in contrast to the highly conserved nature of helodermatid venoms (Figures 4 and 5), taken in concert with the evidence of duplication and diversification of toxin genes such as kallikrein within Varanus (Figures 6 and 7) and other toxin types analysed previously [119]), is suggestive of adaptive evolution as a generator of toxin diversity in Varanus venoms.
The venoms studied contained myriad toxins, although with the exception of kallikrein none of those components were consistently highly expressed in the venom. Thus, a kallikrein-dominated venom is inferred to be reflective of the condition of the venom system of the anguimorph lizard MRCA. These functional and proteomic results are congruent with previous transcriptome studies which recovered kallikrein as the dominant toxin type [23,27,32,33]. Kallikrein enzymes have clearly undergone significant structural and functional diversification within the clade. There was evidence of gene duplication ( Figure 6) and diversification of the molecular surface biochemistry (Figure 7). The cleavage of the two fluorescent substrates were incongruent with each other, thus indicating variations in enzyme specificity ( Figure 9A,B and Figure 10). Similarly the actions upon fibrinogen produced variations in degradation products even when the same chains were targeted, also indicative of diversification of enzyme specificity (Figures 11-13). Thus, there is strong evidence for adaptive evolution of kallikrein toxins in the venoms of anguimorph lizards.
The venoms of larger species of monitor were generally seen to be more complex than that of their smaller congeners (Figures 4 and 5), which may be explained by the broader dietary range in larger species, both of the adult animals and across the life history. Thus, while multifunctional kallikrein enzyme toxins form the venom core, there has been additional diversification in the chemical composition of the venoms. As dietary studies of V. varius show [120], it has the broadest possible dietary range, feeding on everything from small invertebrates and eggs to medium-sized mammals, most likely experiencing ontogenetic niche shifts, which further necessitate adaptation to different prey items. This species also has the most distinct venom of all those in this study, with a secondary downregulation of kallikrein paralleled by an upregulation of other toxins which produce effects ranging from muscle contraction in this study ( Figure 15) through to platelet-inhibition mediated anticoagulation shown previously [33].
Smaller species (members of the Odatria subgenus), on the other hand, predominantly feed on lizards and insects throughout their life [121]. These species retain potent fibrinogen destruction activities in their kallikrein toxins with the exception of V. gilleni (Figures 10-12), the smallest species studied. It is notable that V. gilleni had the least amount of fibrinogen-targeting activity despite being rich in kallikrein toxins and having a stronger action upon substrate S-2302 than some other species which displayed stronger fibrinogen destruction activity (Figure 8). Other species were also rich in kallikrein on the gels but without corresponding activity upon either of the substrates tested. The kallikrein action in these species may instead be linked to high kininogen cleavage activity, with liberated kinins causing a rapid hypotensive effect in envenomed prey/predators. As this was beyond the scope of this study, future work on anguimorph lizard venoms should include assaying the relative release of kinins from kininogen by these venoms and the role this plays in the induction of hypotension.
Previous transcriptome studies revealed the presence of kallikrein transcripts as the dominant types in venom gland transcriptome of varanid lizards and snakes, and these toxins have been inferred as present in the oral glands of the Toxicofera MRCA [23,24,[27][28][29][30][32][33][34]. Characterised toxicoferan venom kallikreins induce fibrinogen depletion (Figures 11-13) in the prey organism and thus promote prolonged bleeding [24] as well as induction of hypotension through the release of kinins from kininogen [45]. Predators likely benefit from inducing blood loss or altering the blood pressure of their prey, as this will increase the chance of successful subjugation by weakening the prey. Previous work has shown that the Type III PLA 2 in V. varius venom block platelet aggregation, thus interfering with blood clotting via the same pathway as do homologous toxins in the venom of Heloderma species [33]. This species also had extremely high levels of PLA 2 activity in comparison to all other species (Figure 16). In the case of some monitor lizards, in particular the larger species such as V. varius or V. komodoensis, this type of toxic action might be beneficial even if a prey manages to escape the initial attack but succumbs to blood loss or shock in the aftermath. Field studies on V. komodoensis have indeed shown such post-bite mortality in a significant percentage (~20%) of prey animals [33,122].
Venom may have functions other than aiding in subjugation of prey. For instance, venom may help in antipredator defence or in intraspecific competition, and may have additional roles in digestion, by providing additional enzymatic activity, or in maintaining oral health by including antimicrobial agents [29,30,41,123,124]. Indeed, venoms are likely to be multifunctional in many (perhaps most) species with either the same or different toxins facilitating multiple ecological functions [26,29,30,41,123,125,126], and the range of roles venom may play in varanids and their relative importance remain uncertain. In addition to possible roles in predation, defence or competitor deterrence, the role of varanid oral secretions in aiding digestion (a typical function of vertebrate oral secretions) has previously been put forward as a testable hypothesis [123] and our data appear to support this additional role, at least for certain species. While the venoms of larger species showed evidence of derived activities, so did those of the dwarf species. Dwarf monitors, unique to Australia, grow no larger than 600 mm SVL-chitinase enzymes are likely helpful in the digestion of the thick exoskeletons of the arthropods that are their primary food source.
One salient hypothesis concerning effects of varanid lizard bites, such as hyperalgesia, inflammation and prolonged bleeding, is that the toxins (particularly those of dwarf species that are frequently predated upon by pythons) responsible for these effects may also serve defensive roles [10,123]. While kallikrein enzymes produce painful swelling, hyperalgesic cramping is a feature of the AVIT toxins type, which is also consistent with the strong muscular contraction observed for V. varius venom (Figure 17). Painful cramping would also aid in prey subjugation through loss of mobility. Another possible function of the secretions is intraspecific competition, particularly male-male combat, similar to the function of venom in platypus ecology [124]. In both such cases, rapid induction of hypotension by kinin release from kininogen would also be beneficial, thus reinforcing the multifunctional use of this functionally flexible enzyme type.
The presence of kallikrein alongside chitinase in the venoms of Odatria suggests that these secretions may serve multiple functional roles for these lizards-multifunctionality of toxic secretions is common in venomous organisms, and where venom is an oral secretion it often preserves its more primitive role in digestion [41]. Some components recovered in the study are very likely to be an adaptation for arthropod-based diets, which strengthens the point that venom glands in varanoid lizards likely have more than a single function. That is, in addition to its potential role in prey subjugation, defence, or intraspecific competition (the classic roles of "venom" sensu stricto [41]), varanoid lizard oral secretions may potentially aid in digestion.
According to our results, varanoid lizard venom is largely based on kallikrein toxins that previous studies have shown to be homologous with those present in the venom of advanced snakes [23,24,[27][28][29][30][32][33][34]. Additional components are present in various species with profile complexity seemingly being a function of size and habitat where larger monitors possess the most complex venom and smaller or aquatic species the least. The high level of variability of varanid venoms relative to the high levels of conservation in helodermatid lizards, points to active evolution under selection pressure. Components such as lysosomal acid lipase (Figure 8) are functionally uncharacterised but underscore the rich biodiscovery potential of lizard venoms. As we were able to examine the venoms of only 16 of the more than 60 extant species of varanid lizard, this study is clearly not a comprehensive investigation of the evolution of venom in Varanus, and future studies are very likely to provide additional novel insights. For example, it will be fruitful to investigate the venom profile of V. salvator, a large but predominantly aquatic lizard, and frugivorous species such as V. olivaceus. Nevertheless, we note that turning away from this fruitful area of research by denying the biochemical reality of lizard venom will hinder progress in this fascinating area. For example, because of their chain-selective fibrinogen-degrading activity, kallikrein enzymes from snake venom have been used as anticoagulant treatment of stroke, heart attack, and deep-vein thrombosis [85,127]. The functional diversity of such enzymes in the lizard venoms in this study underscores the rich biodiscovery and therapeutic potential of these novel natural products.

Species Studied
Australian samples were previously collected at the same time as a transcriptome study [32] under University of Melbourne (2005) approval UM0706247 as part of the long-term cryogenic collection of the Venom Evolution Laboratory, while non-Australian samples were supplied by Alphabiotoxine. Species studied were H. exasperatum (captive bred specimens of unknown founder stock), H. horridum (captive bred specimens of unknown founder stock), H. suspectum (captive bred specimens of unknown founder stock), Lanthanothus borneensis (captive bred specimens of unknown founder stock), Varanus acanthurus (Newman, WA, Australia), V. baritji (Adelaide River, NT, Australia), V. giganteus (Sandstone, WA, Australia), V. gilleni (captive bred specimens of unknown founder stock), V. griseus (captive bred specimen of unknown founder stock), V. jobiensis (captive bred specimens of unknown founder stock), V. melinus (captive bred specimens of unknown founder stock), V. komodoensis (Singapore Zoo captive bred specimens of unknown founder stock), V. melinus (captive bred specimens of unknown founder stock), V. mertensi (Kununurra, WA, Australia), V. prasinus (captive bred specimen of unknown founder stock), V. panoptes rubidus (Sandstone, WA, Australia), V. mitchelli (Kununurra, WA, Australia), V. prasinus (captive bred specimens of unknown founder stock), V. melinus (captive bred specimens of unknown founder stock), V. scalaris (Kununurra, WA, Australia), V. salvadorii (captive bred specimens of unknown founder stock), V. tristis (captive bred specimens of unknown founder stock), and V. varius (Mallacoota, VIC, Australia). 3 adult specimens for each species were pooled to minimize the effects individual variation. In order to remove mucous, all samples were filtered through 0.2 micron syringe filters prior to lyophilisation.

Scanning Electron Microscopy
Teeth were obtained from museum or frozen specimens and dissected out of both the top and bottom jaw. The largest tooth available (usually the ninth) was selected. All teeth were coated with a 20 nm layer of gold and imaged on a Phillips XL 30 scanning electron microscope. All images were taken at 10 kV and a working distance of 88 mm. These images were used for morphological measurements. To provide a comparative estimate of tooth serrations, the total distance from the tip of the tooth to the base was measured along both the anterior and posterior edge. The distance along this edge showing serrations was recorded, and divided by the total distance to calculate the proportion of the tooth edge that was serrated (% serrations).

1D Gel Electrophoresis
In order to establish the proteomics variation, 1D gradient gels were run under both reducing and non-reducing conditions using the manufacturer (BioRad) protocol. Gels were prepared as follows: 0.05 mL deionised H 2 O, 2.5 mL 30% acrylamide mix, 1. . Spreading gel was cast first. After it was set the spacer gel was slowly layered atop of it, and after spacer gel was set the stacking gel was layered atop of it. Running buffers were: 0.2 M Tris-HCl, pH 8.9 (anode buffer); 0.1 M Tris-tricine-HCl pH 8.45. The gels were run at 100 V for three hours at room temperature. 30 µg of venom was reconstituted in Tricine loading buffer (Bio-Rad) with 10 mM DTT added to provide reduced conditions. Gels were stained overnight with colloidal Coomassie brilliant blue G250 (34% methanol, 3% phosphoric acid, 170 g/L ammonium sulphate, 1 g/L Coomassie blue G250). After the staining was complete, water was used to remove excess dye.

2D Gel Electrophoresis
In order to further investigate the proteomics variation, particularly that of isoelectric variation, 2D gels were run using protocols previously optimised in the Fry lab [128][129][130]. 0.3 mg of venom sample were solubilized in 125 µL of rehydration buffer (8 M urea, 100 mM DTT, 4% CHAPS, and 0.5% ampholytes (Biolytes pH 3-10, Bio-Rad Lab)) with 0.01% bromophenol blue. The sample was mixed with shaking and centrifuged for 5 min at 4 • C, 14,000 rpm. This was done to remove any insoluble material. The supernatant was loaded onto IEF strips (Bio-Rad ReadyStrip, non-linear pH 3-10, 7 cm and 17 cm IPG) and left overnight for passive rehydration. Protein focussing was achieved via PROTEAN i12 IEF CELL (Bio-Rad Lab). The IEF running conditions were as follows: 100 V for 1 h, 500 V for 1 h, 1000 V for 1 h and 8000 V until 98,400 V/h. Actual current in the final step of the run varied in accordance to resistance. To each strip, a constant current of 50 µA was applied. After the run, IPG strips were incubated for 10 min in a reducing equilibration buffer (50 mM Tris-HCl, pH 8.8, 6 M urea, 2% SDS, 30% glycerol, 2% DTT) to reduce cysteine bonds. To alkylate reduced bonds, IPG strips were further incubated for 20 min in an alkylating equilibration buffer (50 mM Tris-HCl, pH 8.8, 6 M urea, 2% SDS, 30% glycerol, 2.5% iodoacetamide). After rinsing with SDS-PAGE running buffer, IPG strips were positioned on top of 12% polyacrylamide gels (Protean-II Plus, 18 × 20 cm, Bio-Rad Lab) using 0.5% agarose. Gels were run at 4 • C with a current of 10 mA per gel for 20 min followed by 20 mA per gel for the rest of the run until the bromophenol dye front was within 0.5 cm of the base of the gel. After the run, gels were briefly washed with water and stained with 0.2% colloidal Coomassie brilliant blue G250 overnight. Water was used to remove the excess of the dye after staining was complete. Visible spots were subsequently picked from gels and digested overnight at 37 • C with the use of sequencing grade trypsin (Sigma-Aldrich). Afterwards gel spots were washed with deionised H 2 O, destained (40 mM NH 4 CO 3 /50% acetonitrile (ACN)) and dehydrated (100% ACN), rehydrated in 10 µL of 20 µg/mL TPCK trypsin, and incubated at 37 • C overnight. To elute peptides, the following solutions were used per each gel spot: 20 µL of 1% formic acid (FA), followed by 20 µL of 5% ACN/0.5% FA. Collected peptides were put into MS vials and subjected to LC-MS/MS analysis.

Shotgun Sequencing
In order to identify low molecular weight peptides that do not resolve well on 1D or 2D gels, shotgun sequencing was used. 3 µg of crude venom sample was dissolved in 50 µL of 100 mM ammonium carbonate to reduce and alkylate cysteine bonds with subsequent addition of 50 µL of 2% iodoethanol/0.5% triethylphosphine in acetonitrile. The sample was afterwards resuspended in 20 µL of 40 mM ammonium bicarbonate before overnight incubation (at 37 • C) with 750 ng of sequencing grade trypsin (Sigma-Aldrich). To stop digestion 1 µL of concentrated formic acid was added to each of the samples. Samples were lyophilised then resuspended in 20 µL of 5% ACN/0.5% FA, put into MS vials and subjected to LC-MS/MS analysis.

LC-MS/MS
In order to identify the toxin types present, digested gel spots and digested whole venom (shotgun) samples were processed using an Agilent Zorbax stable bond C18 column (Agilent, Santa Clara, CA, USA) (2.1 mm by 100 mm, 1.8 µm particle size, 300 Å pore size) at a flow rate of 400 µL per minute and a gradient of 1-40% solvent B (90% acetonitrile, 0.1% formic acid) in 0.1% formic acid over 15 min or 4 min for shotgun samples and 2D-gel spots respectively on a Shimadzu Nexera UHPLC (Kyoto, Japan) coupled with an SCIEX 5600 Triple TOF mass spectrometer (Framingham, MA, USA). MS2 spectra are acquired at a rate of 20 scans per second with a cycle time of 2.3 s and optimised for high resolution. Precursor ions were selected between 80 and 1800 m/z with a 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. The isotopes within 2 Da were excluded for MS2. MS2 spectra were searched against known translated transcriptome libraries or UniProt database with Proteinpilot v4.0 (Sciex Framingham, MA, USA) using a thorough identification search, specifying iodoacetamide as an alkylation method, trypsin digestion and allowing for biological and chemical modifications (ethanolyl C or deamidated N in particular) and amino acid substitutions, including artefacts induced by the preparation or analysis processes. This was done to maximize the identification of protein sequences. In order to remove false positives, only 95% confidence value hits were examined.

RDES0011 Substrate
A working stock solution of freeze dried venom was reconstituted in a buffer containing 50% deionised H 2 O/50% glycerol (>99.9%, Sigma-Alrich) at a 1:1 ratio to preserve enzymatic activity and reduce enzyme degradation with the final venom concentration of 0.1 mg/mL, and then stored at

Phsopholipase A 2 Activity
We assessed the continuous PLA 2 activity of the venoms using a fluorescence substrate assay (EnzChek ® Phospholipase A 2 Assay Kit) (ThermoFisher Scientific, Sydney, Australia). A working stock solution of freeze dried venom was reconstituted in a buffer containing 50% deionised H 2 O/50% glycerol (99.9%, Sigma) at a 1:1 ratio to preserve enzymatic activity and reduce enzyme degradation with the final venom concentration of 0.1 mg/mL, and then stored at −20 • C. Venom solution (0.1 µg in dry venom weight) was brought up to 12.5 µL in 1× PLA 2 reaction buffer (250 mM Tris-HCL, 500 mM NaCl, 5 mM CaCl 2 , pH 8.9) and plated out in triplicates on a 384 well plate. Triplicates were measured by adding 12.5 µL quenched 1 mM EnzChek ® Phospholipase A 2 substrate per well (total volume 25 µL/well) over 100 cycles at an excitation of 485 nm and emission of 520 nm, using a Fluoroskan Ascent (ThermoFisher Scientific). The negative control consisted of PLA 2 reaction buffer and substrate only.

Rat ileum Organ Bath Testings
The rat ileum muscle preparations were isolated from adult male rats. The rats were euthanised by CO 2 asphyxiation. The isolated preparations were individually mounted in 15 mL parallel organ baths containing a Krebs solution with the following constituents (mM): NaCl, 118.4; KCl, 4.7; MgSO 4 , 1.2; KH 2 PO 4 , 1.2; CaCl 2 , 2.5; NaHCO 3 , 25 and glucose, 11.1. The Krebs solution was continuously bubbled with carbogen (95% O 2 and 5% CO 2 ) to maintain a pH between 7.2-7.4 at a temperature of 32-34 • C. A resting tension between 1 and 3 g was found to be the optimal starting baseline. Stimulation was performed with 50 µg/mL of crude venom; deionised H 2 O (170 µL) was used as a control. Human fibrinogen was reconstituted to a concentration of 2 mg/mL in isotonic saline solution, flash frozen in liquid nitrogen and stored at −80 • C until use. Freeze-dried venom was reconstituted in deionised H 2 O and concentrations were measured using a Thermo Fisher Scientific™ NanoDrop 2000 (Waltham, MA, USA). Assay concentrations were a 1:10 ratio of venom:fibrinogen, in comparison to 1:5 ratios used in snake venom testing [131]). The following was conducted in triplicate for each venom: Five "secondary" aliquots containing 10 µL buffer (5 µL of 4× Laemmli sample buffer (Bio-Rad, Hercules, CA, USA), 5 µL deionised H 2 O, 100 mM DTT (Sigma-Aldrich, St. Louis, MO, USA)) were prepared. A "primary" aliquot of fibrinogen (volume/concentration as per the above) was warmed to 37 • C in an incubator. 10 µL was removed from the primary aliquot ("0 min incubation" fibrinogen control) and added to a secondary aliquot, pipette mixed, and boiled at 100 • C for 4 min. Assay concentrations were a 1:10 ratio of venom:fibrinogen, in comparison to 1:5 ratios used in snake venom testing [131]). The following was conducted in triplicate for each venom: Five "secondary" aliquots containing 10 µL buffer (5 µL of 4× Laemmli sample buffer (Bio-Rad, Hercules, CA, USA), 5 µL deionised H 2 O, 100 mM DTT (Sigma-Aldrich, St. Louis, MO, USA)) were prepared. A "primary" aliquot of fibrinogen (volume/concentration as per the above) was warmed to 37 • C in an incubator. 10 µL was removed from the primary aliquot ("0 min incubation" fibrinogen control) and added to a secondary aliquot, pipette mixed, and boiled at 100 • C for 4 min. 4 µg (dry weight) of venom was then added to the primary aliquot of fibrinogen (amounting to 0.1 mg/mL of venom and 1 mg/mL of fibrinogen in 40 µL total volume), pipette mixed, and immediately returned to the incubator. At each incubation time period (1 min, 5 min, 20 min, and 60 min), 10 µL was taken from the primary aliquot, added to a secondary aliquot, pipette mixed, and boiled at 100 • C for 4 min. At each incubation time period (1 min, 5 min, 20 min, and 60 min), 10 µL was taken from the primary aliquot, added to a secondary aliquot, pipette mixed, and boiled at 100 • C for 4 min. The secondary aliquots were then loaded into the gels and were run in 1× gel running buffer at room temperature for 20 min at 90 V (Mini Protean3 power-pack from Bio-Rad, Hercules, CA, USA) and then 120 V until the dye front neared the bottom of the gel. Gels were stained with colloidal coomassie brilliant blue G250 (34% methanol (VWR Chemicals, Tingalpa, QLD, Australia), 3% orthophosphoric acid (Merck, Darmstadt, Germany), 170 g/L ammonium sulfate (Bio-Rad, Hercules, CA, USA), 1 g/L coomassie blue G250 (Bio-Rad, Hercules, CA, USA), and destained in deionised H 2 O.

Phylogenetic Comparative Analyses
A phylogeny was assembled using previous genetic studies [65][66][67] and was used for all further analyses conducted in R version 3.2.5 [132] using the ape package [133] for general handling of phylogenetic and trait data. Ancestral states were estimated and reconstructed over the tree in order to investigate the evolutionary history of the traits and consequently their relation to one another over time. The continuous functional traits were reconstructed by maximum likelihood in the contMap function in phytools [134]. We then fit PGLS models [135] in caper [136] to test for relationships.

Phylogenetic Reconstruction
Kallikrein, lysosomal acid lipase and chitinase datasets were analysed using Bayesian inference implemented on MrBayes, version 3.2.1 using lset rates = invgamma with prset aamodelpr = mixed, which enables the program to optimize between nine different amino acid substitution matrices. The analysis was performed by running a minimum of 10 million generations in four chains, and saving every 100th tree. The log-likelihood 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.

Molecular Modelling
Publicly available kallikrein sequences were retrieved from GenBank by using Varanus komodoensis kallikrein sequence as a query for a BLAST search within Anguimorpha [137][138][139]. Sequences with less than 70% coverage were discarded. Sequences were edited to only include the codons for the mature protein using AliView and aligned by using AliView's "Realign everything as translated amino acids" to translate the codons, align the resulting amino acids using MUSCLE, and reverse translating them [140,141]. MrBayes version 3.2 was used to create a phylogenetic tree of the sequences for performing the later selection analyses [142]. Using the closest sequence similarity reptile venom structure (GU441485) as input, the Phyre2 webserver generated a custom protein structure based on the published structure (PDB ID: 3S9C) [143]. This structure has a resolution of 1.8 Å and was 41% identical with the query sequence and with conservation of the cysteine residues. Protein models were rendered in UCSF Chimera version 1.10.2 [144]. Conservation scores were calculated using the UCSF Chimera implementation of AL2CO under the default settings [144,145]. Tests for selection were performed using HyPhy version 2.220150316: overall dN/dS value was calculated using the AnalyzeCodonData method, persistent site-by-site selection was analyzed with the FUBAR method, and episodic site-by-site selection was analyzed with the MEME method [146][147][148]. MEME is used for identifying sites that experience episodic selection pressures, where as FUBAR is improvement on site-wide selection analysis [148].