From birth to bite: the evolutionary ecology of India's medically most important snake venoms

Abstract

Background

Snake venoms can exhibit remarkable inter- and intraspecific variation. While diverse ecological and environmental factors are theorised to explain this variation, only a handful of studies have attempted to unravel their precise roles. This knowledge gap not only impedes our understanding of venom evolution but may also have dire consequences on snakebite treatment. To address this shortcoming, we investigated the evolutionary ecology of venoms of Russell’s viper (Daboia russelii) and spectacled cobra (Naja naja), India’s two clinically most important snakes responsible for an alarming number of human deaths and disabilities.

Methodology

Several individuals (n = 226) of D. russelii and N. naja belonging to multiple clutches (n = 9) and their mothers were maintained in captivity to source ontogenetic stage-specific venoms. Using various in vitro and in vivo assays, we assessed the significance of prey, ontogeny and sex in driving venom composition, function, and potency.

Results

Considerable ontogenetic shifts in venom profiles were observed in D. russelii, with the venoms of newborns being many times as potent as juveniles and adults against mammalian (2.3–2.5 ×) and reptilian (2–10 ×) prey. This is the first documentation of the ontogenetic shift in viperine snakes. In stark contrast, N. naja, which shares a biogeographic distribution similar to D. russelii, deployed identical biochemical cocktails across development. Furthermore, the binding kinetics of cobra venom toxins against synthetic target receptors from various prey and predators shed light on the evolutionary arms race.

Conclusions

Our findings, therefore, provide fascinating insights into the roles of ecology and life history traits in shaping snake venoms.

Background

Venom is an intrinsically ecological trait that enables venom-producing animals to incapacitate prey, outcompete conspecifics, deter predators or a combination thereof [1, 2]. As a result, venom toxins are theorised to coevolve with their interacting partners [1, 2]. Snake venoms, in particular, are well-known to exhibit considerable inter- and intraspecific variation in composition, function and potency that is largely driven by their evolutionary ecology [1]. Numerous factors, including climatic conditions, sexual dimorphism, dietary specialisation, predator pressures and ontogenetic shifts, have been invoked to explain venom variation in snakes [2, 3]. However, studies investigating the importance of these factors in shaping snake venoms have been limited. This is especially true for the medically most important snakes of the Indian subcontinent that are responsible for a vast number of human envenomings, deaths and disabilities [4, 5].

Studies in the past have highlighted the effect of intraspecific venom variation on the clinical manifestations of envenomings and the efficiency of commercial polyvalent antivenoms in treating snakebite [6,7,8,9,10,11,12,13,14]. For instance, Russell’s viper envenomings are often associated with severe local necrosis, venom-induced consumptive coagulopathy (VICC), nephrotoxicity and systemic haemorrhage [15, 16]. However, certain populations of this species have also been reported to inflict peripheral arterial thrombosis [17], neuromuscular paralysis [18] and bilateral blindness [19] in rare instances owing to the compositional variation. While the spectacled cobra envenomings are characterised by systemic neurotoxicity, local dermonecrosis and myonecrosis are also observed in certain regions [6, 20]. These variations in the clinical symptoms, combined with the lack of broadly effective antivenoms, serve as a major hurdle for the clinical management of snakebites. Therefore, understanding the evolutionary ecology of venoms is crucial for the design of a broadly effective snakebite therapy, particularly in regions that suffer the highest burden of snakebite [1].

To address this knowledge gap, we assessed the roles of various ecological and evolutionary processes in shaping the venoms of Russell’s viper (Daboia russelii) and spectacled cobra (Naja naja)—two snake species that collectively contribute the highest number of snakebite-related deaths and disabilities in India. By comparatively assessing their venom proteomes, biochemical and pharmacological functions, receptor-toxin interactions and venom potencies against a variety of prey and predatory animals, we evaluated the roles of sex, ontogeny and evolutionary arms race in shaping the venoms of these snakes. Our findings highlight, for the first time, how distinct ecological and evolutionary forces have influenced the venoms of India’s medically most important snakes.

Results

Proteomic composition of young and adult snake venoms

Electrophoretic profiles of D. russelii and N. naja venoms from various developmental stages provided interesting insights into ontogenetic venom variation (Fig. 1; Additional file 1: Fig S1–2). The venoms of D. russelii adults (mothers and unrelated males) were enriched with high-molecular-weight (HMW; 50–70 kDa) toxins, whereas the venoms of neonates and juveniles were rich in low-molecular-weight (LMW; 10–15 kDa) toxins. Prominent LMW bands (2–5 kDa) were also noted in adults but were less abundant in neonates and juveniles (Fig. 1A; Additional file 1: Fig S1). While the differences in venom sodium dodecyl sulphate–polyacrylamide gel electrophoresis (SDS-PAGE) profiles were more pronounced in D. russelii, only subtle differences were observed in the venoms of juvenile and adult N. naja (Fig. 1B). Interestingly, inter-individual differences in the intensities of bands corresponding to LMW toxins (2–10 kDa) were also noted in the venoms of juvenile individuals of N. naja sampled from the same clutch (J1–J22; Fig. 1B; Additional file 1: Fig S2). A similar pattern was observed in the venoms of neonates and juveniles collected from different clutches of D. russelii (Additional file 1: Fig S1). When the venoms of juvenile D. russelii snakes, periodically collected every 3 months, were analysed, a distinct shift in venom phenotype was recorded between 6 and 9 months (Additional file 1: Fig S1D). The venom profile of 9-month-old juveniles resembled the adult profile, with a prominent HMW band around 70 kDa. The venoms of juvenile N. naja were also collected every 3 months up to a year, but noticeable ontogenetic shifts were not observed in distinct clutches (Additional file 1: Fig S2). To understand the influence of sex, the venom profiles of unrelated adult male and female D. russelii and N. naja were compared (Additional file 1:Fig S3). However, prominent sex-specific differences in the SDS-PAGE banding patterns were not observed.

Fig. 1

SDS-PAGE venom profiles of ontogenetic stages of D. russelii and N. naja. The figures depict the SDS-PAGE profiles of (A) D. russelii and N. naja (B) venoms (replicate n = 1). The highlighted regions represent the major differences between the venom profiles. Relationships of individuals from which the venoms were collected and pooling information are provided in Additional file 2: Tables 1 and 2. M standard protein ladder, ♂ unrelated male, ♀ adult mother, N neonate, J juvenile

