Unexpected lack of specialisation in the flow properties of spitting cobra venom

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 A₂.

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.

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.

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⁻¹ to 10, 000 s⁻¹ (41 steps, 15 s per step), the maximum shear rate achievable by this instrument and geometry combination. Data below 100 s⁻¹ 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.

Assuming that all species have the same velocity at the channel inlet and density, Reynolds numbers are: B. arietans, Re_max.=27.62; N. nigricollis, Re_max.=23.25 and N. nivea, Re_max.=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.

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.

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.

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⁻¹ 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⁻¹ (‘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⁻¹, 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⁻¹ at the same time.

For all Naja venoms tested, the protein concentrations had an average of 132.6 mg ml⁻¹, ranging from 51.11 mg ml⁻¹ (N. nivea) to 159.1 mg ml⁻¹ (N. annulifera). The venoms of B. arietans and H. haemachatus had similar protein concentrations (132.4 and 132.5 mg ml⁻¹, 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).

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⁻¹) (Fig. 2). No significant differences between species or groups were evident (Table 1 and below).

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×10⁶ Pa), while the spitter N. nigricollis presents a lower value (0.17×10⁶ Pa). The viper B. arietans shows the lowest pressure differential (0.10×10⁶ Pa). The pressure differential results for the three snake species are reported in Table 2.

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⁻¹). 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.

Protein concentration, pH and viscosity at 10,000 s⁻¹ 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 λ).

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⁻¹ 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⁻¹. 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⁻¹ 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⁻¹ (Triep et al., 2013), and from 0.01 to 200 s⁻¹ (Young et al., 2011). In these cases, we suggest that, because of the surface tension artefacts, only data from 100 to 200 s⁻¹ (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 (A₁>A₂), 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.

There is pressure loss in fangs associated to converging diameter, which means r_{base of the fang} > r_{end 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).

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.