ABSTRACT

Venom spitting is a defence mechanism based on airborne venom delivery used by a number of different African and Asian elapid snake species (‘spitting cobras’; Naja spp. and Hemachatus spp.). Adaptations underpinning venom spitting have been studied extensively at both behavioural and morphological level in cobras, but the role of the physical properties of venom itself in its effective projection remains largely unstudied. We hereby provide the first comparative study of the physical properties of venom in spitting and non-spitting cobras. We measured the viscosity, protein concentration and pH of the venom of 13 cobra species of the genus Naja from Africa and Asia, alongside the spitting elapid Hemachatus haemachatus and the non-spitting viper Bitis arietans. By using published microCT scans, we calculated the pressure required to eject venom through the fangs of a spitting and a non-spitting cobra. Despite the differences in the modes of venom delivery, we found no significant differences between spitters and non-spitters in the rheological and physical properties of the studied venoms. Furthermore, all analysed venoms showed a Newtonian flow behaviour, in contrast to previous reports. Although our results imply that the evolution of venom spitting did not significantly affect venom viscosity, our models of fang pressure suggests that the pressure requirements to eject venom are lower in spitting cobras than in non-spitting cobras.

INTRODUCTION

A plethora of defensive behaviours can be found across the animal kingdom. Such variety can be explained by natural selection acting more strongly on defence mechanisms than on offence/predation mechanisms, as suggested by the ‘life–dinner principle’ (Dawkins and Krebs, 1979). According to this principle, evolutionary selective pressure on the prey is much stronger than on the predator, because in a predator–prey encounter, the prey may lose its life, while the predator may only lose a meal. Defensive behaviours can be summarised in three main categories: freezing, fleeing and active defence (Eilam, 2005). As part of the latter category, some organisms – for example, hymenoptera, arachnids and venomous snakes – employ venom, defined as an injectable harmful chemical secretion, to mount a more effective defensive attack. The noxious effects of venom increase the dissuading effect of the defence, enabling animals like bees, scorpions and snakes to ward off larger attackers (Schmidt, 2019). Although snake venoms are thought to have mainly evolved for their function in aiding predation (Arbuckle, 2017; Daltry et al., 1996), it is their use in defensive behaviour that makes them relevant to human health (Gutiérrez et al., 2017).

Snake venom consists of a complex mixture of peptides and proteins, small organic molecules and salts in an aqueous medium (Chan et al., 2016). The high peptide and protein content makes venom more viscous than water (Young et al., 2011), and it has been previously identified as a non-Newtonian shear-thinning fluid (Triep et al., 2013; Young et al., 2011). Venomous snakes (superfamily Colubroidea) inject venom into the body of their prey, or defensively into the body of their attackers, through specialised fangs or grooved teeth (Broeckhoven and du Plessis, 2017; Vonk et al., 2008). Members of the families Viperidae, Elapidae and Atractaspididae use an advanced front-fanged venom delivery system (Kerkkamp et al., 2017). In these snakes, the venom originates from the primary venom gland, and is expelled by the pressure of a skeletal muscle (referred to as m. compressor glandulae in viperids or m. adductor mandibulae externus superficialis in elapids; Haas, 1973) through the primary duct, the secondary (accessory) gland and into the fang, which acts like a hypodermic needle (Fransen et al., 1986; Jackson, 2003; Young and Kardong, 2007; Young et al., 2001). Once injected, venom toxins become systemic via dispersal by the bloodstream and lymphatic system, interacting with the prey/attacker's physiological proteins and receptors, ultimately disrupting the nervous system, the blood coagulation cascade, the cardiovascular and neuromuscular system, and/or homeostasis in general (Kerkkamp et al., 2017).

The Elapidae family of snakes includes taipans, mambas, coral snakes, kraits and cobras. Snakes of this family inject their venom through short, fixed fangs located in the frontal part of the upper jaw, as opposed to the movable front fangs of the Viperidae and Atractaspididae (Bogert, 1943; Vitt and Caldwell, 2013). Cobra species of the genus Naja Laurenti 1768 possess venoms with neurotoxic and/or cytotoxic properties, which they use to rapidly immobilise their prey for consumption, or to dissuade predators (Petras et al., 2011; Vitt and Caldwell, 2013). Members of this genus are present in both Africa and Asia (Vitt and Caldwell, 2013; Wüster, 1996; Wüster et al., 2007), and cobras from these two continents form phylogenetically distinct groups, which are thought to have separated about 16 Mya (Wüster et al., 2007).

Several Naja species are well known for their peculiar ability to spit venom as an exclusively defensive mechanism, expelling it as pressurised jets or sprays at their attackers (Berthé et al., 2009; Bogert, 1943; Panagides et al., 2017; Rasmussen et al., 1995; Westhoff et al., 2005; Wüster and Thorpe, 1992). These spits are generally aimed at the face and eyes of an aggressor (Westhoff et al., 2005), and once in contact with the eyes, can cause severe pain and inflammatory pathology (Chu et al., 2010; Westhoff et al., 2005). The ability to spit venom likely evolved from non-spitting ancestors on three independent occasions, once in African cobras and once in Asian cobras, and on a third occasion in the closely related rinkhals Hemachatus haemachatus (Kazandjian et al., 2021; Panagides et al., 2017; Slowinski et al., 1997; Wüster et al., 2007).

The venom delivery system of spitting cobras possesses several subtle morphological adaptations that enable them to eject their venom over long distances, and which distinguish them from non-spitting cobras. The discharge orifice, for example, has a more circular shape (Bogert, 1943; Wüster and Thorpe, 1992), and is directed more anteriorly, creating a 90 deg bend in the venom channel inside the fang (Balmert et al., 2011; Triep et al., 2013). This channel has internal ridges unique to spitting cobras (Berthé, 2011; Triep et al., 2013) that reduce the pressure loss by about 30% compared with an identical channel without ridges, thus helping to achieve a longer reach of the jet (Triep et al., 2013). Furthermore, spitting cobras actively displace the fang sheath (removing a physical barrier to venom expulsion), unlike other venomous snakes, where displacement of the fang sheath is passive (Young et al., 2004). Additional behavioural adaptations found in African spitting Naja species include adjusting head movements to distance from target to optimise the spread of venom (Berthé et al., 2009), and tracking and anticipating target movements to improve accuracy (Westhoff et al., 2010). Spitting cobras also show a certain degree of variation in their spitting modes: as demonstrated by previous studies (Rasmussen et al., 1995; Westhoff et al., 2005), some specialised spitters eject their venom in streams (e.g. Naja pallida) while others produce a fine mist (e.g. Naja nigricollis). The combination of morphological and behavioural adaptations allows most spitting cobras to eject venom up to at least 1 m, with some species (e.g. Naja mossambica) able to spit up to about 3 m (Rasmussen et al., 1995).

To date, considerable research effort has been focused on the anatomical features of the specialised venom delivery apparatus of spitting cobras (Bogert, 1943; Triep et al., 2013; Wüster and Thorpe, 1992; Young et al., 2004; 2009), and on their associated peculiar defensive behaviour (Berthé et al., 2009; Westhoff et al., 2005; 2010). In contrast, the possibility of changes in the composition of the venom itself, as an adaptation for its new role as a venom applied outside of the body, or toxungen (Nelsen et al., 2014), has remained largely neglected. Two recent studies have suggested that the venom of spitting species may have evolved for increased effectiveness when applied externally. Panagides et al. (2017) showed that African spitting cobras have venom that is more potently cytotoxic than that found in African non-spitters. Kazandjian et al. (2021) demonstrated that all three spitting lineages independently evolved venoms with more potent pain-inducing effects. These determine enhanced activation of sensory neurons through synergy between the ancestral cytotoxins widespread among cobras and phospholipases A2.

However, in addition to new selective pressures relating to its function as a toxungen, venom spitting may also have changed the mechanical demands of the venom, but so far this has not been studied. Since the venom has to pass through the narrow ducts of the venom apparatus, we expect that a lower venom viscosity (i.e. resistance to flow) would serve to reduce pressure loss during venom expulsion, thereby reducing the energetic requirements of ejection. Furthermore, for a given ejection force, venom projection distance would also be aided by more rapid expulsion, obtainable with a less viscous venom. On the other hand, in spitting cobras, a higher viscosity would aid jet cohesion after venom ejection, keeping the jet of venom from breaking up into droplets for longer, thus improving spitting distance and accuracy by reducing air drag. The reported strong shear-thinning, non-Newtonian behaviour of snake venom (Triep et al., 2013; Young et al., 2011) would result in a reduced viscosity in the high-shear environment of the venom channel, but a high viscosity in the low-shear environment of an airborne jet, and would thus likely aid in meeting these two seemingly conflicting demands.