Additionally, D. russelii venoms were also subjected to RP-HPLC (Fig. 2). Considering the lack of sufficient amounts of venoms from the early developmental stages, and the high replicability of Reversed-phase high-performance liquid chromatography (RP-HPLC) experiments (Additional file 1: Fig S4), we generated a single fractionation profile for each venom. Here, neonate, juvenile and adult venoms separated into 12, 14 and 16 fractions, respectively. These fractions were subsequently subjected to SDS-PAGE (Additional file 1: Fig S5), followed by in-gel digestion and mass spectrometry. In RP-HPLC profiles, the prominent peaks that eluted around the 30th minute (peaks 1–3) in the adult venom were negligible in both neonate and juvenile D. russelii venoms. When we subjected these fractions from adult venoms to LC–MS/MS, we detected C-type lectins (CTL), cytotoxic three-finger toxins (C-3FTx) and disintegrins (Additional file 2: Table S2A; Additional file 4). Similarly, the other major peak around the 40th minute (peaks 4–6) in the adult venom, which contained Kunitz-type serine protease inhibitors (Kunitz) and Phospholipases (PLA₂s), was negligible in the venoms of the younger stages. Likewise, the 84th-minute peak 14 corresponding to snake venom metalloproteinase (SVMP) was also negligible in the younger stages. Venoms were separated into four fractions between 60 and 80 min for neonates (peaks 7–11) and juveniles (peaks 8–11), whereas a single prominent peak (#9) was detected in the adult venom. We recovered PLA₂, snake venom serine protease (SVSP) and vascular endothelial growth factor (VEGF) toxin families from these peaks of neonates and juveniles (Additional file 2: Table S2B and S2C; Additional file 4). The same toxin families were also recovered from the corresponding peak (peak #9) in adult venoms. Overall, in terms of the number of RP-HPLC fractions recovered, the venom of the adult D. russelii appeared to be more complex than that of neonate and juvenile individuals. The major toxin families detected in adult venoms included PLA₂ (20.4%), CTL (20.2%), SVSP (15%), Kunitz (13.3%) and SVMP (11.4%), whereas neonates and juveniles had PLA₂ (~ 39–44.2%) and SVSP (~ 39%). Minor toxin families identified in adult venoms included cysteine-rich secretory proteins (CRISP) (6%), VEGF (5.1%), C-3FTx (4.2%) and L-amino acid oxidase (LAAO; 3.9%). In the case of neonates and juveniles, minor toxin families included SVMP (~ 4.8–7.2%), CTL (~ 1%) and Kunitz (~ 1%). The abundance of neurotoxic three-finger toxins (N-3FTx), 5’-nucleotidase (5'-NT), phosphodiesterase (PDE), nerve growth factor (NGF) and hyaluronidase was less than 1% in adult venoms. Similarly, NGF, C-3FTx and 5’-NT were detected in both neonates and juveniles, and their proportions were lower than 1% of the whole venom. Additionally, PDE, hyaluronidase, natriuretic peptides (NP) and phospholipase B (PLB) were only detected in neonates in negligible amounts. Overall, the venom profile of D. russelii was consistent with previous reports and was largely dominated by PLA₂, SVMP and SVSP toxins [21,22,23,24,25,26]. Previously, N-3FTx and C-3FTx transcripts have been identified from the venom gland transcriptomes of viperid snakes [27,28,29,30]. Similarly to a few reports ([7, 31, 32]), we were also able to recover these 3FTx variants from the venom proteome of D. russelii. While it is unlikely that 3FTxs play any major role in Daboia envenoming/prey capture, their venom gland expression reflects on the recruitment of this toxin type at the base of Caenophidia and Henophidia [33].

Fig. 2

Ontogenetic shift in the venom proteomic composition of Russell’s viper. This figure depicts the A relative proteomic abundance and B RP-HPLC profiles corresponding to the ontogenetic stages (neonate, juvenile and adult) of D. russelii. The relative abundance of toxin families, determined by tandem mass spectrometry, is shown as parts of the whole graph on the left. RP-HPLC profiles (n = 1) corresponding to each developmental stage, where the absorbance at 215 nm (y-axis) is plotted against elution time (x-axis) in minutes, are shown on the right. The supporting data are provided in Additional file 2: Table S2A-S2C

In contrast to the stark ontogenetic shifts observed in D. russelii, RP-HPLC venom profiles of N. naja juveniles and their adult mother revealed lesser differences as they resolved into 14 and 12 peaks, respectively. A major difference was seen in the peak that eluted at the 44th minute (peaks 4–6) in the adult N. naja venom, where a higher peak area was documented compared to corresponding elution in juveniles (Fig. 3). Although the elution profiles of fractions between 60 and 70 min were similar in juvenile (peaks 10–12) and adult venoms (peaks 9 and 10), their proteomic composition varied (Fig. 3, Additional file 2: Table S2D and S2E; Additional file 4). The difference in the proteomic profiles was majorly due to the abundance of CRISPs in juveniles and C-3FTx in adults. While the major components in juvenile N. naja venoms were CRISPs (30%), N-3FTx (24.4%) and C-3FTx (16%), the venom of their adult mother was chiefly composed of C-3FTx (64%) and N-3FTx (26.5%). Minor toxin families documented in juveniles were Kunitz (7.4%), PLA₂ (6.3%), LAAO (6%), SVSP (4.4%) and SVMP (1.8%), whereas Kunitz (4.9%) and CRISP (1%) were identified in the adult female. LAAO, PLA₂, cobra venom factor (CVF), SVMP, CTL, PDE, acetylcholinesterase (ACE) and 5'-NT were less than 1% in adult venoms, whereas in the case of juveniles, NGF, Vespryn and CVF constituted less than 1%.

Fig. 3

Ontogenetic shift in the venom proteomic composition of spectacled cobra. This figure depicts the A relative proteomic abundance and B RP-HPLC profiles corresponding to the ontogenetic stages (juvenile and adult) of N. naja. The relative abundance of toxins, determined by tandem mass spectrometry, is shown as parts of the whole graph on the left. RP-HPLC profiles (n = 1) corresponding to each developmental stage, where the absorbance at 215 nm (y-axis) is plotted against elution time (x-axis) in minutes, are shown on the right. The supporting data are provided in Additional file 2: Table S2D-S2E

As we could only source very small amounts of venoms (5 µl) from neonates of N. naja, the RP-HPLC, mass spectrometry and other functional and toxinological comparisons could not be performed. However, marked differences were not observed in the venom SDS-PAGE profiles of neonates and juveniles of this species (Additional file 1: Fig S2).

Biochemical assessment of toxin functions

To assess the impact of compositional variation of D. russelii venom on its function, several biochemical assays were conducted. The venoms of younger stages exhibited relatively higher PLA₂ activity (~ 400 nmol/mg/min) than the venoms of adults (~ 100 nmol/mg/min; p < 0.0001; Fig. 4A). When the proteolytic activities of young and adult D. russelii venoms were evaluated, the adult venoms exhibited higher proteolysis (> 30%) compared to the younger individuals (0 to 1%; p < 0.0001; Fig. 4B). Preincubation of adult venoms with ethylenediaminetetraacetic acid (EDTA) inhibited the cleavage of azocasein, but phenylmethylsulfonyl fluoride (PMSF) treatment did not affect the venom activity. These results are in line with the relative proteomic abundance of PLA₂s and SVMPs in the venoms of younger and adult individuals. The venoms of younger D. russelii also induced higher direct haemolysis in comparison to the adult venoms that lacked quantifiable haemolytic activity (p < 0.0001; Fig. 4C). In addition, when tested for their ability to cause deamination of L-amino acid substrates, the adult D. russelii venoms exhibited relatively higher activity in comparison to the venoms of younger stages (p < 0.0001; Fig. 4D). Similar patterns were observed for all the neonate and juvenile venoms collected from various clutches (Additional file 1: Fig S6). Consistent with the SDS-PAGE profile, venoms collected from the 9-month-old juvenile D. russelii snakes exhibited functional activities similar to that of the adults. The venoms of both the younger and adult stages resulted in the cleavage of the Aɑ human fibrinogen chain, which was, in turn, inhibited by the addition of either EDTA or PMSF (Additional file 1: Fig S7A). The Clauss fibrinogen quantification assay revealed that the fibrinogen concentrations in plasma (control—251.85 mg/dl) treated with the venoms of younger stages were significantly lower than the plasma treated with the adult venoms (p < 0.0001; Fig. 4E). Interestingly, the juvenile D. russelii venom reduced the fibrinogen concentration to as low as 30 mg/dl from an initial concentration of 252 mg/dl (control plasma), while the venoms of neonate and adult stages reduced the fibrinogen concentration to 70 mg/dl and ~ 160 mg/dl, respectively (Fig. 4E).

Fig. 4

Enzymatic activities of venoms of young and adult stages of D. russelii and N. naja. This figure depicts PLA₂ (A and F), proteolytic (with and without SVMP and SVSP inhibitors) (B), haemolytic (C), LAAO (D and G), and fibrinogenolytic activities (E) of venoms of young and adult D. russelii (A–E) and N. naja (F–G). In these figures, the mean activities of replicates (n = 3) are provided, and the standard deviation is represented as error bars. NC negative control, PC positive control, N neonate, J juvenile, and ♀ adult mother. The activities of the young and adult venoms are compared using one-way ANOVA. The statistical significance is represented as follows: p < 0.01 and p < 0.0001 are indicated as ** and ****, respectively. ns indicates statistical insignificance. The supporting data are provided in Additional File 3

When the venoms of juvenile and adult N. naja were tested, their functional profiles did not exhibit prominent differences, which was consistent with their venom proteomic compositions (Figs. 4F–G). While the venoms of juvenile N. naja exhibited slightly lower PLA₂ activity (34–40 nmol/mg/min) than that of the adults (81–104 nmol/mg/min; p < 0.0001), the LAAO activity of juvenile venoms (1762–2389 nmol/mg/min) was comparatively higher than adults (775–849 nmol/mg/min; p < 0.0001). Additionally, notable differences were not observed in the fibrinogenolytic activities of N. naja venoms (Additional file 1: Fig S7B).

Acetylcholine receptor specificity of venoms

Since non-enzymatic 3FTxs constitute a major proportion of N. naja venoms, their binding kinetics against synthetic nicotinic acetylcholine receptors (nAChRs; ɑ1-subtype) from various organisms, including prey and predators of N. naja, was assessed. Here, the wavelength shift or the increase in thickness of the molecular layer on streptavidin sensor tips, resulting from the interaction between the venoms and biotinylated nAChR mimotopes, was recorded and plotted against time (Fig. 5; Additional file 1: Fig S8). The association phase is defined as the time interval during which the analyte interacts with the mimotope molecules. In contrast, the dissociation phase represents the time when the neurotoxin-receptor complexes bound to streptavidin sensor tips dissociate into the assay buffer in the absence of the analyte. The area under the curve (AUC) comparison of the association and dissociation kinetics demonstrated that Naja venoms, irrespective of whether they were from juveniles or adults, exhibited increased binding towards rat, frog and lizard nAChRs in comparison to snake and chicken receptors (p < 0.05; Additional file 1: Fig S8). However, venoms of both juvenile and adult N. naja bound the least to the receptors of mongooses and humans. Further, in the case of juvenile J2 and its adult female, the binding against the mongoose receptor was also similar to that of chicken and snake receptors (p > 0.05), with the binding against human receptors being the lowest. N. naja venoms exhibited 25% higher binding towards amphibian nAChRs than other receptors, while they exhibited 66% lower binding towards mongoose nAChRs (Fig. 5; Additional file 1: Fig S8).

Fig. 5

Binding kinetics of N. naja venoms towards ɑ1-nAChRs. In this figure, the binding kinetics of N. naja venoms from young and adult stages are depicted as line graphs (replicate n = 1). Here, the thickness of the molecular layer corresponding to the binding strength of venom toxins to ɑ1-nAChRs from various organisms, including prey and predators of N. naja, is plotted against time (s). The dotted line at 400 s represents the demarcation between the association and dissociation phases. J₁ and J₂ juvenile venoms, ♂ unrelated male, ♀ adult mother of the juveniles. The AUC comparisons for association and dissociation kinetics are provided in Additional file 1: Fig S8. The supporting data are provided in Additional File 3

Lethal potencies of venoms across developmental stages

The median lethal doses (LD₅₀) of D. russelii venoms varied across their developmental stages (Fig. 6A; Additional File 2: Table S5A). The venoms of adult (0.166 mg/kg) and juvenile (0.176 mg/kg) D. russelii from Karnataka exhibited potencies similar to that of adults from other populations across India (0.1–0.3 mg/kg; [7, 34]). However, the venoms of neonates exhibited extremely high venom potency against mice (LD₅₀ of 0.072 mg/kg), which was comparable to that of certain krait venoms [5, 10, 35]. In contrast, venoms of N. naja juvenile (0.380 mg/kg) and adult (0.363 mg/kg) stages had a venom potency that was similar to the pan-Indian N. naja populations Fig. 6A; Additional File 2: Table S6A; [6, 9, 36]. Interestingly, a similar trend was observed in assays involving house geckos, where neonate D. russelii venoms were found to be extremely potent (2.024 mg/kg). The venom of juveniles killed 50% of the test population at ~ 4 mg/kg. In both cases, geckos were paralysed and lost the righting reflex within 30 min of venom injection. On the contrary, though the adult Daboia venoms induced partial paralysis during the initial hours of experimentation (2–4 h), venom-injected geckos eventually recovered completely, even at doses as high as 20 mg/kg (Fig. 6B; Additional File 2: Table S5B). Not surprisingly, N. naja venoms exhibited higher lethal potencies against geckos in comparison to D. russelii. LD₅₀ of 0.510 and 1.265 mg/kg were recorded for adult and juvenile venoms, respectively (Fig. 6B; Additional File 2: Table S6B). However, unlike ontogenetic shifts in venom potency documented in D. russelii, a significant shift was not observed in N. naja. When the D. russelii venoms were tested on crickets, neonate and juvenile venoms did not exhibit any toxicity at concentrations as high as 16 mg/kg (Fig. 6C; Additional File 2: Table S5C). Surprisingly, though the adult venom did not induce instant paralysis, it was lethal to crickets at a concentration of 6.375 mg/kg several hours post-venom injection. When similar experiments were conducted with scorpion venoms, crickets were paralysed within 5 min of injection and lethality was observed at doses as low as 2.828 mg/kg (Fig. 6C; Additional File 2: Table S5C). Both adult and juvenile N. naja venoms did not induce paralysis in crickets, and only the adult venoms killed crickets (LD₅₀ of 5.657 mg/kg) several hours post-venom-injection (Fig. 6C; Additional File 2: Table S6C).

Fig. 6

Toxicity profiles of ontogenetic stages of D. russelii and N. naja. Bar graphs in this figure represent the LD₅₀ (mg/kg) of the young and adult D. russelii and N. naja venoms against A mouse, B gecko and C cricket. N neonate, J juvenile, ♂ unrelated male, ♀ adult mother, PC positive control (H. tamulus venom); asterisk indicates absence of lethality at the highest tested concentration. The error bars indicate 95% confidence interval (CI) calculated using Probit statistics. Representative photographs of neonate D. russelii and N. naja are also shown. Detailed information on venom doses (µg), number of test animals, survival patterns and LD₅₀ with 95% CI are provided in Additional file 2: Tables S5 and S6

Ontogenetic variation in haemorrhagic potential

Given the ontogenetic variation in the proteomic composition and function of D. russelii venoms, the variation in haemorrhagic potential was estimated in the murine model. The mean diameter of the haemorrhagic lesion induced was plotted against various concentrations of venoms tested and simple linear regression was performed. The slopes of the curves obtained for the neonate, juvenile and adult venoms were compared. The slopes were significantly different from one another (p = 0.0096). The minimum amount of venom required to induce a 10-mm haemorrhagic lesion was interpolated from these curves. While only 1.23 μg of the adult D. russelii venom was sufficient to induce a haemorrhagic lesion of 10 mm, the MHD of the neonate (5 μg/mouse) and juvenile (6.5 μg/mouse) venoms was four to five times higher (Additional file 1: Fig S9 and S10; Additional file 2: Table S7).

In vitro binding potential of antivenoms

Indirect enzyme-linked immunosorbent assays (ELISAs) were performed to assess the binding efficacy of various Indian antivenoms against the venoms of adult and younger individuals. Here, the titre values of various antivenoms were calculated by comparing them against the naive horse IgG negative control. These titre values indicate the highest dilution at which the antivenoms exhibit specificity towards the target venom. The negative control cutoff was also used as the baseline the half maximal inhibitory concentration (IC₅₀) estimates. Amongst the four polyvalent antivenom products tested (Additional File 2: Table S4), VINS and Premium Serums performed better against a majority of D. russelii and N. naja venoms as reported previously [6, 7]. These antivenoms also recognised neonate and adult D. russelii venoms similarly, with a titre of 1:12,500 (IC₅₀: 0.6–1 µg/mL; p < 0.0001) but exhibited a relatively lower titre towards juvenile venoms [(1:2500); IC₅₀: 4 and 3 µg/mL for VINS and Premium Serums, respectively; p < 0.0015; Fig. 7; Additional file 1: Fig S11A]. However, in the case of N. naja, antivenoms exhibited slightly higher binding towards the venoms of juveniles and adult females (1:12,500) in comparison to the venom of adult males [a titre of 1:2500, as reported previously [9]]. Though negligible differences were observed in the IC₅₀ of various antivenoms towards the venom of N. naja juvenile 02 (p = 0.1240), the Bharat Serums and Haffkine antivenoms exhibited significantly lower binding towards the juvenile 01 with IC₅₀ estimates of 3.5 and 4.2 µg/mL, respectively (p < 0.0001; Additional file 1: Fig S11B). Overall, IC₅₀ estimates also highlighted that the Premium Serums antivenom exhibited better recognition potential against the N. naja and D. russelii venoms (Additional file 1: Fig S11). Therefore, based on the titre values and the IC₅₀ estimates, the best binding antivenom (Premium Serums) was further downselected for animal protection studies. This downselection was performed with ethical considerations to reduce the total number of mice sacrificed.

Fig. 7

In vitro binding of antivenoms against the venoms of various ontogenetic stages of D. russelii and N. naja. Line graphs in this figure depict the in vitro binding potentials of Indian polyvalent antivenoms against the young and adult [1,2,3,4] D. russelii and [5,6,7,8] N. naja venoms (replicate n = 3). The binding was assessed by plotting absorbance at 405 nm against various dilutions of the antivenom (x-axis) in indirect ELISA. Dotted vertical lines indicate the titre values of respective antivenoms calculated using the negative control cutoff. B Bharat Serums and Vaccines Ltd, H Haffkine Bio-Pharmaceutical Corporation Ltd, P Premium Serums & Vaccines Pvt. Ltd, V VINS Bioproducts Ltd. The following symbols represent various developmental stages: N neonate, J juvenile, ♂ unrelated male, ♀ adult mother. The IC₅₀ estimates are provided in the Additional file: Fig S11. The supporting data are provided in Additional File 3

In vivo neutralisation potency of polyvalent antivenom

While the in vitro binding assay suggested that Indian antivenoms recognise the venoms of younger and adult individuals with similar efficiency, the in vivo mice toxicity assays revealed the shortcomings of the marketed antivenoms. In line with previous findings, the Premium Serums antivenom, which, despite being one of the best binding antivenoms under in vitro conditions in this study, exhibited lower neutralisation potencies in comparison to the market claim against all tested venoms (Fig. 8). While the marketed neutralisation potency for D. russelii and N. naja venoms is 0.60 mg/mL, the neutralisation potencies towards the venoms of younger and adult stages of D. russelii from Karnataka ranged between 0.319 mg/mL and 0.468 mg/mL (Additional File 2: Table S8A). Similarly, the neutralisation potency of this antivenom against the juvenile (0.221 mg/mL) and adult (0.321 mg/mL) N. naja venom was critically lower than the marketed efficiency (Additional File 2: Table S8B).

Fig. 8

Neutralisation potencies of commercial polyvalent antivenom against D. russelii and N. naja venoms from different ontogenetic stages. The bar graphs show the neutralisation potencies of polyvalent antivenom manufactured by Premium Serums. In these graphs, the x-axis shows different samples—N neonate, J juvenile 01, ♀ adult mother of D. russelii and N. naja venoms. The y-axis represents neutralisation potency calculated using the LD₅₀ of the venoms and the effective dose (ED₅₀) corresponding to that venom in milligrammes per millilitres. The dotted lines parallel to the x-axis at 0.6 mg/mL indicate the marketed neutralisation potency of antivenom, and the error bars indicate 95% CI calculated using Probit statistics. Detailed information on the antivenom dose (µl), the survival patterns and ED₅₀ values with 95% CI are provided in Additional File 2: Table S8

Discussion

Ontogenetic shift in Russell’s viper venoms

Venom is a highly adaptive trait that is crucial for prey capture, defence, competition or a combination thereof [1]. It is, therefore, not surprising that stark inter and intraspecific differences in snake venom composition, activity and toxicity have been reported in various venomous organisms [6,7,8, 37,38,39,40,41]. Several proteomic and toxinological studies conducted across the geographic distribution of these species have revealed remarkable variations in the abundances of toxins that inflict severe neurotoxicity and coagulopathies [21,22,23,24,25,26]. In a complete contrast, a recent study has reported the role of certain factors (e.g. dietary specialisation) in promoting the conservation of venom phenotypes across wide geographic distributions [42]. Moreover, various ecological, environmental and life history traits, such as sex, season, and ontogeny, may shape venoms [43,44,45,46,47]. For instance, given the gradual increase in gape size associated with the developmental stage of snakes, variations in diet has been reported [48, 49]. In turn, to facilitate such shifts in diet, ontogenetic shifts in venoms have also been reported [50, 51]. For example, in the case of certain Crotalus spp. (e.g. Crotalus simus simus, C. molossus nigrescens, C. oreganus helleri and C. o. oreganus), the relative abundance of venom toxins shifts from neurotoxic PLA₂s in juveniles to haemotoxic SVMPs in adults to facilitate a corresponding shift in diet from lizards to mammals [52,53,54,55]. Thus, developmental stage-specific venom compositions can underpin ontogenetic shifts in diet. However, they remain poorly investigated in snakes around the world.

The influence of ecology and environment in shaping the venoms of Indian snakes, in particular, remains elusive. For instance, while considerable geographical venom variation is well-documented in D. russelii and N. naja [6, 7, 12, 56, 57], the underlying factors responsible for this variation remain unclear. Similarly, ontogenetic shifts in the venoms of Indian snakes have never been investigated. Thus, to address these shortcomings, and to understand the ontogenetic venom variation in Indian snakes, we investigated 194 (from 7 clutches) and 32 (2 from clutches) individuals of D. russelii and N. naja, respectively. Neonates and juveniles examined in this study were captive-born offsprings of adult females caught in the wild. All adult females and males were caught from the same geographic location (Hunsur, Karnataka, India).

Our findings provide interesting insights into ontogenetic shifts, or lack thereof, in the venoms of D. russelii and N. naja—India’s two clinically most important snake species. Proteomics experiments revealed that the venoms of neonate and juvenile D. russelii were predominantly composed of PLA₂ and SVSP. However, the adult venoms transitioned into a more complex repertoire, constituted by CTL, Kunitz and SVMP in nearly equal proportions (Fig. 2; Additional file 2: Tables S2A-C). Enzyme kinetic assays supported these observations, wherein caseinolytic activity, which is largely driven by venom proteases, was recorded only for the venoms of adults (Fig. 4B). Repetition of this assay with EDTA and PMSF inhibitors confirmed that the proteolytic activity of adult venom was driven by SVMP and not by SVSP. While PLA₂s were abundant in the venoms of both younger and adult stages, higher PLA₂ activity was recorded in the former (Fig. 4A; Additional file 2: Tables S2A-C). This could be a result of a relatively higher proportion of catalytically active PLA₂s in the venoms of younger stages. Previous reports have highlighted lower PLA₂ activities in adults of certain populations of D. russelii across India [7, 34]. Therefore, it may be interesting to investigate whether adults in these populations lose catalytic PLA₂s due to ontogenetic shifts. Further differences in venom functions were offered by direct haemolysis and haemorrhagic dose experiments. Certain snake venom PLA₂s are known to disrupt the phospholipid bilayer of erythrocytes, resulting in haemolysis [58]. On the contrary, SVMPs are predominantly attributed to the haemorrhagic manifestations observed in viper envenomings [59]. Consistent with the proteomic and functional profile of D. russelii venoms, only neonate and juvenile venoms exhibited haemolysis (Fig. 4C), while adults induced haemorrhage at venom concentrations lower than that of younger stages (1.23 µg vs 5–6.5 µg; Additional file 1: Fig S9-S10; Additional File 2: Table S7).

D. russelii venoms are dominated by toxins that inflict intravascular coagulopathy and systemic haemorrhage. SVMPs and SVSPs in Daboia venoms activate blood coagulation factor X and factor V, respectively [60, 61]. In addition, SVMPs can induce local and systemic haemorrhage, damage blood vessels and result in profuse internal bleeding and, ultimately, death of bite victims [62]. While these manifestations help in the rapid immobilisation of mammalian prey, human snakebite victims suffer VICC and immutable injuries [15]. Venoms of all developmental stages of Daboia exhibited varying levels of fibrinogenolysis, with the neonate and juvenile individuals reducing the plasma fibrinogen concentration significantly more than the adults (30–70 mg/dl vs 160 mg/dl; p < 0.0001; Fig. 4E). Such differences could be attributed to the higher proportions of SVSPs in the venoms of younger stages (Fig. 4B; [43]). LAAO, a toxin attributed to platelet aggregation, apoptosis and necrosis in human bite victims [63], was also found in higher proportions in adult venoms (2352–2851 nmol/mg/min) compared to neonates and juveniles (52–147 nmol/mg/min; p < 0.05; Fig. 2).

Further evidence of a notable ontogenetic shift in venom phenotypes of D. russelii was provided by toxicity assays. When LD₅₀ experiments were conducted on rodents, the venoms of D. russelii neonates were 2.3 and 2.5 times more potent than their adult mother and juvenile counterparts, respectively. When the venom potency was assessed against lizards, adult Daboia did not exhibit toxicity, even up to a dose of 20 mg/kg, whereas neonates and juveniles induced lethality at doses as low as 2–4 mg/kg. We further assessed ontogenetic shifts in toxicities towards insects using crickets (0.5–1 g) as a model system. Here, the venoms of neonates and juveniles did not exhibit toxicity towards insects, even up to doses equivalent to 16 mg/kg (Fig. 6C). In contrast, the venom of adult Daboia eventually killed crickets at a test dose of 6 mg/kg, but the observed toxicity was, perhaps, not intended. Experimental insects injected with these venoms did not exhibit any signs of toxicity but ultimately succumbed to the effects of envenoming 24 h post venom injection. In contrast, the venom of the Indian red scorpion (Hottentotta tamulus)—a major insect predator—killed crickets within a minute of injection at concentrations as low as 2.8 mg/kg. These observations suggest that the venoms of D. russelii across all developmental stages may not play a role in capturing arthropod prey (Fig. 6C, Additional file 2: Tables S5C and S6C).

In summary, a shift from a higher phospholipolytic neonate and juvenile venoms to a higher proteolytic adult venom was documented in D. russelii. Periodic evaluation of venoms of various developmental stages of D. russelii in this study reveals that this shift in venom phenotype occurs between 6 and 9 months (Additional file 1: Fig S1 and S4). A similar ontogenetic shift in the venoms of Crotalus spp. (Crotalinae) has been reported to occur gradually over a period of 2 years [64]. Interestingly, this is the first documentation of ontogenetic shift in venom within the Viperinae clade.

A lack of ontogenetic shift in spectacled cobra venoms

In contrast to the considerable ontogenetic shift in venoms of D. russelii, developmental stages of N. naja were found to produce similarly functional venoms. Except for the marked variation in the proportion of C-3FTx and CRISPs between juvenile and adult venoms (Fig. 3; Additional file 2: Tables S2D-E), the overall composition, biochemical activity and potency of the venom did not significantly change across developmental stages (Figs. 5 and 6). Differences were not observed in the relative abundance of N-3FTx—a major Naja toxin superfamily—between the venoms of juveniles and adults (24.4% and 26.5%, respectively; Fig. 3). In vitro binding assays involving N. naja venoms and nAChR mimotopes from various animals, including prey and predators of this species, corroborated these findings. Both N. naja juvenile and adult venoms exhibited increased binding towards reptilian, amphibian and rodent nAChRs (Fig. 5). This is consistent with field documentation of N. naja feeding on these animals [65]. Lower binding was observed towards mongoose, snake and human mimotopes. While the ɑ1 nAChR of the Egyptian mongoose has been previously identified to be resistant to ɑ-neurotoxins [66], the lower binding towards snake receptors could be due to substitutions that prevent autotoxicity [67]. Similarly, relatively weaker binding towards human nAChRs is consistent with the fact that we are neither their prey nor predator. In agreement with the proteomic composition, the venoms of juveniles exhibited higher LAAO activity (775–849 nmol/mg/min) than adults (1762–2389 nmol/mg/min; Fig. 4F). Further, despite the presence of higher amounts of PLA₂s in juveniles (6.3% vs 0.25%; Fig. 3), adults exhibited PLA₂ activity (81–104 nmol/mg/min) that was 2 × higher than the former (34–44 nmol/mg/min; p < 0.05; Fig. 4G). However, haemolytic activity was not observed in either of these venoms. In toxicity assays, N. naja juvenile and adult venoms did not exhibit a marked variation in lethal potencies against mammalian (0.380 mg/kg vs 0.363 mg/kg) and reptilian (0.510 mg/kg vs 1.265 mg/kg) prey. Similarly to the venom of D. russelii, N. naja venom, too, did not exhibit toxicity towards crickets, suggesting that they do not target arthropods. While ontogenetic shifts were not documented in N. naja, inter-individual venom variation was observed in juveniles (n = 15; Additional file 1: Fig S2). Notable differences in abundances of LMW toxins, such as 3FTxs (< 15 kDa), were recorded in their SDS-PAGE profiles. Such differences observed in closely related individuals may contribute to considerable intrapopulation venom diversity documented in snakes [8].

Many ways to skin a cat: deployment of broadly effective versus prey-specific toxins across developmental stages

Differing predatory strategies of D. russelii and N. naja may have shaped their venom compositions and potencies across developmental stages. Spectacled cobra is known to be an opportunistic feeder that actively hunts various prey animals, including lizards, frogs, toads, snakes and rodents [65]. In contrast, Russell’s viper is predominantly an ambush predator that probably hunts specific prey at various developmental stages. For instance, we documented caudal luring by neonates in the presence of the common Indian cricket frog (Minervarya agricola). During this behaviour, neonates wiggled their tails to mimic worms—a trait that, perhaps, helps them attract prey (Additional file 5). Consistently, D. russelii neonate individuals readily accepted frogs in captivity. However, adults have not been documented to exhibit such a behaviour. It is also very clear that gape size limitations prevent neonate D. russelii from feeding on larger rodents, further pointing towards a diet shift with age. Consistent with this hypothesis, D. russelii was documented to produce distinct venom cocktails across developmental stages that exhibited varying prey specificities. Toxicity experiments revealed that neonate venoms exhibit increased potency towards mice, in comparison to juveniles and adults (Fig. 6A; Additional File 2: Table S5A). However, only neonate and juvenile venoms were lethal to lizards, and adult venoms did not induce lethality even up to a dose of 20 mg/kg (Fig. 6B; Additional File 2: Table S5B).

Several hypotheses could be put forward to explain the observed temporal variation in composition and potency of D. russelii venoms. Ontogenetic shifts in venom could relate to the developmental stage-specific changes in the gape size of snakes [48, 49]. Given the limited gape size, the younger stages may be only able to capture smaller reptilian or amphibian prey. Eventually, as adults, with an increased gape size, they could shift to feeding on larger mammals. Field observations in crotaline snakes have also identified a diet shift from ectothermic (arthropod, reptilian or amphibian) to endothermic (mammalian) prey with age, which is underpinned by a corresponding shift in venom profiles [50, 55, 68, 69]. Behavioural observations in this study point to a similar shift in diet from ectotherms (lizards and frogs) to endotherms (rodents) in D. russelii. Here, the HMW SVMPs have been implicated in the rapid immobilisation and digestion of large mammalian prey [70, 71], which could explain the shift towards an SVMP-rich venom in adult D. russelii. Moreover, an increase in the complexity of venom with age has been previously reported in several species from the Crotalinae clade [43, 47, 72]. Consistent with this, we find that the venoms of neonate and juvenile D. russelii transition into a more complex phenotype in adults, where additional toxin classes, including Kunitz and CTL, are expressed. However, their roles in adult venoms remain unclear and require further investigation.

Another hypothesis that could explain ontogenetic shifts in the composition and potency of D. russelii venom pertains to resource optimisation. We find that the venoms of neonate and juvenile D. russelii are rich in LMW toxins. While adults produce both LMW and HMW toxins, their LMW components do not appear to be the same as the ones expressed by younger stages (Fig. 2). This, perhaps, suggests that the younger stages produce certain LMW toxins that are responsible for their increased venom potency. Since younger stages also produce relatively smaller amounts of venom, the increased venom potency may compensate for the lack of sufficient yield. Moreover, the expression of highly potent LMW toxins could be metabolically more energy efficient for the newborns, as opposed to the production of large multi-domain SVMPs.

In contrast to stark shifts in Daboia venoms, we find that the venoms of N. naja are rich in 3FTxs across developmental stages (Fig. 3). Similarly, monocled cobras (N. kaouthia) have also been reported to produce 3FTx-rich venoms across their life stages [73]. Structurally and functionally diverse 3FTxs in elapid venoms can exhibit a broad range of prey specificity [74, 75]. In support of this, the results of nAChR binding and toxicity experiments in this study suggest that N. naja 3FTxs may facilitate the effective capture of various prey species, including rodents, lizards and amphibians (Figs. 5 and 6). A lack of ontogenetic shift in the venoms of elapids has also been reported in coastal taipans [76]. However, age-associated shifts in biochemical, coagulopathic and toxicity profiles have been documented in others (N. atra, N. nigricolis and Pesudonaja spp.) [77,78,79].

Studies on ontogenetic shifts in venoms have often invoked paedomorphism—retention of newborn characters into adulthood—to explain similarities between the venoms of younger and adult stages [80, 81]. This theory, first proposed by Walter Garstang in 1922, suggests that competition and predation pressure drive paedomorphism [82]. While contradictory theories have been proposed to explain the evolutionary role of this trait [83], the absence of an ontogenetic shift in N. naja venoms may be associated with phylogenetic inertia and not paedomorphy. The expression of diverse 3FTxs in venoms is an ancestral trait observed across the Elapidae clade and, hence, the increased abundance of this toxin family in both juvenile and adult venoms cannot be termed paedomorphy.

The impact of ontogenetic shifts in venom on the effectiveness of snakebite therapeutics

Considerable differences in venom compositions have been shown to reduce the efficiency of marketed antivenoms [39, 81, 84]. Despite the remarkable inter- and intraspecific venom variation, the Indian antivenoms are manufactured by sourcing venom from a point location in southern India and by using a century-old manufacturing strategy. This polyvalent antivenom manufactured against the ‘big four’ snake venoms has been shown to be ineffective/poorly effective in countering bites from the pan-Indian ‘big four’ snake populations [6,7,8,9]. In addition to geographic variation, studies have also suggested that the ontogenetic shift in venoms may also affect the efficacy of antivenoms [85]. While considerable research has been conducted towards understanding geographic variation in Indian snake venoms and its repercussion on the effectiveness of antivenoms, ontogenetic shifts in venom profiles and their clinical impact have never been documented.

In this study, in vitro binding and in vivo neutralisation potencies of commercial antivenoms from leading Indian manufacturers against the venoms of the developmental stages of D. russelii and N. naja were assessed. While the in vitro binding assay suggested that these antivenoms recognised the venoms of younger and adult snakes with similarly decent efficiencies (Fig. 7), in vivo neutralisation experiments involving mice revealed their shortcomings (Fig. 8; Additional file 2: Table S8). Our results showed that Premium Serums antivenom, which was one of the best-binding antivenoms under in vitro conditions, exhibited poor neutralisation against the venoms of all developmental stages of both D. russelii and N. naja. However, whether the observed lack of antivenom efficacy results from ontogenetic shifts in venoms remains unclear. Snake venoms investigated in this study were collected in Hunsur, Karnataka, a high-altitude Deccan Plateau zone, which is starkly different from the coastal regions of Tamil Nadu, where the venoms are sourced for the commercial production antivenoms. Therefore, the observed lack of neutralisation could also stem from the geographic variation in venom. Assessing the neutralisation potencies of Indian antivenoms on the venoms of the developmental stages of snakes from Tamil Nadu will reveal the impact of ontogenetic shifts on snakebite therapeutics. In any case, even when such an impact is documented, understanding the burden of snakebite incidents by younger individuals is imperative for suggesting policy changes for the inclusion of newborn snake venoms in the immunisation mixture.

However, based on these results, it is clear that the biogeographic venom variation plays a greater role in determining the neutralisation potency and, hence, requires immediate attention. Previous studies have also emphasised the poor performance of Indian antivenoms against the venoms from several N. naja populations, including Karnataka [6, 9]. While we note a poor neutralisation potency against D. russelii venom from Karnataka, Indian antivenoms have been reported to sufficiently neutralise the lethal effects of Daboia venoms of several other populations [7]. The only other exception to this was the northern semi-arid D. russelii population [7]. Overall, our findings collectively emphasise the need for the development of region-specific antivenoms that target the venoms of medically important snakes by the region.

Limitations

Due to sampling restrictions by the forest department and insufficient venom yield from the neonate and juvenile Russell’s vipers, technical replicates for proteomic analyses, such as RP-HPLC and LC–MS/MS, were not performed. This limitation also hindered our ability to investigate the variations in minor enzymatic components, including 5’-NTD and PDE. While these toxins could contribute to the overall compositional and functional variations between ontogenetic stages, given their limited abundance (< 1%) in the venom proteome, it is unlikely that they play an important ecological or clinical role. This limitation further prevented us from performing comparative nAChR binding assays with D. russelii venoms. While we do not expect to find any nAChR-specific toxins in D. russelii venom, these assays may shed light on the mild neurotoxic symptoms associated with envenoming from this species in certain parts of India. Additionally, we chose to analyse venoms from a single population in Karnataka to minimise the effects of regional venom variation. While we have documented identical venom profiles for neonates, juveniles and adults from other distant populations, dedicated research can elucidate the extent of ontogenetic venom variation across India. Moreover, our study did not distinguish whether the observed antivenom inefficacy is attributable to regional variation or ontogenetic differences. Future investigations should encompass venom samples from diverse ontogenetic stages collected, particularly from regions where acceptable antivenom efficacies have been demonstrated. These studies are essential to comprehensively understand the ontogenetic and regional variations in venom composition and their implications for antivenom development.

Conclusion

In conclusion, our study provides compelling evidence that ecological and life-history factors significantly influence the composition, function and potency of D. russelii and N. naja venoms across their developmental stages. While marked ontogenetic shifts in venom profiles were observed in D. russelii, the venoms of N. naja demonstrated a lack of functional variation across developmental stages. Moreover, the receptor-binding assays revealed that the venoms of N. naja exhibit a broad range of specificities towards various prey and predatory animals, thereby explaining the lack of ontogenetic shifts in the venom of this species. These findings unravel the complex interplay between ecological and evolutionary processes that shape venom diversity across different life stages of an individual. Understanding these dynamics not only enhances our knowledge of venom evolution but also has important implications for the development of more effective antivenoms and medical treatments for snakebite victims.

Methods

Venoms evaluated in this study

Venoms of N. naja and D. russelii individuals of various developmental stages, including neonates (< 30 days), juveniles (between 1 to 12 months) and mature individuals (> 36 months), were sampled with permission from the state forest department of Karnataka (PCCF(WL)/E2/CR-06/2018–19 and PCCF(W1)/C1(C3)/CR-09/2017–18); 226 individuals from 9 clutches were examined. Of these, 194 individuals belonged to seven clutches of D. russelii, and 32 belonged to two clutches of N. naja. Venom was sourced every 3 months from N. naja and D. russelii individuals to track ontogenetic changes across time. Freshly collected venoms were flash-frozen in liquid nitrogen, lyophilised and stored at − 80° C. A detailed list of venom samples analysed in this study is provided in Additional file 2: Tables S1A-B. Adult male and female snakes rescued from the Hunsur City in the Mysore district of Karnataka were housed in the same locality, and their venoms were collected periodically. Gravid females gave birth to the neonates used in the study under captive conditions. This strategy was adopted to negate the effect of intrapopulation venom variation on our findings. As venom could exhibit significant intrapopulation venom variation, even at shorter geographic scales [8], the distribution of the samples was restricted to a few kilometres around the housing site. Throughout the experimental time, all snakes were raised on a rodent diet.

Reversed-phase high-performance liquid chromatography (RP-HPLC)

Snake venoms were fractionated on a Shimadzu LC-20AD series HPLC system (Kyoto, Japan), following a modified version of the previously published protocol [12, 86]; 1 mg of venom diluted in water was injected into a reversed-phase C18 Shim-Pack GIST column [4.6 × 250 mm, 5 μm particle size and 100 Å pore size (Shimadzu, Japan)], previously equilibrated with 95% buffer A [0.1% trifluoroacetic acid (TFA) in water (v/v)] and 5% buffer B [0.1% TFA in acetonitrile (v/v)]. Venom fractions were eluted at a steady flow rate of 1 mL/min under a gradient of buffer A (0.1% TFA in HPLC grade water), buffer B (0.1% TFA in acetonitrile) and isocratic (5% buffer B) for 5 min, followed by 5–15% buffer B for 10 min, 15–45% buffer B for 60 min and 45–70% buffer B for 10 min and held at 70% buffer B for 9 min, followed by 70–98% buffer B for 6 min and again held at 98% buffer B for 5 min and completed with 5% buffer B for 5 min. The column was equilibrated at 5% buffer B for 10 min, completing a total run time of 120 min. Fractions were collected separately by monitoring the absorbance peaks at 215 nm, and the areas under the peaks were estimated to calculate their relative abundances.

Protein estimation and gel electrophoresis

The above fractions were concentrated in a vacuum centrifuge (Thermofisher Scientific, USA), followed by resuspension in ultrapure water. These fractions were subjected to 12.5% SDS-PAGE under reducing conditions at a constant voltage of 80 V [87]. Following this, the gels were stained overnight with Coomassie Brilliant Blue R-250 (Sisco Research Laboratories Pvt. Ltd, India) and visualised using an iBright CL1000 gel documentation system (Thermo Fisher Scientific, MA, USA). The densitometry analysis of individual bands was then performed using the ImageJ software [88] and the bands were excised for mass spectrometric analyses.

In-solution and in-gel trypsin digestion

A combined analysis strategy was adopted for mass spectrometry, wherein RP-HPLC fractions with higher peak area were subjected to reducing SDS-PAGE and smaller peaks that could not be visualised on the gel were directly subjected to in-solution digestion. The excised gel bands were destained using 30% acetonitrile in 50 mM ammonium bicarbonate buffer. The gel pieces were dehydrated using 100% acetonitrile, followed by reduction using 10 mM dithiothreitol (DTT) at 56 °C for 45 min and alkylation using 55 mM iodoacetamide for 35 min at room temperature. The samples were then digested using sequencing-grade trypsin and incubated overnight at 37 °C. The following day, the supernatant was collected and desalted in a Pierce C18 spin column (Thermofisher Scientific, USA) following the manufacturer’s protocol before being subjected to mass spectrometric analyses. Similarly, for in-solution digestion, fractions were reduced and alkylated using DTT (100 mM) and iodoacetamide (100 mM), respectively, followed by tryptic digestion.

Liquid chromatography-tandem mass spectrometry (LC–MS/MS)

The digested peptides were subjected to nano-liquid chromatography (nano-LC) and electrospray ionisation tandem mass spectrometry (ESI–MS/MS). The samples were injected into a PepMap C18 nano-LC column (50 cm × 75 μm, 2 μm particle size and 100 Å pore size) mounted on the Thermo EASY nLC Ultimate 3000 series system (Thermo Fisher Scientific, MA, USA), following a gradient elution of buffer A (0.1% formic acid in MS grade water) and buffer B (0.1% formic acid in 80% acetonitrile) at a constant flow rate of 250 nL/min for 90 min. An 8–35% gradient of buffer B was used over the first 70 min for elution, followed by 35–95% over the next 5 min, and finally 95% for the last 15 min. Subsequently, the fractions from the nano-LC were injected into a Thermo Orbitrap Fusion Mass Spectrometer (Thermo Fisher Scientific, MA, USA). The following protocol was adopted to perform the MS scans: scan range (m/z) of 300–2000 with a resolution of 120,000 and maximum injection time of 100 ms. To perform the precursor (MS) and fragment (MS/MS) scans, an orbitrap detector with high collision dissociation (HCD) fragmentation (30%) was used with the following parameters: scan range (m/z) of 110–2000 and maximum injection time of 50 ms [89].

Raw MS/MS spectra were searched against the National Center for Biotechnology Information non-redundant (NCBI-NR) Serpentes database (taxid: 8570; with 549,650 entries as of September 2023) to identify toxins that constitute the venom fractions. The search was performed in PEAKS Studio X Plus software (Bioinformatics Solutions Inc., ON, Canada) by setting parent and fragment mass error tolerance limits to 10 ppm and 0.6 Da, respectively. A ‘monoisotopic’ precursor ion search type and semispecific trypsin digestion were specified as parameters. Fixed and variable modifications were set as carbamidomethylation (+ 57.02) and oxidation (+ 14.99), respectively. The false discovery rate (FDR) was set to 0.1%, and PEAKS Studio automatically determined the corresponding − 10lgP cutoff value. Hits with at least one unique matching peptide were considered for protein identification, and redundant hits from each protein family were manually removed.

The proportions of toxins in the venom were determined by implementing a multipronged normalisation strategy [12, 89,90,91,92]. First, the area under the curve of spectral intensities of a toxin hit (AUC_t) was normalised within the fraction using the formula below.

$${(\text{AUC})}_{\text{MS}}= \frac{{\text{AUC}}_{ t}}{\sum_{t=1}^{N}{(\textrm{AUC})}_{t}}$$

Here, AUC_MS represents the normalised AUC_t, while ‘t’ represents toxins in a given fraction and N denotes the total number of toxins in that fraction.

To calculate the proportion of the ith toxin (T_i) in the venom, the obtained AUC_MS was further normalised by multiplying with the relative RP-HPLC area (AUC_HPLC) of the respective fraction and the density (D) of the SDS-PAGE band.

$${T}_{i}= \sum_{f=1}^{N}{({\text{AUC}}_{\text{HPLC}} \times D \times {\text{AUC}}_{\text{MS}})}_{f})$$