Here, we measured and compared the rheological properties of the venoms of 12 spitting and non-spitting cobra species of the genus Naja from Africa and Asia, the only known ‘non-Naja’ species of spitting elapid, H. haemachatus, and the African non-spitting viperid Bitis arietans (used as an outgroup). We also compared the protein concentration and pH of the studied venoms, two properties known to play an important role in the stability of some snake venom components (Kurt and Aurich, 1976) and often directly correlated to the severity of the envenomation (Bon, 2003; Ribeiro et al., 2016; Sanhajariya et al., 2018).

Given the morphological differences between the fangs of spitting and non-spitting cobras (Bogert, 1943; Triep et al., 2013; Wüster and Thorpe, 1992; Young et al., 2004, 2009), we hypothesised that the two venom delivery mechanisms (i.e. spitting and biting) might be associated with different pressure requirements for venom ejection. Furthermore, we hypothesised that the venom of spitting cobras has a more pronounced shear-thinning behaviour than the venom of non-spitting cobras, in order to reduce pressure loss inside the venom duct and to increase jet cohesion in the airborne venom. To test this, we calculated and compared the pressure needed for venom to flow through the fang channel of one spitting and one non-spitting cobra species (Naja nigricollis and Naja nivea, respectively), using previously available microCT scanning data and our rheological data.

MATERIALS AND METHODS

Venom extraction

In total, venom samples of 30 snakes were used in this study. Venom was extracted from 28 cobras belonging to 13 different species of the genus Naja, namely: Naja annulifera Peters 1854, Naja atra Cantor 1842, Naja haje (Linnaeus 1758), Naja kaouthia Lesson 1831, Naja mossambica Peters 1854, Naja naja (Linnaeus 1758), Naja nigricollis Reinhardt 1843, Naja nivea (Linnaeus 1758), Naja nubiae Wüster & Broadley 2003, Naja pallida Boulenger 1896, Naja philippinensis Taylor 1922, Naja siamensis Laurenti1768 and Naja subfulva Laurent 1955. Venom was also extracted from one rinkhals, Hemachatus haemachatus Bonnaterre 1790, and one puff adder, Bitis arietans Merrem 1820, used for comparative analyses, as a ‘non-Naja’ venom spitter and non-spitter, respectively. Twelve of the specimens were captive bred (CB), while the remaining 18 were collected in the wild (see Table 1 for details). All snakes were maintained in individual cages within the temperature, humidity and light-controlled environment of the herpetarium at the Centre for Snakebite Research and Interventions, Liverpool School of Tropical Medicine. This facility and its protocols for the expert husbandry of the snakes are inspected and approved by the UK Home Office and the LSTM Animal Welfare and Ethical Review Board. Before the beginning of the experiments, none of the specimens considered for this study had been milked for at least 4 weeks. After milking, the snakes were immediately returned to their enclosures and the venom transferred into 2 ml low-protein binding cryotubes (Simport Scientific, Beloeil, Canada) using a pipette. Table 1 shows the average mass of fresh venom extracted from each specimen. The tubes were then transferred on ice to the laboratory of the Department of Materials Science and Engineering of the University of Sheffield for rheological, pH and concentration measurements on the same day.

Table 1.

Properties of the venom samples per specimen

Properties of the venom samples per specimen
Properties of the venom samples per specimen

Rheological tests