Here, ‘N’ denotes the total number of fractions.

The toxins were further categorised into their respective families, and the relative abundance (%) of each family was calculated by summing up the proportions (T_i) of all the toxins within that family. Mass spectrometry data is deposited to the ProteomeXchange Consortium via the PRIDE partner repository with data identifiers PXD046781 [93]. The list of toxins identified through tandem MS and their relative proportions are tabulated in Additional file 2: Tables S2A-E.

Biochemical and pharmacological analyses

Differences in biochemical and pharmacological activities of venoms of mother snakes and their offspring were evaluated using various assays, including PLA₂, SVSP, SVMP, LAAO and fibrinogenolysis. All concentrations of venoms and reagents used in the assays were based on dose–response experiments performed as previously described.

Colourimetric PLA2 assay

A chromogenic lipid substrate, 4-nitro-3-[octanoyloxy] benzoic acid (NOB; Enzo Life Sciences, New York, NY, USA), was used to determine the PLA₂ activity of venoms. A known amount of venom (5 µg) was incubated at 37° C for 40 min with 500 mM substrate in a total reaction volume of 200 µl made up using the NOB buffer (10 mM Tris–HCl, 10 mM CaCl₂, 100 mM NaCl, pH 7.8). The assay kinetics was recorded by measuring the absorbance at 425 nm every 10 min using an Epoch 2 microplate spectrophotometer (BioTek Instruments, Inc., USA). Further, a standard curve was plotted with 4 M NaOH and varying concentrations of the NOB substrate (4 to 130 nmol) following an identical protocol. The amount of the phospholipid substrate (nmol) cleaved per minute per mg of the venom was calculated from the standard curve [94, 95].

LAAO assay

The LAAO activity of venoms was evaluated using a modified colourimetric assay described previously [95, 96]; 200 µl of the L-leucine substrate solution (600 µl of 5 mM L-leucine, 120 µl of 5 IU/mL horseradish peroxidase, 1200 µl of 2 mM o-phenylenediamine dihydrochloride, 6000 µl of 50 nM Tris–HCl buffer at pH 8) and 0.5 µg of the crude venom were preincubated separately at 37 °C for 10 min. Following this, the venom and substrate solution was mixed and incubated at 37 °C. After 10 min, the reaction was terminated by the addition of 2 M H₂SO₄ solution and the amount of H₂O₂ released was quantified by measuring the absorbance of the solution at 492 nm using an Epoch 2 microplate spectrophotometer (BioTek Instruments, Inc., USA). A standard curve of H₂O₂ with concentrations between 0 and 4.895 nmol/min was constructed using a similar protocol. The LAAO activity of venom samples was estimated as the nmol of H₂O₂ released per minute per mg of venom using this standard curve.

Protease assay

The snake venom protease activity was assessed using an azocasein substrate following a protocol described earlier [97]. To assess the relative contribution of SVMP and SVSP to the overall proteolytic activity, 10 µg of the crude venom was incubated at 37 °C for 15 min with 0.1 M EDTA or PMSF, respectively [98, 99]. Venoms were also incubated in the presence and absence of both inhibitors as control. Following treatment with the inhibitor, samples were incubated with the azocasein substrate (400 µg) at 37 °C for 90 min. After stopping the reaction with trichloroacetic acid (200 µl), the reaction mixture was centrifuged at 1000 × g for 5 min to remove the cleaved products. Equal volumes of 0.5 M NaOH were added to the supernatant, and the absorbance was measured at 440 nm in an Epoch 2 microplate spectrophotometer (BioTek Instruments, Inc., USA). The relative protease activity of venoms was calculated against a purified bovine pancreatic protease (Sigma-Aldrich, USA).