Shear viscosity measurements were performed in the Department of Materials Science and Engineering of the University of Sheffield, using a DHR-2 (TA Instruments, USA) rheometer, equipped with a cone-plate geometry (20 mm diameter, 1 deg angle cone, 27 μm truncation gap, 36 μl to fill), and subjecting samples to a shear rate ramp from 1.0 s−1 to 10, 000 s−1 (41 steps, 15 s per step), the maximum shear rate achievable by this instrument and geometry combination. Data below 100 s−1 were not included in later analysis as the apparent shear thinning observed is most likely attributed to surface tension effects and artefacts (see Fig. S1 and Ewoldt et al., 2015). Unless otherwise stated, all samples were tested at a room temperature of 25°C. This temperature was selected as it falls within the range of body temperatures of active snakes (El-Deib, 2005; Lillywhite, 2014) and approximates the temperature at which spitting was elicited from specimens of N. nigricollis, N. pallida, N. mossambica and H. haemachatus in previous studies (Westhoff et al., 2005; Young and O'Shea, 2005). Only species where sufficient venom was obtained to perform at least two replicates are shown. We were able to achieve up to three replicates for 19 of the 30 specimens included in this study. Venom samples that were not sufficient included H. haemachatus (African ‘non-Naja’ spitter), N. subfulva (African non-spitter) and N. naja (Asian non-spitter). In order to control for the potential presence of intraspecific variation in the considered rheological properties, all measurements were carried out on the venoms of single individuals, without pooling them.

Calculating fang venom shear rate

To support the range of shear rates tested and their biological relevance, it is necessary to calculate the natural range of shear rates encountered by venom. If venom is considered to be flowing down a channel, assuming all species spit in the same time and produce the same volume, the maximum shear strain rate at the fang wall is given by:
formula
(1)
where Q is the volumetric flow in m3 s−1, R is the radius of the venom channel in m and is the shear rate in s−1. According to data on N. pallida obtained by Triep et al. (2013) and du Plessis et al. (2018), the values considered during the venom spitting process are: volume of a single spitting event, Vsingle spit=1.0×10−8 m3; time for a single spitting event, tsingle spit=4×10−2 s; for B. arietans, R=3.8×10−4 m (du Plessis et al., 2018); for N. nigricollis, R=2.2×10−4 m (du Plessis et al., 2018); and for N. nivea, R=2.0×10−4 m (du Plessis et al., 2018). Therefore Q=2.5×10–7 m3 s−1 and using Eqn 1, for B. arietans, =5801 s−1; for N. nigricollis, =29,894 s−1 and for N. nivea, =38,051 s−1, which from a rheological perspective is in broad agreement with the 10,000 s−1 shear rate applied in this study.

Calculating the pressure needed to eject venom

If the venom is considered to be flowing down a venom channel of converging radius from R1 to R2, the pressure drop will be the result of the radius reduction from the fang base to the end of the fang where the exit orifice of the venom channel is located, plus the losses due to the viscous material (i.e. venom) flowing in the venom channel (Synolakis and Badeer, 1989). In order to corroborate if the flow is laminar or turbulent for the appropriate use of equations, the Reynolds number for the three species considered needs to be determined. The maximum Reynolds number defined for a Newtonian fluid can be calculated with the following equation:
formula
(2)
where Remax. is the maximum Reynolds number; ρ is the density of the venom=1084 kg m−3 (Triep et al., 2013); u1 is the venom velocity at the channel inlet=1.33 m s−1 (calculated with information from Triep et al., 2013); D1 is the diameter at the channel inlet=7.6×10−4 m for B. arietans (du Plessis et al., 2018), 4.4×10−4 m for N. nigricollis (du Plessis et al., 2018) and 4.0×10−4 m for N. nivea (du Plessis et al., 2018); μmin. is the dynamic viscosity of venom from our own data at 10,000 s−1=0.026±8.5×10−4 Pa s for B. arietans; 0.031±8.6×10−3 Pa s for N. nigricollis and 0.170±0.079 Pa s for N. nivea.

Assuming that all species have the same velocity at the channel inlet and density, Reynolds numbers are: B. arietans, Remax.=27.62; N. nigricollis, Remax.=23.25 and N. nivea, Remax.=4.24. All these Reynolds numbers are <100, corresponding to a laminar flow (in line with the predictions made by Triep et al., 2013), which is below the critical Reynolds number of 2300, above which turbulent flow is observed.

As the flow is in the laminar region, then the following equation, which corresponds to an Extended Generalised Bernoulli Equation, will be used to calculate the total pressure differential in the venom channel (see Appendix for detailed deduction of this equation):
formula
(3)
where ΔP is the pressure differential in the venom channel, in Pa; P1 and P2 are the pressures at the inlet and outlet points of the venom channel, respectively, in Pa; u1 and u2 are the velocities at the inlet and outlet points of the venom channel, respectively, in m s−1; ρ is the density of the venom, in kg m−3; A1 and A2 are the cross-section areas at the inlet and outlet points, in m2; Re is the Reynolds number; L is the length and D the average diameter of the venom channel, both in m; is the average velocity of the venom in the venom channel, in m s−1.

To directly relate these calculations to the natural system and the measured rheological data, microCT scans from du Plessis et al. (2018), and available at the GigaScience Database (http://dx.doi.org/10.5524/100389), were used to calculate venom channel length and radius. Fang morphology data was available for three species included in this study: B. arietans (viper), N. nigricollis (African spitting cobra) and N. nivea (African non-spitting cobra). MicroCT image stacks were imported into Amira (Thermo Fisher Scientific, v.2019.4) and 10 evenly spaced measurements were taken along the length of the venom channel (l), from the end of the entry orifice into the channel at the base of the fang to the opening point of the exit orifice at the tip of the fang. Of the 10 measurements per species, the average diameter was obtained (D) for input into Eqn 3. The values used for each variable for the three snake species are reported in Table 2.

Table 2.

Parameters used to calculate the pressure differential in the venom channel of the fang of Bitis arietans, Naja nigricollis and Naja nivea

Parameters used to calculate the pressure differential in the venom channel of the fang of Bitis arietans, Naja nigricollis and Naja nivea
Parameters used to calculate the pressure differential in the venom channel of the fang of Bitis arietans, Naja nigricollis and Naja nivea

Protein concentration

Protein concentration was measured for each venom sample using a UV300 Thermo Spectronic spectrometer (Unicam/Thermo, UK). All samples (dilutions consisting of 1.5 µl of fresh venom+1 ml of water) were analysed at room temperature in 1 cm path-length polystyrene cuvettes from 200 to 500 nm wavelength. Double distilled water was used as a blank and for all dilutions. Protein concentration was estimated as follows, using absorbance at 260 and 230 nm (Aitken and Learmonth, 2009):
formula
(4)
where A260 nm and A230 nm correspond to absorbance at 260 and 230 nm, respectively.

pH measurements

A Sentron pH meter (Netherlands) equipped with a cupFET pH probe was used to make pH measurements at room temperature. Two 3 μl droplets from each undiluted venom sample were measured individually and averaged to generate a pH measurement.

Phylogenetic comparative methods

The aim of the analyses reported here was to test for patterns in the measured parameters between spitting and non-spitting cobra venoms across the sampled species. All the analyses were performed using R 3.6.1 (https://www.r-project.org/) implemented using RStudio 1.2.1335, always taking the species phylogeny into account. We used the species tree reported in Kazandjian et al. (2021). This tree contained 46 elapid species belonging to 11 different genera and was generated using a multispecies coalescent model based on DNA sequence alignments of both mitochondrial (partial cytb and ND4 gene sequences) and nuclear genes (CMOS, NT3, PRLR, UBN1 and RAG1). For the analyses in the current study, we pruned the original tree and retained only the species used in the venom rheology tests (i.e. H. haemachatus and the various Naja species). The viper B. arietans was added manually to the tree as an outgroup, with branch lengths adjusted manually to reflect previous research suggesting that viperids separated from elapids about 61 Mya (Zheng and Wiens, 2016).

Within spitting cobras, a further division can be made in the different ways venom is ejected, which likely require different rheological properties of the venom. Following previous studies (Rasmussen et al., 1995; Westhoff et al., 2005), we divided the modes of venom ejection into three categories: (i) ‘streams’ where venom is ejected in the form of more or less continuous jets; (ii) ‘mist’ where venom is ejected in the form of a fine spray; (iii) ‘mixed’ where venom is ejected in a form in between the other two categories (see Table 1). Information about the venom spitting modes of seven species of spitting elapids considered in this study (N. atra, N. kaouthia, N. mossambica, N. nigricollis, N. pallida, N. siamensis and H. haemachatus) was gathered from the literature (Paterna, 2019; Rasmussen et al., 1995; Santra and Wüster, 2017; Westhoff et al., 2005). The spitting mode category for N. nubiae and N. philippinensis was assigned based on the authors' personal observations. The category ‘non-spitter’ was assigned to the non-spitting cobras N. annulifera, N. haje, N. naja, N. nivea and N. subfulva. The spitting mode category assigned to each studied species is reported in Table 1.

To first test if there was a difference between spitting and non-spitting cobras and/or between Asian and African cobras across all the measured physical properties, we performed a MANOVA using spitting behaviour (defined in the analysis as ‘spit’) as a binary factor (spitter or non-spitter), and the data about protein concentration and viscosity at 10,000 s−1 as multivariate dependent variables. We considered spitting behaviour as a binary trait only in this analysis. After this preliminary MANOVA, we performed the same test considering the three different spitting mode categories, in order to look for possible correlation between differences in spitting modes and the measured physical properties of the venoms.

To test if there was a difference in venom viscosity due to spitting behaviour, protein concentration or pH we performed an ANCOVA using viscosity at 10,000 s−1 (‘visc10000’) as dependent variable and ‘spit’, protein concentration (‘ProtConc’) and pH as independent variables.

To test if there was a difference in protein concentration due to spitting behaviour, we performed an ANCOVA using protein concentration as dependent variable and ‘spit’ as independent variable. We looked for possible presence of phylogenetic signal for pH, protein concentration and viscosity at 10,000 s−1, calculating both Blomberg's K (Blomberg et al., 2003) and Pagel's λ (Pagel, 1999), using the R packages caper (https://cran.r-project.org/web/packages/caper/), geomorph (https://cran.r-project.org/web/packages/geomorph/) and phytools (https://cran.r-project.org/web/packages/phytools/). Finally, we calculated Blomberg's K for protein concentration and viscosity at 10,000 s−1 at the same time.

RESULTS

Physical properties of the venom

For all Naja venoms tested, the protein concentrations had an average of 132.6 mg ml−1, ranging from 51.11 mg ml−1 (N. nivea) to 159.1 mg ml−1 (N. annulifera). The venoms of B. arietans and H. haemachatus had similar protein concentrations (132.4 and 132.5 mg ml−1, respectively). No significant differences were found between species or groups (Fig. 1B, see also Table 1 and more details below). The same was also true following quantification of venom pH, where the average pH of the Naja venoms was 5.77, ranging from 5.49 (N. kaouthia) to 6.02 (N. pallida). The pH of H. haemachatus venom was 5.76, and finally the pH of B. arietans venom was the lowest at 5.43 (Fig. 1C).

Fig. 1.

Physical properties ofsnakevenoms. (A) Cladogram of the elapid species analysed, extrapolated from the phylogenetic analyses performed (following Zheng and Wiens, 2016, viperids separated from elapids ∼61 Mya, therefore Bitis arietans has not been included in the cladogram). (B) Box plot of protein concentration for venoms extracted for each species. (C) Box plot of venom pH, where each datapoint represents the mean of two individual measurements. Triangles represent African Naja spp., diamonds represent Asian Naja spp. Venom-spitting species are in blue, non-spitting species in violet. The green circle and the green star represent, respectively, the spitting elapid Hemachatus haemachatus and the non-spitting viper outgroup B. arietans.

Fig. 1.

Physical properties ofsnakevenoms. (A) Cladogram of the elapid species analysed, extrapolated from the phylogenetic analyses performed (following Zheng and Wiens, 2016, viperids separated from elapids ∼61 Mya, therefore Bitis arietans has not been included in the cladogram). (B) Box plot of protein concentration for venoms extracted for each species. (C) Box plot of venom pH, where each datapoint represents the mean of two individual measurements. Triangles represent African Naja spp., diamonds represent Asian Naja spp. Venom-spitting species are in blue, non-spitting species in violet. The green circle and the green star represent, respectively, the spitting elapid Hemachatus haemachatus and the non-spitting viper outgroup B. arietans.

Rheological tests demonstrated that, contrary to our starting hypothesis, the venoms of both spitting and non-spitting cobras show a Newtonian behaviour, at least over the range reported here (i.e. 100 to 10,000 s−1) (Fig. 2). No significant differences between species or groups were evident (Table 1 and below).

Fig. 2.

Rheological properties of snake venoms. (A) Box plot of viscosity at 10,000 s−1 for venoms extracted from each species. (B) Viscosity as a function of shear rate for each species. Note the absence of data for H. haemachatus, N. subfulva and N. naja (venom volume insufficient to run the experiments). The same colour code used in Fig. 1 has been applied. Data are means±s.e.m. from at least two experiments per specimen.

Fig. 2.

Rheological properties of snake venoms. (A) Box plot of viscosity at 10,000 s−1 for venoms extracted from each species. (B) Viscosity as a function of shear rate for each species. Note the absence of data for H. haemachatus, N. subfulva and N. naja (venom volume insufficient to run the experiments). The same colour code used in Fig. 1 has been applied. Data are means±s.e.m. from at least two experiments per specimen.

Combining rheological and morphological data to determine the pressure required for venom to flow down the venom channel, Fig. 3 shows the results for the African non-spitting cobra N. nivea, the African spitting cobra N. nigricollis and the viper B. arietans. MicroCT scans obtained from du Plessis et al. (2018) indicate two different types of fangs, closed fused (B. arietans) and non-fused (N. nigricollis and N. nivea, Fig. 3A), and subsequent measurements provide information as to the fang length/diameter ratio (Fig. 3B). The results of fang pressure calculations shown in Fig. 3C report that the highest value corresponds to the non-spitter N. nivea (2.8×106 Pa), while the spitter N. nigricollis presents a lower value (0.17×106 Pa). The viper B. arietans shows the lowest pressure differential (0.10×106 Pa). The pressure differential results for the three snake species are reported in Table 2.

Fig. 3.

Fang pressure prediction for Naja nivea, Naja nigricollis and Bitis arietans. (A) MicroCT images showing fang types (data analysed from du Plessis et al., 2018; available at GigaScience Database, http://dx.doi.org/10.5524/100389) in the three species. (B) Fang length/diameter ratio. (C) ΔP in the fang venom channel, calculated using representative rheological data for each species.

Fig. 3.

Fang pressure prediction for Naja nivea, Naja nigricollis and Bitis arietans. (A) MicroCT images showing fang types (data analysed from du Plessis et al., 2018; available at GigaScience Database, http://dx.doi.org/10.5524/100389) in the three species. (B) Fang length/diameter ratio. (C) ΔP in the fang venom channel, calculated using representative rheological data for each species.

Phylogenetic comparative methods

The results of both MANOVAs showed no significant relationships between spitting behaviour and the multivariate combination of the measured physical properties of the venom (protein concentration, viscosity at 10,000 s−1). An additional MANOVA including pH among the variables was also performed, but then discarded because of the non-significance of the added variable and to simplify the model. The results of the ANCOVAs also showed no significant effect of spitting behaviour, protein concentration or pH on viscosity, or of spitting behaviour on protein concentration. Results of the statistical analyses performed considering the three spitting mode categories are reported in Table 3.

Table 3.

Results of statistical testing

Results of statistical testing
Results of statistical testing

Protein concentration, pH and viscosity at 10,000 s−1 show both Blomberg's K and, particularly, Pagel's λ close to 0 (Table 4), indicating phylogenetic independence (Karatzas and Shreve, 1998). The same can be said for the multivariate analysis, which takes into account both protein concentration and viscosity, and for which only Blomberg's K has been calculated. None of these results was significant, with P-values always higher than 0.05 (between 0.276 and 0.707 for K and equal to 1 for λ).

Table 4.

Results of phylogenetic signal testing

Results of phylogenetic signal testing
Results of phylogenetic signal testing

DISCUSSION

Young's study on venom gland pressure in spitting cobras suggested that the force required by the m. adductor mandibulae externus superficialis to expel venom would be reduced if a highly shear-thinning venom was present (Young et al., 2004). The sudden increase in shear rate upon entering the venom channel would cause a decrease in the viscosity of the venom, which could therefore be pushed through the fang more easily and thus at the higher velocities that are required to increase the reach of the venom jet (Triep et al., 2013). However, upon exiting the fang, the effective shear rate in the airborne venom jet ejected by a spitting cobra would be dramatically reduced, and as such, a higher viscosity in the jet would reduce internal flow, thus slowing down the breaking up of the jet into separate droplets. This provides the advantage of a more coherent jet of venom, resulting in less drag and thus a longer reach. Given that non-spitting cobras do not eject their venom, they presumably have less need for a higher venom ejection speed, and hence less need for a highly shear-thinning venom. In light of these biomechanical considerations, we expected a more pronounced shear-thinning behaviour in spitting cobras than in non-spitting cobras, in order to reduce pressure loss inside the venom duct and to increase jet cohesion.

Thus, when considering the above and the specific morphological adaptations to spitting in spitting cobras, such as the ridges present along the channel inside their fangs (Berthé, 2011; Triep et al., 2013), the more circular and anteriorly oriented discharge orifice of their fangs (Bogert, 1943; Wüster and Thorpe, 1992; Young et al., 2004) and the apparently higher algesic activity of venoms of the three spitting lineages (Kazandjian et al., 2021), we expected the rheological properties of the venom between spitting and non-spitting cobras to also be different. Hence, in light of our findings, it is surprising to find no systematic differences in venom viscosity between spitting and non-spitting species. However, it is worth noting that this result might be influenced by the small number of rheological tests performed for most of the analysed snakes, owing to the relatively small amount of venom a single cobra specimen produces.

Nevertheless, we did find differences in viscosity between and within species, suggesting that there is enough variability for natural selection to potentially act on. Between species, we found that the average venom viscosities at 10,000 s−1 went from a minimum of 0.0103 Pa s (N. naja) to a maximum of 0.1709 Pa s (Naja nivea) (see Fig. 2A and Table 1). Similarly, we found that viscosity could vary greatly even among specimens of the same species. For instance, the average venom viscosities measured for the three N. nubiae specimens (NajNubCB001, NajNubCB003 and NajNubCB004) were, respectively, 0.0064, 0.0252 and 0.0790 Pa s (Table 1 and Table S1). These results suggest that the venom of all the elapid species we analysed may vary in its viscosity owing to functional or other non-flow related requirements. We speculate that, within the range of rheological variability we recovered here for spitting cobras, other selective pressures may dictate the observed rheological properties. Although protein concentration and pH have been previously shown to vary and be of influence in snake venoms (Takahashi and Ohsaka, 1970) and in other secreted protein systems (e.g. silk, Holland et al., 2007; Terry et al., 2004), these two parameters did not vary significantly in our study.

Snake venom is known to vary in composition depending on different factors, such as diet (Daltry et al., 1996; Gibbs et al., 2011), ontogeny (Alape-Girón et al., 2008; Cipriani et al., 2017; Mackessy et al., 2006) and, potentially, local adaptation driven by relatively small changes in the physical environment (Zancolli et al., 2019). Compositional alterations in snake venom likely influence its rheology. Environmental changes determined by captivity (e.g. food supply restricted to a single type of prey) can also result in modifications of venom composition. However, most of the evidence produced so far suggests that the effect of captivity on snake venom composition is minimal (Farias et al., 2018; Freitas-de-Sousa et al., 2015; McCleary et al., 2016). In light of this, and considering that all venom samples analysed here were sourced from adult snakes fed on the same diet and kept under the same enclosure conditions, age, diet and ecology-related sources of variability have been minimised as much as possible, and thus seem unlikely to play a major role in the findings of this study. Thus, we suspect inherited differences in molecular venom composition (Mukherjee and Maity, 2002; Silva-de-França et al., 2019; Tan and Tan, 1988) to be the primary influence for any rheological differences. However, considering that both Petras et al. (2011) and Kazandjian et al. (2021) found the venoms of African spitting cobras (N. katiensis, N. mossambica, N. nigricollis, N. nubiae, N. pallida) to show similar compositional patterns in terms of proteins, we speculate that long chain (high molecular weight) non-protein molecules present in snake venom, such as carbohydrates (Bieber, 1979; Gowda and Davidson, 1992; Nawarak et al., 2004; Soares and Oliveira, 2009), could be responsible for the detected variation in rheological properties.

Surprisingly, our rheological testing showed Newtonian behaviour for all analysed snake venoms across the shear rates presented. This appears to be in direct contrast to previous studies where snake venoms have been classified as non-Newtonian (Balmert et al., 2011; Triep et al., 2013; Young et al., 2011). For example, Triep et al. (2013) suggested that N. pallida venom had non-Newtonian behaviour in the range of 1 to 37 s−1. However, upon closer inspection of the data within this range, we conclude that the apparent shear-thinning behaviour of N. pallida venom could be attributed to surface tension effects (Ewoldt et al., 2015). As a result, through comparison of our findings to previous studies, and accounting for the potential confounding influence of surface tension artefacts, we propose that any venom rheological data obtained below 100 s−1 presented to date should not be considered when determining if a venom is Newtonian or non-Newtonian (see Fig. S1). Previous studies have interpreted the rheological behaviour of snake venom based on experimental shear rate values ranging from 1 to 100 s−1 (Triep et al., 2013), and from 0.01 to 200 s−1 (Young et al., 2011). In these cases, we suggest that, because of the surface tension artefacts, only data from 100 to 200 s−1 (indicating a Newtonian flow behaviour) should be considered.

To explore the delivery mechanism and pressure requirements of venom ejection, we combined our rheology data with microCT scans of snake fangs reported by du Plessis et al. (2018). For the corresponding calculations, since fang venom channels are typically slightly curved and may have additional pressure-increasing features such as internal ridges (Berthé, 2011; Triep et al., 2013) and pressure losses due to viscosity, an extended generalised Bernoulli equation (Eqn 3) was used. We were able to model the pressure required for venom to flow through the fang for three of the species we studied: N. nivea (African non-spitting cobra), N. nigricollis (African spitting cobra) and B. arietans (viper). Despite the limited number of species investigated, there are clear differences in the pressure required to move venom down the fang. The spitter N. nigricollis has a smaller fang length/diameter ratio and a lower pressure requirement, whereas the non-spitter N. nivea has a larger fang length/diameter ratio and a higher pressure requirement. Interestingly, the viper B. arietans displayed both the largest fang length/diameter ratio and the lowest pressure requirement overall (Fig. 3), likely related to the relatively larger absolute diameter and/or curvature of the fang channel in this species. We found that the effect of viscosity and friction of the fluid in the venom channel (which is included in the Reynolds number; see Appendix for details) represents 5% of the pressure loss in B. arietans; 17% in N. nigricollis (spitter); and 9% in N. nivea (non-spitter). It appears that with this approach neither density nor viscosity contributes significantly to pressure losses, and that the major influence is the cross-section area variations along the venom channel (A1>A2), which represent between 83 and 95% of the total pressure loss. In light of this, we conclude that for all the viscosities observed, and for all the snake species analysed in this study, venom viscosity does not strongly influence the pressure requirements of venom ejection, and that what most defines such requirements are the morphological adaptations of the venom delivery systems (i.e. tapering of the fang venom channel).

Considering the ‘life–dinner principle’ (Dawkins and Krebs, 1979), which suggests that selection for defensive strategies should take precedence over selection for predatory efficiency, the lack of significant signs of adaptation of venom rheological properties to spitting behaviour is unexpected. In fact, if the principle is true, considering the lack of consistent differences in venom rheology between spitting and non-spitting cobras, and that venom spitting is an unambiguously defensive behaviour, it is interesting to question why selective pressures have not favoured the emergence of venom spitting in all cobras.

A recent study investigating patterns of venom-induced pain across snake species and time has suggested that the common ancestor of all elapids might have possessed early-pain-inducing venom (Ward-Smith et al., 2020). With the rapid infliction of pain being a requirement of defensive venoms (Eisner and Camazine, 1983; Ward-Smith et al., 2020), this could indicate that the use of venom for defensive purposes appeared early in elapid evolution, before the evolution of spitting behaviour. While a trend towards loss of rapidly painful venom is common in snakes (Ward-Smith et al., 2020), venom spitting, coupled with enhanced algesic activity (Kazandjian et al., 2021) could be an extension of this basic defensive strategy (i.e. injection of early-pain-inducing venom), which allows contactless defence at a distance, and of shorter duration and higher accuracy than striking/biting (Kardong and Bels, 1998; Westhoff et al., 2010; Young et al., 2001). In this scenario, spitting behaviour is probably the evolutionary response to specific selective pressures. Exposure to agile vertebrates (including visually acute primates, as suggested by Kazandjian et al., 2021), likely attacking from an elevated position, and for which a defensive strategy involving striking/biting could be hazardous and/or ineffective, could have been one of the drivers of spitting behaviour evolution. It is therefore possible that spitting behaviour would not emerge in the absence of this kind of selective pressures, thus offering a conjecture for why not all cobra species are able to spit venom. Alternatively, the existence of yet unidentified constraints preventing the evolution of spitting in non-spitting species is not to be excluded a priori.

Spitting behaviour has been recently documented for two species of Asian cobras that are generally considered non-spitters and that display very limited modification of their fangs, namely N. kaouthia and N. atra (Paterna, 2019; Santra and Wüster, 2017; Wüster and Thorpe, 1992). These reports suggest that venom spitting can evolve in the presence of very limited adaptation of the dentition, without the greater level of morphological adaptation and precision documented for specialised spitters (Triep et al., 2013; Young et al., 2004). The reason why these species have not evolved the more specialised venom spitting apparatus that other species possess (e.g. N. mossambica, N. nigricollis, N. pallida), may be due to differences in selective pressures, as outlined above, or perhaps the more recent origin of spitting in Asian cobras (Kazandjian et al., 2021). In light of these findings, spitting behaviour in cobras should probably not be seen as a binary trait, but may vary continuously in prevalence among the species of the genus Naja. Understanding the evolution, or lack of evolution, of specialised spitting behaviour and associated physical adaptations would likely require studying the efficacy and prevalence of spitting behaviour as a defence against natural predators, an under-documented aspect in the literature on this adaptation.

Although, perhaps surprisingly, our results did not show any clear adaptation of the rheological properties of venom to spitting behaviour, we demonstrated that both spitting and non-spitting cobra venoms are Newtonian fluids over a biologically relevant shear rate range, in contrast to previous literature reports. In order to gain a more comprehensive understanding of the mechanics behind venom spitting in cobras, we suggest considering the continuous nature of the prevalence of spitting behaviour and spitting modes, fang morphology and parts of the cobra venom delivery system at play in venom spitting but not included in this study (e.g. m. adductor mandibulae externus superficialis). Furthermore, future studies should increase the sample size in terms of both venom samples, specimens and species, in order to more comprehensively address the remarkably high variability in viscosity we detected in the present work. We hope our findings will stimulate further comparative study of the rheology of venom spitting across the genus Naja.

APPENDIX

Delta pressure equation

There is pressure loss in fangs associated to converging diameter, which means rbase of the fang > rend of the fang and close to the exit orifice, which is in line with our fang measurements using microCT data (data analysed from du Plessis et al., 2018; available at GigaScience Database, http://dx.doi.org/10.5524/100389). However, that is not the only effect in pressure loss, because there is the effect of venom flowing in the venom channel, i.e. viscous pressure loss. Therefore, Poiseuille's law is not correct in this case because the diameter of the venom channel is not constant, and Bernoulli's equation is only accepted if there is no viscous pressure loss. Therefore, an extended generalised Bernoulli equation must be used in order to have an approximation of the pressure loss in the venom channel considering radius variations and viscosity (Synolakis and Badeer, 1989).

If the venom channel is considered as a converging radius pipe (see Fig. S2), then the generalised Bernoulli's equation considered for the venom channel can be written as:
formula
(A1)
where P1 and P2 are the pressures at the inlet and outlet points, in Pa; u1 and u2 are the velocities at the inlet and outlet points, in m s−1; ρ is the density of venom, in kg m−3; hf corresponds to losses due to viscosity, in m; g is the acceleration of gravity, 9.81 m s−2.
hf can be expressed as defined by Soares and Santos (2013), as follows:
formula
(A2)
where f is the friction factor; l is the length of the venom channel, in m; D is the mean diameter of the venom channel, in m; is the mean velocity of the venom in the venom channel, in m s−1, and can be calculated with the following equation:
formula
(A3)
where Q is the volumetric flow in the venom channel, in m3 s−1; is the mean cross section area of the venom channel, in m2. The friction factor, for laminar flow, can be expressed as:
formula
(A4)
If we combine Eqns A2 and A4, we obtain:
formula
(A5)
And combining Eqns A1 and A5:
formula
(A6)
From the continuity equation (Munson et al., 2006):
formula
(A7)
where A1 and A2 are the cross-section areas at the inlet and outlet points, in m2.
Rearranging Eqn A7:
formula
(A8)
If we define P1P2P, rearrange Eqn A6, and combine with Eqn A8, we obtain Eqn 3:
formula
(A9)
which is the equation used to calculate the pressure loss in the venom channel.

Acknowledgements

The authors thank Andreas Koeppel (University of Sheffield, Department of Materials Science and Engineering, Sheffield, UK) for his help with some of the rheological measurements, and Professor Anton du Plessis (Research group 3D Innovation, Stellenbosch University, Stellenbosch, South Africa) for providing microCT details of his published paper, which were used to analyse microCT images included in this study. IA thanks Drs Bart Hallmark, Simon Butler and Ian Wilson (Department of Chemical Engineering and Biotechnology, University of Cambridge, UK) for their help during a pilot study on cobra venom rheology, and Pedro Coelho and Yuri Simone (CIBIO/InBIO, Centro de Investigação em Biodiversidade e Recursos Genéticos da Universidade do Porto, Vairão, Portugal) for providing useful comments. The authors thank two anonymous reviewers for their insightful comments which improved the manuscript.

Footnotes

Author contributions

Conceptualization: I.A., C.H., A.v.d.M.; Methodology: I.A., E.B.-L., W.W., C.H., A.v.d.M.; Validation: E.B.-L., C.H., A.v.d.M.; Formal analysis: I.A., E.B.-L., W.W., C.H., A.v.d.M.; Investigation: I.A., E.B.-L., C.H.; Resources: E.B.-L., N.R.C., R.A.H., P.D.R., E.C., W.W., C.H., A.v.d.M.; Data curation: I.A., E.B.-L., P.D.R., E.C.; Writing - original draft: I.A., E.B.-L., N.R.C., R.A.H., W.W., R.C., C.H., A.v.d.M.; Writing - review & editing: I.A., E.B.-L., N.R.C., R.A.H., P.D.R., E.C., W.W., R.C., C.H., A.v.d.M.; Visualization: I.A., E.B.-L., N.R.C., R.A.H., W.W., R.C., C.H., A.v.d.M.; Supervision: C.H., A.v.; Project administration: I.A., C.H., A.v.d.M.; Funding acquisition: I.A., A.v.d.M.,C.H.

Funding

I.A. gratefully acknowledges the financial support provided by The Company of Biologists in the form of a Travelling Fellowship (JEBTF 171127). I.A. is supported by a Fundação para a Ciência e a Tecnologia, Portugal (FCT) Doctoral Fellowship (SFRH/BD/137797/2018). A.v.d.M. is financed through FCT under contract number DL57/2016/CP1440/CT0009. E.B.L. acknowledges the CONACYT (Consejo Nacional de Ciencia y Tecnologia, Mexico) for its financial support (472130). N.R.C., R.A.H. and W.W were supported in part by Leverhulme Trust grant RPG-2012-627.

References

Aitken
,
A.
and
Learmonth
,
M. P.
(
2009
).
Protein determination by UV absorption
. In
The Protein Protocols Handbook
(ed.
J. M.
Walker
), pp.
3
-
6
.
Totowa
:
Humana Press
.
Alape-Girón
,
A.
,
Sanz
,
L.
,
Escolano
,
J.
,
Flores-Diaz
,
M.
,
Madrigal
,
M.
,
Sasa
,
M.
and
Calvete
,
J. J.
(
2008
).
Snake venomics of the lancehead pitviper Bothrops asper: geographic, individual, and ontogenetic variations
.
J. Proteome Res.
7
,
3556
-
3571
.
Arbuckle
,
K.
(
2017
).
Evolutionary context of venom in animals
. In
Evolution of Venomous Animals and Their Toxins
(ed.
P.
Gopalakrishnokone
and
A.
Malhotra
), pp.
3
-
31
.
Dodrecht
:
Springer
.
Balmert
,
A.
,
Hess
,
D.
,
Brücker
,
C.
,
Bleckmann
,
H.
and
Westhoff
,
G.
(
2011
).
Spitting cobras: fluid jets in nature as models for technical applications
.
Bioinspiration, Biomimetics, and Bioreplication. Proc. SPIE
7975
,
797514
.
Berthé
,
R. A.
(
2011
).
Spitting behaviour and fang morphology of spitting cobras
.
PhD thesis
,
Rheinische Friedrich-Wilhelms Universität Bonn
,
Bonn
,
Germany
.
Berthé
,
R. A.
,
de Pury
,
S.
,
Bleckmann
,
H.
and
Westhoff
,
G.
(
2009
).
Spitting cobras adjust their venom distribution to target distance
.
J. Comp. Phsiol. A
195
,
753
-
757
.
Bieber
,
A. L.
(
1979
).
Metal and nonprotein constituents in snake venoms
. In
Snake Venoms
(ed.
C. Y.
Lee
), pp.
295
-
306
.
Heidelberg
:
Springer-Verlag
.
Blomberg
,
S. P.
,
Garland
,
T.
and
Ives
,
A. R.
(
2003
).
Testing for phylogenetic signal in comparative data: behavioral traits are more labile
.
Evolution
57
,
717
-
745
.
Bogert
,
C. M.
(
1943
).
Dentitional phenomena in cobras and other elapids with notes on adaptive modifications of fangs
.
Bull. Am. Museum Nat. Hist.
81
,
260
-
285
.
Bon
,
C.
(
2003
).
Pharmacokinetics of venom toxins and their modification by antivenom therapy
.
J. Toxicol. Toxin Rev.
22
,
129
-
138
.
Broeckhoven
,
C.
and
du Plessis
,
A.
(
2017
).
Has snake fang evolution lost its bite? New insights from a structural mechanics viewpoint
.
Biol. Lett.
13
,
20170293
.
Chan
,
Y. S.
,
Cheung
,
R. C. F.
,
Xia
,
L. X.
,
Wong
,
J. H.
,
Ng
,
T. B.
and
Chan
,
W. Y.
(
2016
).
Snake venom toxins: toxicity and medicinal applications
.
Appl. Microbiol. Biotechnol.
100
,
6165
-
6181
.
Chu
,
E. R.
,
Weinstein
,
S. A.
,
White
,
J.
and
Warrell
,
D. A.
(
2010
).
Venom ophthalmia caused by venoms of spitting elapid and other snakes: report of ten cases with review of epidemiology, clinical features, pathophysiology and management
.
Toxicon
56
,
259
-
272
.
Cipriani
,
V.
,
Debono
,
J.
,
Goldenberg
,
J.
,
Jackson
,
T. N. W.
,
Arbuckle
,
K.
,
Dobson
,
J.
,
Koludarov
,
I.
,
Li
,
B.
,
Hay
,
C.
,
Dunstane
,
N.
et al. 
(
2017
).
Correlation between ontogenetic dietary shifts and venom variation in Australian brown snakes (Pseudonaja)
.
Comp. Biochem. Phys. C
197
,
53
-
60
.
Daltry
,
J. C.
,
Wüster
,
W.
,
Thorpe
,
R. S.
(
1996
).
Diet and snake venom evolution
.
Nature
379
,
537
-
540
.
Dawkins
,
R.
and
Krebs
,
J. R.
(
1979
).
Arms races between and within species
.
Proc. Royal Soc. Lond. B
205
,
489
-
511
.
du Plessis
,
A.
,
Broeckhoven
,
C.
and
le Roux
,
S. G.
(
2018
).
Snake fangs: 3D morphological and mechanical analysis by microCT, simulation, and physical compression testing
.
GigaScience
7
,
1
-
8
.
Eilam
,
D.
(
2005
).
Die hard: a blend of freezing and fleeing as a dynamic defense - implications for the control of defensive behavior
.
Neurosci. Biobehav. Rev.
29
,
1181
-
1191
.
Eisner
,
T.
and
Camazine
,
S.
(
1983
).
Spider leg autotomy induced by prey venom injection: an adaptive response to “pain”?
Proc. Natl. Acad. Sci. USA
80
,
3382
-
3385
.
El-Deib
,
S.
(
2005
).
Serum catecholamine and hormonal titers in the hibernating snake Naja haje haje, with reference to the annual climatic cycle
.
J. Therm. Biol.
30
,
580
-
587
.
Ewoldt
,
R. H.
,
Johnston
,
M. T.
and
Caretta
,
L. M.
(
2015
).
Experimental challenges of shear rheology: how to avoid bad data
. In
Complex Fluids in Biological Systems. Experiment, Theory, and Computation
(ed.
S.
Spagnolie
), pp.
207
-
241
.
New York
:
Springer
.
Farias
,
I. B.
,
Morais-Zani
,
K.
,
Serino-Silva
,
C.
,
Sant'Anna
,
S. S.
,
Rocha
,
M. M. T. D.
,
Grego
,
K. F.
,
Andrade-Silva
,
D.
,
Serrano
,
S. M. T.
and
Tanaka-Azevedo
,
A. M.
(
2018
).
Functional and proteomic comparison of Bothrops jararaca venom from captive specimens and the Brazilian bothropic reference venom
.
J. Proteom.
174
,
36
-
46
.
Fransen
,
J. A. M.
,
Kardong
,
K. V.
and
Dullemeijer
,
P.
(
1986
).
Feeding mechanism in the rattlesnake Crotalus durissus
.
Amphib-Reptilia
7
,
271
-
302
.
Freitas-de-Sousa
,
L. A.
,
Amazonas
,
D. R.
,
Sousa
,
L. F.
,
Sant'Anna
,
S. S.
,
Nishiyama
,
M. Y.
Jr.
,
Serrano
,
S. M. T.
,
Junqueira-de-Azevedo
,
I. L. M.
,
Chalkidis
,
H. M.
,
Moura-da-Silva
,
A. M.
and
Mourao
,
R. H. V.
(
2015
).
Comparison of venoms from wild and long-term captive Bothrops atrox snakes and characterization of batroxrhagin, the predominant class PIII metalloproteinase from the venom of this species
.
Biochimie
118
,
60
-
70
.
Gibbs
,
H. L.
,
Sanz
,
L.
,
Chiucchi
,
J. E.
,
Farrell
,
T. M.
and
Calvete
,
J. J.
(
2011
).
Proteomic analysis of ontogenetic and diet-related changes in venom composition of juvenile and adult Dusky Pigmy rattlesnakes (Sistrurus miliarius barbouri)
.
J. Proteomics
74
,
2169
-
2179
.
Gowda
,
D. C.
and
Davidson
,
E. A.
(
1992
).
Structural features of carbohydrate moieties in snake venom glycoproteins
.
Biochem. Bioph. Res. Co.
182
,
294
-
301
.
Gutiérrez
,
J. M.
,
Calvete
,
J. J.
,
Habib
,
A. G.
,
Harrison
,
R. A.
,
Williams
,
D. J.
and
Warrell
,
D. A.
(
2017
).
Snakebite envenoming
.
Nat. Rev. Dis. Primers
3
,
17063
.
Haas
,
G.
(
1973
).
Muscles of the jaw and associated structures in the Rhyncocephalia and Squamata
. In
Biology of the Reptilia
, Vol.
4
(ed.
C.
Gans
and
T.
Parsons
), pp.
285
-
490
.
New York
:
Academic Press
.
Holland
,
C.
,
Terry
,
A. E.
,
Porter
,
D.
and
Vollrath
,
F.
(
2007
).
Natural and unnatural silks
.
Polymer
48
,
3388
-
3392
.
Jackson
,
K.
(
2003
).
The evolution of venom-delivery system in snakes
.
Zool Jour. Linn. Soc.
137
,
337
-
354
.
Karatzas
,
I.
and
Shreve
,
S. E.
(
1998
).
Brownian motion
. In
Brownian Motion and Stochastic Calculus
, pp.
47
-
127
.
New York
:
Springer
.
Kardong
,
K. V.
and
Bels
,
V.
(
1998
).
Rattlesnake strike behavior: kinematics
.
J. Exp. Biol.
201
,
837
-
850
.
Kazandjian
,
T. D.
,
Petras
,
D.
,
Robinson
,
S. D.
,
van Thiel
,
J.
,
Greene
,
H. W.
,
Arbuckle
,
K.
,
Barlow
,
A.
,
Carter
,
D. A.
,
Wouters
,
R. M.
,
Whiteley
,
G.
et al. 
(
2021
).
Convergent evolution of pain-inducing defensive venom components in spitting cobras
.
Science
371
,
386
-
390
. science.abb9303
Kerkkamp
,
H. M. I.
,
Casewell
,
N. R.
and
Vonk
,
F. J.
(
2017
).
Evolution of the snake venom delivery system
. In
Evolution of Venomous Animals and Their Toxins
(ed.
P.
Gopalakrishnakone
and
A.
Malhotra
), pp.
303
-
316
.
Dordrecht
:
Springer
.
Kurt
,
J.
and
Aurich
,
H.
(
1976
).
The effect of pH value and temperature on the stability of L-aminoacidoxidase from the venom of the sand viper
.
Acta Biol. Med. Ger.
35
,
175
-
182
.
Lillywhite
,
H.
(
2014
).
Temperature and ectothermy
. In
How Snakes Work. Structure, Function and Behavior of the World's Snakes
, pp.
103
-
116
.
Oxford
:
Oxford University Press
.
Mackessy
,
S. P.
,
Sixberry
,
N. M.
,
Heyborne
,
W. H.
and
Fritts
,
T.
(
2006
).
Venom of the brown treesnake, Boiga irregularis: ontogenetic shifts and taxa-specific toxicity
.
Toxicon
47
,
537
-
548
.
McCleary
,
R. J.
,
Sridharan
,
S.
,
Dunstan
,
N. L.
,
Mirtschin
,
P. J.
and
Kini
,
R. M.
(
2016
).
Proteomic comparisons of venoms of long-term captive and recently wild-caught Eastern brown snakes (Pseudonaja textilis) indicate venom does not change due to captivity
.
J. Proteom.
144
,
51
-
62
.
Mukherjee
,
A. K.
and
Maity
,
C. R.
(
2002
).
Biochemical composition, lethality and pathophysiology of venom from two cobras Naja naja and N. kaouthia
.
Comp. Biochem. Physiol. B Biochem. Mol. Biol.
131
,
125
-
132
.
Munson
,
B. R.
,
Young
,
D. F.
,
Okiishi
,
T. H.
and
Huebsch
,
W. W.
(
2006
).
Fundamentals of Fluid Mechanics
.
New Jersey
:
John Wiley and Sons, Inc
.
Nawarak
,
J.
,
Phutrakul
,
S.
and
Chen
,
S. H.
(
2004
).
Analysis of lectin-bound glycoproteins in snake venom from the Elapidae and Viperidae families
.
J. Proteome Res.
3
,
383
-
392
.
Nelsen
,
D. R.
,
Nisani
,
Z.
,
Cooper
,
A. M.
,
Fox
,
G. A.
,
Gren
,
E. C. K.
,
Corbit
,
A. G.
and
Hayes
,
W. K.
(
2014
).
Poisons, toxungens, and venoms: redefining and classifying toxic biological secretions and the organisms that employ them
.
Biol. Rev.
89
,
450
-
465
.
Pagel
,
M.
(
1999
).
Inferring the historical patterns of biological evolution
.
Nature
401
,
877
-
884
.
Panagides
,
N.
,
Jackson
,
T. N. W.
,
Ikonomopoulou
,
M. P.
,
Arbuckle
,
K.
,
Pretzler
,
R.
,
Yang
,
D. C.
,
Ali
,
S. A.
,
Koludarov
,
I.
,
Dobson
,
J.
,
Sanker
,
B.
et al. 
(
2017
).
How the cobra got its flesh-eating venom: cytotoxicity as a defensive innovation and its co-evolution with hooding, aposematic marking, and spitting
.
Toxins
9
,
103
.
Paterna
,
A.
(
2019
).
Spitting behaviour in the Chinese cobra Naja atra
.
Herpetol. Bull.
148
,
22
-
25
.
Petras
,
D.
,
Sanz
,
L.
,
Segura
,
A.
,
Herrera
,
M.
,
Villalta
,
M.
,
Solano
,
D.
,
Vargas
,
M.
,
Leon
,
G.
,
Warrell
,
D. A.
,
Theakston
,
R. D. G.
et al. 
(
2011
).
Snake venomics of African spitting cobras: toxin composition and assessment of congeneric cross-reactivity of the pan-African EchiTAb-Plus-ICP antivenom by antivenomics and neutralization approaches
.
J. Proteome Res.
10
,
1266
-
1280
.
Rasmussen
,
S.
,
Young
,
B.
and
Krimm
,
H.
(
1995
).
On the ‘spitting’ behaviour in cobras (Serpentes: Elapidae)
.
J. Zool., Lond.
231
,
27
-
35
.
Ribeiro
,
P. H.
,
Zuliani
,
J. P.
,
Fernandes
,
C. F.
,
Calderon
,
L. A.
,
Stábeli
,
R. G.
,
Nomizo
,
A.
and
Soares
,
A. M.
(
2016
).
Mechanism of the cytotoxic effect of L-amino acid oxidase isolated from Bothrops alternatus snake venom
.
Int. J. Biol. Macromol.
92
,
329
-
337
.
Sanhajariya
,
S.
,
Duffull
,
S.
and
Isbister
,
G.
(
2018
).
Pharmacokinetics of snake venom
.
Toxins
10
,
73
.
Santra
,
V.
and
Wüster
,
W.
(
2017
).
Naja kaouthia (monocled cobra). Behavior/spitting
.
Herpetol. Rev.
48
,
455
-
456
.
Schmidt
,
J. O.
(
2019
).
Pain and lethality induced by insect stings: an exploratory and correlational study
.
Toxins
11
,
427
-
441
.
Silva-de-França
,
F.
,
Villas-Boas
,
I. M.
,
de Toledo Serrano
,
S. M.
,
Cogliati
,
B.
,
de Andrade Chudzinski
,
S. A.
,
Lopes
,
P. H.
,
Shigueo Kitano
,
E.
,
Kimori Okamoto
,
C.
and
Tambourgi
,
D. V.
(
2019
).
Naja annulifera snake: new insights into the venom components and pathogenesis of envenomation
.
PLoS Negl. Trop. Dis.
13
,
e0007017
.
Slowinski
,
J. B.
,
Knight
,
A.
and
Rooney
,
A. P.
(
1997
).
Inferring species trees from gene trees: A phylogenetic analysis of the Elapidae (Serpentes) based on the amino acid sequences of venom proteins
.
Mol. Phylogenet. Evol.
8
,
349
-
362
.
Soares
,
S. G.
and
Oliveira
,
L. L.
(
2009
).
Venom-sweet-venom: N-Linked glycosylation in snake venom toxins
.
Protein Peptide Lett.
16
,
913
-
919
.
Soares
,
T.
and
Santos
,
E.
(
2013
).
Preliminary design of fuel filling systems applying the extended Bernoulli equation on numerical calculation tools
.
SAE Tech. Pap.
13
,
1
-
6
.
Synolakis
,
C. E.
and
Badeer
,
H. S.
(
1989
).
On combining the Bernoulli and Poiseuille equation—A plea to authors of college physics texts
.
Am. J. Phys.
57
,
1013
-
1019
.
Takahashi
,
T.
and
Ohsaka
,
A.
(
1970
).
Purification and some properties of two hemorrhagic principles (HR2a and HR2b) in the venom of Trimeresurus flavoviridis; complete separation of the principles from proteolytic activity
.
BBA Protein Struct.
207
,
65
-
75
.
Tan
,
N. H.
and
Tan
,
C. S.
(
1988
).
A comparative study of cobra (Naja) venom enzymes
.
Comp. Biochem. Physiol. B Comp. Biochem.
90
,
745
-
750
.
Terry
,
A. E.
,
Knight
,
D. P.
,
Porter
,
D.
and
Vollrath
,
F.
(
2004
).
pH induced changes in the rheology of silk fibroin solution from the middle division of Bombyx mori silkworm
.
Biomacromolecules
5
,
768
-
772
.
Triep
,
M.
,
Hess
,
D.
,
Chaves
,
H.
,
Brücker
,
C.
,
Balmert
,
A.
,
Westhoff
,
G.
and
Bleckmann
,
H.
(
2013
).
3D Flow in the venom channel of a spitting cobra: do the ridges in the fangs act as fluid guide vanes?
PLoS ONE
8
,
1
-
11
.
Vitt
,
L. J.
and
Caldwell
,
J. P.
(
2013
).
Herpetology: an Introductory Biology of Amphibians and Reptiles
.
Academic press
.
Vonk
,
F. J.
,
Admiraal
,
J. F.
,
Jackson
,
K.
,
Reshef
,
R.
,
De Bakker
,
M. A. G.
,
Vanderschoot
,
K.
,
Van Den Berge
,
I.
,
Van Atten
,
M.
,
Burgerhout
,
E.
,
Beck
,
A.
et al. 
(
2008
).
Evolutionary origin and development of snake fangs
.
Nature
454
,
630
-
633
.
Ward-Smith
,
H.
,
Arbuckle
,
K.
,
Naude
,
A.
and
Wüster
,
W.
(
2020
).
Fangs for the memories? A survey of pain in snakebite patients does not support a strong role for defense in the evolution of snake venom composition
.
Toxins
12
,
201
.
Westhoff
,
G.
,
Tzschätzsch
,
K.
and
Bleckmann
,
H.
(
2005
).
The spitting behavior of two species of spitting cobras
.
J. Comp. Physiol. A Neuroethol. Sensory, Neural, Behav. Physiol.
191
,
873
-
881
.
Westhoff
,
G.
,
Boettig
,
M.
,
Bleckmann
,
H.
and
Young
,
B. A.
(
2010
).
Target tracking during venom ‘spitting’ by cobras
.
J. Exp. Biol.
213
,
1797
-
1802
.
Wüster
,
W.
(
1996
).
Taxonomic changes and toxinology: systematic revisions of the Asiatic cobras (Naja naja species complex)
.
Toxicon
34
,
399
-
406
.
Wüster
,
W.
and
Thorpe
,
R. S.
(
1992
).
Dentitional phenomena in cobras revisited: spitting and fang structure in the Asiatic species of Naja (Serpentes: Elapidae)
.
Herpetologica
48
,
424
-
434
.
Wüster
,
W.
,
Crookes
,
S.
,
Ineich
,
I.
,
Mané
,
Y.
,
Pook
,
C. E.
,
Trape
,
J. F.
and
Broadley
,
D. G.
(
2007
).
The phylogeny of cobras inferred from mitochondrial DNA sequences: evolution of venom spitting and the phylogeography of the African spitting cobras (Serpentes: Elapidae: Naja nigricollis complex)
.
Mol. Phylogenet. Evol.
45
,
437
-
453
.
Young
,
B. A.
and
Kardong
,
K. V.
(
2007
).
Mechanisms controlling venom expulsion in the western diamondback rattlesnake, Crotalus atrox
.
J. Exp. Zool. Part A Ecol. Genet. Physiol.
307A
,
18
-
27
.
Young
,
B. A.
and
O'Shea
,
M.
(
2005
).
Analyses of venom spitting in African cobras (Elapidae: Serpentes)
.
Afr. Zool.
40
,
71
-
76
.
Young
,
B. A.
,
Blair
,
M.
,
Zahn
,
K.
and
Marvin
,
J.
(
2001
).
Mechanics of venom expulsion in Crotalus, with special reference to the role of the fang sheath
.
Anat. Rec.
264
,
415
-
426
.
Young
,
B. A.
,
Dunlap
,
K.
,
Koenig
,
K.
and
Singer
,
M.
(
2004
).
The buccal buckle: the functional morphology of venom spitting in cobras
.
J. Exp. Biol.
207
,
3483
-
3494
.
Young
,
B. A.
,
Boetig
,
M.
and
Westhoff
,
G.
(
2009
).
Functional bases of the spatial dispersal of venom during cobra “spitting.” Physiol
.
Biochem. Zool.
82
,
80
-
89
.
Young
,
B. A.
,
Herzog
,
F.
,
Friedel
,
P.
,
Rammensee
,
S.
,
Bausch
,
A.
and
Van Hemmen
,
J. L.
(
2011
).
Tears of venom: hydrodynamics of reptilian envenomation
.
Phys. Rev. Lett.
106
,
198103
.
Zancolli
,
G.
,
Calvete
,
J. J.
,
Cardwell
,
M. D.
,
Greene
,
H. W.
,
Hayes
,
W. K.
,
Hegarty
,
M. J.
,
Herrmann
,
H. W.
,
Holycross
,
A. T.
,
Lannutti
,
D. I.
,
Mulley
,
J. F.
et al. 
(
2019
).
When one phenotype is not enough: divergent evolutionary trajectories govern venom variation in a widespread rattlesnake species
.
Proc. Royal Soc. Lond. B
286
,
20182735
.
Zheng
,
Y.
and
Wiens
,
J. J.
(
2016
).
Combining phylogenomic and supermatrix approaches, and a time-calibrated phylogeny for squamate reptiles (lizards and snakes) based on 52 genes and 4162 species
.
Mol. Phylogenet. Evol.
94
,
537
-
547
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information