Haemolytic assay

The direct haemolytic effect was assessed by treating known concentrations of venoms (5 to 40 µg) with red blood cells (RBCs) isolated from healthy human volunteers. Briefly, blood was collected in a 3.2% sodium citrate tube and centrifuged at 3000 × g for 5 min at 4 °C. Post the separation of the platelet-poor plasma (PPP), RBCs were washed with PBS to remove residual proteins. Following this, venoms were added to the RBC suspension in a ratio of 1:9 and incubated at 37 °C for 24 h. The amount of heme released during the breakdown of RBCs was then monitored by measuring the absorbance of the supernatant at 540 nm [100]. The relative haemolytic activity was calculated using Triton X as a positive control.

Fibrinogen quantification assay

The ability of venoms to reduce the total active fibrinogen content in plasma was assessed through von Clauss fibrinogen assay using a commercially available Fibroquant kit (Tulip Diagnostics, Goa, India). A standard curve was constructed using various dilutions of a fibrinogen calibrator solution in Owren’s buffer (Undiluted, 1:5 and 1:10 v/v) following the manufacturer’s protocol. The time to form the first fibrin clot was plotted against the total fibrinogen concentration in mg/dl. PPP separated from human blood, collected after the removal of RBCs, was diluted in Owren’s buffer (1:10 v/v), and 40 µg of venom was added. The mixture was incubated at 37 °C for 180 s. In parallel, bovine thrombin was activated by separately incubating at 37 °C for 60 s. The time taken for fibrin clot formation following the addition of thrombin to the plasma–venom mixture was measured using the optomechanical principle in a Hemostar XF 2.0 (Tulip Diagnostics, Goa, India).

Fibrinogenolytic assay

The fibrinogenolytic activity of venom was electrophoretically assessed using a previously described methodology [101]. A predefined quantity of venom (1.5 µg) was dissolved in phosphate-buffered saline (PBS; pH: 7.4) and incubated with 1 µl of SVMP- (0.1 M EDTA) and/or SVSP- (0.04 M PMSF) inhibitors at 37 °C for 15 min [102]. Following this, 15 µg of human fibrinogen (Sigma-Aldrich, USA) was added to the treatment groups and the samples were further incubated at 37 °C for 60 min. After stopping the reaction by the addition of equal volumes of sample loading buffer (1 M Tris–HCl, pH 6.8; 50% Glycerol; 0.5% Bromophenol blue; 10% SDS; and 20% β-mercaptoethanol) and heat treatment at 70 °C for 10 min, samples were subjected to a 15% SDS-PAGE. Densitometric analyses were performed using the ImageJ software to assess the cleavage pattern of the three fibrinogen subunits in the treatment groups compared to the untreated control [88].

Nicotinic acetylcholine receptor specificity

The binding kinetics between cobra venom neurotoxins and nAChRs from various prey/predatory organisms were determined using a Biolayer interferometry (BLI) assay in an Octet Red 96 system (ForteBio, Sartorius, Fremont, CA). The amino acid sequence of this 14 residue-long orthosteric binding site of the ɑ-neuroreceptor was obtained from the Uniprot database (Additional File 2: Table S3) and chemically synthesised by Genpro Biotech (Uttar Pradesh, India). Following previously described protocols [103], the mimotopes were synthesised with two molecules of aminohexanoic acid (Ahx) spacers to separate the amino acid sequence from biotin to avoid steric hindrance. Additionally, disulphide bridges were replaced with Ser-Ser to prevent thiol oxidation without affecting the binding efficiency, as demonstrated previously [104, 105]. Lyophilised stocks of the mimotopes were reconstituted in 100% dimethyl sulfoxide (DMSO) and stored at − 80 °C until use. A final working concentration of 1 µg/mL was prepared by diluting the mimotope stock in the running buffer [1 × Dulbecco's phosphate-buffered saline (DPBS) with 1% BSA and 0.05% Tween-20]. Similarly, the reconstituted venom samples were also diluted to a concentration of 25 µg/mL in the running buffer. The assay was performed with 5 µg of venom and 0.2 µg of the mimotope solution per well in a black 96-well microtiter plate.

The streptavidin sensors were hydrated in the running buffer for 30 min before the beginning of the experiment. The following data acquisition programme was set up: 60 s baseline, 600 s loading, 100 s baseline, 400 s association and 200 s dissociation. The running buffer and mimotopes were introduced, following the loading of venoms onto the sensor tip. A running buffer control was included along with each venom sample, and the thickness of the molecular layer or wavelength shift (in nm) was recorded and plotted against time (sec) to visualise the association and dissociation kinetics. The exported data was inter-step corrected by aligning the baseline to the Y-axis of the initial baseline step, and the dissociation curve was aligned with the association curve. The high-frequency noise in the data was also removed using the Savitsky–Golay filter. The area under the association and dissociation curve was calculated using GraphPad Prism 8.4.3 (GraphPad Software Inc., La Jolla, CA, USA).

In vivo venom toxicity and morbidity

The variation and similarities in the venom lethalities and haemorrhage-inducing potencies were evaluated in the mouse model.

The median lethal dose

Variation in lethal potencies of venoms across developmental stages was determined by estimating the median lethal dose using the mouse model (Mus musculus) of envenoming [106]. The LD₅₀ value indicates the minimum amount of venom required to kill 50% of the test population injected with venom. Five graded concentrations of crude venoms were intravenously injected into the tail vein of five male CD-1 mice (18–22 g) per venom dose group. Subsequently, the mortality in each group was recorded post a 24-h time window. These values were plotted against venom doses to estimate the LD₅₀ using Probit statistics [107].

The lethal potencies of venoms against house geckos (Hemidactylus frenatus) and house crickets (Acheta domesticus) were determined by intraperitoneally injecting five graded concentrations of venoms into three individuals per dose group as described previously [108]. Venom doses ranging from 1 to 20 mg/kg and 2 to 16 mg/kg were injected into geckos and crickets, respectively. The mortalities were recorded over a period of 24 h. Scorpion venoms were tested against crickets as a positive control. Geckos and crickets injected with 0.9% physiological saline and Ringer’s insect saline [109] served as a negative control, respectively.

The minimum haemorrhagic dose (MHD)

The minimum haemorrhagic dose of a venom (MHD) is defined as the amount of venom (in μg) which induces a 10-mm haemorrhagic lesion within 2–3 h of intradermal injection in mice [106]. Five graded venom concentrations were dissolved in normal saline (50 µl) and intradermally injected into each mouse. Five male CD-1 mice (18–22 gm) were injected per concentration to obtain a statistically significant value. After 3 h, the mice were humanely euthanised, and the dorsal skin patch was examined to assess the diameter of the haemorrhagic lesion. Measurements of lesion diameter from individual mice were taken using a vernier calliper to obtain the mean diameter at a given concentration. Subsequently, a dosage at which the lesion diameter is 10 mm was considered the MHD of a given venom by interpolating from a simple linear equation. A group of five mice injected with normal saline served as a negative control.

In vitro and in vivo antivenomics

Enzyme-linked immunosorbent assay

The in vitro venom recognition potential of commercial antivenoms was evaluated using the previously described ELISA protocol [110]. The details of the antivenoms tested are mentioned in Additional File 2: Table S4. Venom samples (100 ng) diluted in a carbonate buffer (pH 9.6) were coated onto 96-well plates and incubated overnight at 4° C. Following a brief washing step (six times) with Tris-buffered saline with 1% Tween 20 (TBST) to remove the unbound proteins, a blocking buffer (5% skimmed milk in TBST) was added to reduce non-specific binding. After incubating the plate for 3 h at room temperature, the plate was rewashed, and serially diluted antivenoms (Premium Serums, VINS Bioproducts, Bharat Serums, and Haffkine Biopharma at 1 mg/mL initial concentration) were added. The plate was subsequently incubated overnight at 4 °C. This step was followed by brief washing and incubation with horseradish peroxidase (HRP)-conjugated rabbit anti-horse secondary antibody (Sigma-Aldrich, USA) diluted in PBS (1:1000) for 2 h. Post incubation, ABTS (2,2’-azino-bis (3-ethylbenzothiazoline-6-sulfonic acid)) substrate solution (Sigma-Aldrich, USA) was added. The absorbance values were measured at a wavelength of 405 nm for 40 min in the Epoch 2 microplate reader and plotted against the antivenom dilutions to determine the antibody titres. Immunoglobulin (IgG) from unimmunised (naive) horses (Biorad) was used for determining the extent of non-specific binding exhibited by equine antibodies [5]. The non-specific binding cut-off was calculated as the mean absorbance of the negative control plus two times the standard deviation. Further, the titre for a particular antivenom was considered as the first dilution where the binding (mean absorbance) was above this cutoff value. The binding efficacy of the antivenom is considered better if the titre (or dilution) is higher. Additionally, the IC₅₀ concentrations of various antivenoms were also calculated. The non-specific binding cut-off calculated above was considered as the baseline value for calculating the IC₅₀ and the 95% CI.

Median effective dose (ED50)

In vivo neutralisation potency of commercial antivenoms was determined by estimating their ED₅₀ values using WHO-recommended preclinical assays [106]. These values indicate the minimum volume of the reconstituted antivenom required to protect 50% of mice injected with a 5 × LD₅₀ challenge dose of venom. Venoms were preincubated with various concentrations of antivenom at 37 °C for 30 min. These mixtures with five graded concentrations of antivenoms were intravenously injected into the tail vein of five male CD-1 mice (18–22 g) per dose group. Subsequently, the mortality and survival rate in each group were recorded post a 24-h time window. These values were plotted against antivenom doses to estimate the ED₅₀ and 95% CI using Probit statistics [107]. ED₅₀ values were transformed to neutralisation potency after considering the venom dosage, based on the following Eq. [74]. Here, n represents the number of LD₅₀ used as the challenge dose.

$$\text{Neutralisation potency }\left(\frac{\text{mg}}{\textrm{mL}}\right)= \frac{\left(n-1\right)\times {\text{LD}}_{50}\textrm{ of venom }(\frac{\text{mg}}{\textrm{mouse}})}{{\textrm{ED}}_{50}\textrm{ of antivenom }(\textrm{mL})}$$

Statistical analysis

Statistical comparisons in various biochemical assays were performed using one-way ANOVA with Tukey’s and Dunnett’s multiple comparison tests. Probit statistics were used for LD₅₀ and ED₅₀ experiments to calculate the median doses and the 95% CI. Simple linear regression curve fitting was performed for the MHD experiment, and the doses were interpolated from the equation. Two-way ANOVA was used to compare the area under the association and dissociation curves in BLI. For calculating the in vitro IC₅₀ values of antivenoms and assessing if the values are significantly different between various antivenoms, the non-linear regression curve fitting and the sum of squares F test were used. All statistical analyses were performed in GraphPad Prism (GraphPad Software 8.0, San Diego, California USA, www.graphpad.com).