Skip to main content
Download PDF

Venom-gland transcriptomic, venomic, and antivenomic profiles of the spine-bellied sea snake (Hydrophis curtus) from the South China Sea

BMC Genomics volume 22, Article number: 520 (2021) Cite this article

Abstract

Background

A comprehensive evaluation of the -omic profiles of venom is important for understanding the potential function and evolution of snake venom. Here, we conducted an integrated multi-omics-analysis to unveil the venom-transcriptomic and venomic profiles in a same group of spine-bellied sea snakes (Hydrophis curtus) from the South China Sea, where the snake is a widespread species and might generate regionally-specific venom potentially harmful to human activities. The capacity of two heterologous antivenoms to immunocapture the H. curtus venom was determined for an in-depth evaluation of their rationality in treatment of H. curtus envenomation. In addition, a phylogenetic analysis by maximum likelihood was used to detect the adaptive molecular evolution of full-length toxin-coding unigenes.

Results

A total of 90,909,384 pairs of clean reads were generated via Illumina sequencing from a pooled cDNA library of six specimens, and yielding 148,121 unigenes through de novo assembly. Sequence similarity searching harvested 63,845 valid annotations, including 63,789 non-toxin-coding and 56 toxin-coding unigenes belonging to 22 protein families. Three protein families, three-finger toxins (3-FTx), phospholipase A2 (PLA2), and cysteine-rich secretory protein, were detected in the venom proteome. 3-FTx (27.15% in the transcriptome/41.94% in the proteome) and PLA2 (59.71%/49.36%) were identified as the most abundant families in the venom-gland transcriptome and venom proteome. In addition, 24 unigenes from 11 protein families were shown to have experienced positive selection in their evolutionary history, whereas four were relatively conserved throughout evolution. Commercial Naja atra antivenom exhibited a stronger capacity than Bungarus multicinctus antivenom to immunocapture H. curtus venom components, especially short neurotoxins, with the capacity of both antivenoms to immunocapture short neurotoxins being weaker than that for PLA2s.

Conclusions

Our study clarified the venom-gland transcriptomic and venomic profiles along with the within-group divergence of a H. curtus population from the South China Sea. Adaptive evolution of most venom components driven by natural selection appeared to occur rapidly during evolutionary history. Notably, the utility of commercial N. atra and B. multicinctus antivenoms against H. curtus toxins was not comprehensive; thus, the development of species-specific antivenom is urgently needed.

Background

Sea snakes, the largest group of current marine reptiles, are widely distributed in many tropical and subtropical waters throughout the Indo-Pacific Ocean [1, 2]. Sea snakes display a diverse array of morphological, physiological, and behavioural traits modified for secondary adaptation to the marine environment; e.g., the venom has evolved to be biochemically simple despite comprising a pharmacologically toxic arsenal [3, 4]. It is generally believed that the adaptive evolution of snake venom composition and function is heavily propelled by natural selection from dietary shift, which might be accompanied by gene birth-and-death processes at the molecular level. A typical case is the shift of the marbled sea snake (Aipysurus eydouxii) to trophic specialization for fish eggs, which involved the dinucleotide deletion of neurotoxin genes and consequently the loss of primary neurotoxic activity in the venom [5]. Defining the venom profiles at the transcriptomic and proteomic levels is helpful for understanding the adaptive evolution and functional complexity of sea snake venom, as well as managing the envenomation caused by sea snakes and, exploring the medically important components and designing the antivenoms against sea snake venom. However, owing to the relatively ambiguous taxonomy of sea snakes, little attention has been paid to snakebites caused by sea snakes and few investigations have been conducted to uncover the profiles of sea snake venom at the ‘-omics’ level. Sea snake venom was previously shown to be mainly composed of three-finger toxins (3-FTx; 26.3–86.9%) and phospholipase A2s (PLA2s; 10.9–66.7%) that vary between and within species [3, 4, 6,7,8,9,10]. Generally, the diversity of the venom proteome is closely related to the profiles of toxin-coding genes, with the sequence conservation of these genes varying during evolutionary history [11]. Therefore, detecting the strength of natural selection on toxin-coding genes together with their homologs would be helpful in explaining the diversity of the venom proteome. Moreover, transcriptional and translational regulation along with post-translational modification is also considered to affect the diverse evolution of venom components. In particular, the abundance and type of toxin transcripts from the venom gland are highly concordant with the venom proteins in several snakes but extremely divergent in others [11,12,13,14,15,16,17,18]. For example, the toxin-coding genes detected in the venom gland of the spine-bellied sea snake (Hydrophis curtus; formerly Lapemis curtus) are not well correlated with the secreted venom proteins [6, 8, 19]. Nevertheless, such a divergence is deduced according to the venom-gland transcriptomic and venom proteomic data on H. curtus specimens from different regions, and still requires in-depth verification using specimens from the same region.

As most sea snakes are usually benign unless they are provoked, the incidence of sea snake bites is much lower than that caused by terrestrial venomous snakes [4, 20, 21]. Actually, over a half of the sea snake bites are related to professional fishing because sea snakes constitute a consistent by-catch of tropical trawl fisheries and the fishermen can get bitted by sea snakes when handling the fishing net [22,23,24]. It is estimated that significant envenomation and death occur in approximately 20 and 3% of sea snake bites, respectively [24], with the death rate reaching 50% if the victims are not properly treated [25]. Generally, the toxin constituents of sea snake venom are predominated by PLA2, 3-FTx and CRISP families, and express various biological/toxic activities, such as myotoxicity (PLA2) [26], postsynaptic neurotoxicity [27], and inhibit smooth muscle contraction and cyclic nucleotide-gated ion channels (CRISP) [28]. Nevertheless, the clinical manifestations of sea snake envenomation are mainly attributed to the 3-FTxs and PLA2s, which can induce systematic symptoms (e.g. paralysis, dysphagia, respiratory failure, and myonecrosis) [24, 25]. Only a single sea snake antivenom is currently commercially available, which is prepared against Hydrophis schistosa (formerly Enhydrina schistosa) venom by Seqirus (a division of Australian company CSL Limited) and is effective in the treatment of envenomation caused by a wide variety of sea snakes [29, 30]. In practice, however, not all victims envenomed by sea snakes receive an injection of sea snake antivenom. In China, for example, sea snakebites are predominantly caused by H. curtus and Hydrophis cyanocinctus, but the patients are either treated with traditional Chinese and western medicine or receive injections of many-banded krait Bungarus multicinctus and Chinese cobra Naja atra antivenoms [31, 32]. Nevertheless, the efficacy and safety of such therapeutic schedules based on antivenoms from phylogenetically related species have yet to be estimated thoroughly.

Hydrophis curtus is the most important of the 16 species of sea snakes in China in terms of the quantity and distribution [33]. This species potentially causes severe illness or death even though it is not defined as a medically important venomous snake [26]. Actually, WHO and the toxicologists are always vigorously appealing people and governments around the world to pay broad attention on snakebites, a neglected tropical disease, including those caused by H. curtus. Therefore, clarifying the global profiles of H. curtus venom would be helpful to understand the clinical symptoms of snakebites caused by H. curtus. Moreover, it is also necessary to qualitatively and quantitatively evaluate the availability of heterologous antivenoms in treatment of snakebites caused by H. curtus. In this study, we therefore applied a strategy combining next-generation sequencing technology and proteomic analysis to unveil the toxin profiles in the venom-gland transcriptome and venom proteome of H. curtus from the South China Sea. The strength of natural selection experienced during the evolutionary history of toxins was tested based on the identified toxin-coding unigenes. Moreover, antivenomic analysis, ELISA and western blotting were used to evaluate and compare the capacity of commercial B. multicinctus and N. atra antivenoms to immunocapture H. curtus venom components.

Results

Sequencing and de novo assembly

A total of 90,909,384 pairs of clean reads were filtered from 93,226,644 pairs of raw reads generated from Illumina sequencing of the venom gland cDNA library of H. curtus, derived from six H. curtus individuals captured as by-catch from the South China Sea, and assembled into 284,495 transcripts (N50/N90 = 1763/290) using Trinity. Subsequently, 148,121 unigenes (N50/N90 = 2229/580) were clustered from these transcripts using Corset, among which 88,806 unigenes passed the quality filter (FPKM > 1). Following similarity alignment between unigenes and sequences in NCBI non-redundant protein (NR)/nucleotide (NT) and Uniprot databases (strict to Serpentes) using Diamond and BLAST, the assemblies were finally categorized into toxin (56), non-toxin (63,789), and unidentified (24,961) groups. The toxin group was expressed at markedly higher redundancies (50,456.48 FPKM/unigene) than those of the non-toxin group (18.02 FPKM/unigene). Moreover, the toxin group comprised only 56 unigenes but accounted for 68% of the total abundance expressed in FPKM of the H. curtus transcriptome, whereas the non-toxin and unidentified groups accounted for 27.7 and 4.3%, respectively (Fig. 1).

Fig. 1
figure 1

Venom-gland transcriptomic profiles of H. curtus. The details are listed in Additional Tables S1 and S2. FPKM, expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced. The percentage of annotated (toxin and non-toxin) and unannotated (unidentified) unigenes in the whole transcriptome was listed in the inserted pie graph. 3-FTx, three finger toxin; PLA2, phospholipase A2; CRISP, cysteine-rich secretory protein; CTL, C-type lectin; SVMP, snake venom metalloproteinase; SVSP, snake venom serine proteinase; PLB, phospholipase B; PDE, phosphodiesterase; HA, hyaluronidase; 5′NT, 5′ nucleotidase; CREGF, cysteine-rich EGF-like domain; VEGF, vascular endothelial growth factor; VF, venom factor; AP, aminopeptidase; NGF, nerve growth factor; AchE, acetylcholinesterase; QC, glutaminyl-peptide cyclotransferases; LAAO, l-amino acid oxidase

Full size image

Venom-gland transcriptomic profile

A total of 22 protein families were derived from 56 toxin-coding unigenes with partial and complete CDS using bioinformatics analyses (Fig. 1 and Additional Table S1). Three toxin families including PLA2 (59.71%), 3-FTx (27.15%), and cysteine-rich secretory protein (CRISP, 12.82%) appeared to dominate the toxin constituents. The remaining 19 toxin families were only expressed at low abundance with a total FPKM of 0.32%, including snake venom metalloproteinase (SVMP), C-type lectin (CTL), snake venom serine proteinase (SVSP), cysteine-type inhibitor (cystatin), phospholipase B (PLB), phosphodiesterase (PDE), hyaluronidase (HA), 5′ nucleotidase (5′NT), PLA2 inhibitor, cysteine-rich with EGF-like domain (CREGF), vascular endothelial growth factor (VEGF), venom factor, aminopeptidase, nerve growth factor (NGF), acetylcholinesterase (AchE), metalloproteinase inhibitor (MP), waprin, glutaminyl-peptide cyclotransferases (QC), and l-amino acid oxidase (Fig. 1). The relative abundance and sequences of toxin unigenes, arranged according to toxin family, are listed in Additional Tables S1 and S2.

Venom proteomic profile

Eleven chromatographic fractions (peaks) were resolved in crude H. curtus venom by RP-HPLC and 32 protein bands from gels were identified by MS analysis (Fig. 2 and Table 1). The proteomic analysis showed that H. curtus venom was dominated by three toxin families including PLA2 (49.36%), 3-FTx (41.94%), and CRISP (8.70%) (Fig. 3 and Table 1). Specifically, chromatographic peaks 1–4 contained 3-FTx accounting for 40.37% of the total venom protein, peaks 5–10 contained PLA2 (49.36%) at high abundance and 3-FTx (1.57%) at very low abundance, and peak 11 contained only CRISP. Additionally, acidic PLA2 (36.08%) and short-neurotoxin (SNX, 33.13%) were the predominant components in PLA2 and 3-FTx families, respectively.

Fig. 2
figure 2

Characterization of the venom proteins of H. curtus from South China Sea. Three milligrams of total venom were applied to a C18 column, and separated as described on Materials and methods. Fractions were collected manually and submitted to molecular determination by SDS-PAGE under reduced conditions (original images were listed in Figure S1). Protein bands were excised, tryptic digested and analyzed by MALDI-TOF/TOF or nESI-MS/MS for their assignment to known protein families. The results are shown in Table 1

Full size image
Table 1 Assignment of the chromatographic fractions and electrophoretic bands from H. curtus venom to protein families
Full size table
Fig. 3
figure 3

Venom proteomic profiles of H. curtus. The details are listed in Table 1. 3-FTx, three finger toxin; PLA2, phospholipase A2; CRISP, cysteine-rich secretory protein; LNX, long chain α-neurotoxin; SNX, short chain α-neurotoxin

Full size image

Moreover, as the proteins identified by MS could be assigned to 11 toxin-coding transcripts of three families (Table 1), the correlation between transcript and protein abundances of individual genes for each toxin family was further evaluated. Non-parametric correlation analysis revealed a weak correlation between both abundances of these 11 toxins with a low Spearman’s rank correlation coefficient (ρ = 0.30), whereas linear regression analysis indicated no correlation (F1,9 = 0.49, P = 0.51) with very low Pearson’s correlation coefficient (R = 0.23). In particular, extreme divergence between transcript and protein abundances could be observed for SNX (Fig. 4 and Table S1). When the analyses were conducted again by excluding SNX, a strong correlation between the abundances of the remaining 10 toxins was detected using non-parametric correlation analysis with ρ = 0.73 in addition to linear regression analysis (F1,8 = 21.8, P = 0.002) with R = 0.86.

Fig. 4
figure 4

Correlation between mRNA and protein abundances of individual gene for each toxin family. SNX was excluded from the analysis. All data were centered log-ratio (clr) transformed. N, number of toxin transcripts; ρ, Spearman’s rank correlation coefficient; R, Pearson’s correlation coefficient

Full size image

Adaptive molecular selection

Considering that not all full-length toxin-coding unigenes could be aligned with multiple homologous sequences with full CDS, we used either codeml or yn00 in PAML 4.8 to evaluate the potential for adaptive natural selection. Furthermore, one unigene (AMP) was excluded from the selection analysis because of the lack of homologs with full CDS. As excessive sequence divergence might reduce the power of the likelihood-ratio test [34], the unigenes from the same toxin family were divided into several groups at a threshold of 10% nucleotide sequence divergence and codeml analysis was conducted separately with their homologs. In total, 23 tests were performed using the codeml program on 15 toxin families including 30 full-length unigenes, with a significance level of 0.002 following Bonferroni correction (Table 2 and Additional Table S3). The selection of three full-length unigenes together with their homologs from the three toxin families was then directly analyzed using yn00.

Table 2 Summary of codeml tests for positive selection of toxins from venom-gland transcriptome in H. curtus
Full size table

To analyze the 3-FTxs, five unigenes were first divided into three groups: 3-FTx(1, 2, 3), 3-FTx(4), and 3-FTx(5). The results indicated that the null hypothesis, the nearly neutral model (M1), could be easily rejected in favour of the positive selection model (M2) in the three tests, with all P < 0.001 following Bonferroni correction (Table 2). Moreover, 28–45% codon sites were found to have experienced positive selection with 6.13 ≤ ω ≤ 9.59. The single-ratio model (M0) indicated that all sites and branches of sequence experienced an average strength of selection with 2.81 ≤ ω ≤ 3.48.

The nucleotide sequence divergence among the four CTL unigenes was > 10%; therefore, these unigenes were analyzed separately with their homologs. In CTL(2) and CTL(3), M1 could be easily rejected in favour of M2, with all P < 0.01 following Bonferroni correction (Table 2). In addition, 16–60% codon sites experienced positive selection with 3.50 ≤ ω ≤ 4.12. However, M1 could not be rejected in CTL(1) and CTL(4) with P > 0.05 in all cases. M0 indicated that all sequence sites and branches experienced an average strength of selection with 0.21 ≤ ω ≤ 2.43.

In PLA2, the null hypothesis, M1, could only be rejected in one group, PLA2(2, 3), in favour of M2 with P < 0.001 following Bonferroni correction (Table 2). In addition, 40% of codon sites experienced positive selection with ω = 6.87. The null hypothesis could not be rejected in PLA2(1) with P = 0.21 following Bonferroni correction (Table 2). M0 showed that all sequence sites and branches experienced an average strength of selection of 0.86 ≤ ω ≤ 2.08.

Unigenes from both SVMP and SVSP families all exhibited > 30% divergence from one another; thus, they were analyzed separately along with their homologs using codeml. For all four tests, M1 could be readily rejected in favour of M2 with all P < 0.001 following Bonferroni correction (Table 2). Approximately 15–18% codon sites experienced positive selection with 5.53 ≤ ω ≤ 5.56 in the SVMP family, with estimated values in the SVSP family of 4–23% codon sites with 7.78 ≤ ω ≤ 9.49. M0 indicated that all sites and branches experienced an average strength of selection with 1.08 ≤ ω ≤ 1.30 in SVMP, whereas that estimated in the SVSP family was 0.73 ≤ ω ≤ 2.25.

The 5′NT, CRISP, and NGF families comprised 2–3 unigenes, with sequence divergence within each toxin family below 10%. Thus, the unigenes from the same family were combined for codeml analysis. The results indicated that the null hypothesis, M1, could be rejected in favour of M2 with all P < 0.05, following Bonferroni correction (Table 2). Additionally, 5% (5′NT), 16% (CRISP), and 10% (NGF) codon sites experienced positive selection with 4.35 ≤ ω ≤ 6.69. M0 showed that all sequence sites and branches experienced an average strength of selection with 0.40 ≤ ω ≤ 1.63.

Of the remaining families subjected to condeml tests, four (Cystatin, PDE, PLB, and VEGF families) easily rejected the null model, M1, and accepted M2 with P < 0.05 in all cases following Bonferroni correction (Table 2). For these four tests, an estimated 7, 12, 2, and 5% codon sites, respectively, experienced positive selection with 3.57 ≤ ω ≤ 7.57. M0 implied that all sites and branches of the four unigenes on average experienced a selection strength of 0.58 ≤ ω ≤ 0.94. Alternatively, tests for HA, PLA2 inhibitor, and QC unigenes could not reject the null hypothesis, M1, with all P > 0.05 following Bonferroni correction (Table 2); M0 indicated ω of 0.54, 0.48, and 0.22, respectively.

As only one homolog with high similarity could be obtained from the database, maximum-likelihood estimation based on dN and dS of the remaining three unigenes was conducted using the yn00 program. Comparing each pair of sequences, we estimated dN of 0.11, 0.04, and 0.003 for AchE (homologous sequence AB852000 from Ovophis okinavensis), CREGF (HQ414087 from Crotalus adamanteus), and MP inhibitor (MG132025 from Bothrops moojeni) of H. curtus, and dS of 0.18, 0.12, and 0.12, with a standard error < 0.04. Thus, we calculated dN/dS = 0.61, 0.31, and 0.02 for these three toxins, respectively, suggesting that these sequences probably experienced purifying selection.

Antivenomic, ELISA, and western blotting evaluation of antivenom efficiency

Using third-generation antivenomic analysis [35] to evaluate the efficacy of commercial monovalent B. multicinctus and N. atra antivenoms in capturing H. curtus venom components, we found that the capacity of N. atra antivenom to immunocapture the H. curtus venom was stronger than that of B. multicinctus antivenom, with increasing abundance of snake venom able to be captured by antivenom as the ratio of antivenom to venom increased (Fig. 5 and Table 3). Based on the ratio of antivenom (20 mg) to venom (50 μg) corresponding to approximately 48-fold molar excess per “10 kDa of toxins”, 53.5% of venom components, comprised mainly of long-neurotoxin (LNX), PLA2, and CRISP, could be immunocaptured by B. multicinctus antivenom, with the non-immunocaptured components including 94.3% of the SNX being minimally recognized. Comparatively speaking, 90.6% of the total venom components could be immunocaptured by N. atra antivenom (Fig. 5A-E and Table 3). However, when the cross-reaction was conducted between 600 μg venom and 20 mg antivenom with a ratio corresponding to approximately 4-fold molar excess per “10 kDa of toxins” (assuming an average molecular mass of 13 kDa for H. curtus venom components), the relative abundance of venom components immunocaptured by the antivenoms decreased dramatically. Specifically, only 17.7% of total venom components could be recognized by B. multicinctus antivenom with the abundance of non-immunocaptured components accounting for 99.1% SNX (peak 1), 84.8% PLA2 and long-neurotoxin (peak 5), 83.7% PLA2 (peak 6), 89.3% PLA2 (peak 8), and 50.4% CRISP (peak 11), respectively (Fig. 5U-W and Table 3). Similarly, only 27.9% of the total venom components could be recognized by N. atra antivenom, with the abundance of non-immunocaptured components accounting for 84.3% SNX (peak 1), 83.6% PLA2 and long-neurotoxin (peak 5), 75.3% PLA2 (peak 6), 76.8% PLA2 (peak 8), and 75.0% CRISP (peak 11), respectively (Fig. 5U, X, Y, and Table 3). In the other three reactions, the antivenoms could immunocapture 43.1% (100 μg), 35.0% (150 μg), and 30.6% (300 μg) of whole venom by B. multicinctus antivenom, and 83.0, 69.6, and 44.0% of that by N. atra antivenom (Fig. 5F-T and Table 3).

Fig. 5
figure 5

Antivenomics analysis of H. curtus venom using commercial antivenom by RP-HPLC. Chromatographic profiles of 50, 100, 150, 300 and 600 μg whole venom components (panels A, F K, P and U) and, of immunocaptured venom components (panels B, D, G, I, L, N, Q, S, V and X) and non-immunocaptured venom components (panels C, E, H, J, M, O, R, T, W and Y) recovered from the affinity columns after incubation with the corresponding amounts of venom. Bm, B. multicinctus; Na, N. atra

Full size image
Table 3 Total and concentration-dependent immunoretained H. curtus venom proteins by heterologous antivenoms
Full size table

We also found that the capacity of the antivenoms to immunocapture the peak fractions differed and that the immunocapture efficiency changed nonsynchronously as the amount of venom increased (Fig. 5 and Table 3). Additionally, the maximal immunocapturing capacity of N. atra antivenom for peaks 1, 6, and 8–11 was higher than that of B. multicinctus antivenom, whereas the effects were opposite for peaks 2–4 and comparable for peaks 5 and 7. It could be calculated that 20 mg of immobilized B. multicinctus and N. atra antivenom presented a maximal immunocapturing capacity of 120.7 and 172.9 μg of H. curtus venom, respectively.

A significant cross-reaction between B. multicinctus/N. atra antivenom and H. curtus venom could be detected by ELISA analysis, with the cross-reactivity showing a progressive increase as the concentration of antivenom increased; moreover, the N. atra antivenom presented relatively higher cross-reactivity than that of the B. multicinctus antivenom (Fig. 6). A similar capacity was demonstrated using western blotting, another conventional protocol that is commonly used to assess the preclinical efficacy of antivenom. The electrophoretic profile revealed that five significant bands with molecular masses ranging from 8 to 21 kDa could be detected in H. curtus venom, of which three (8, 9, and 11 kDa) were present high abundance in venom whereas two bands with molecular masses of 16 and 21 kDa were expressed at relatively low abundance; and another two bands with relatively higher molecular masses of 40 and 60 kDa were expressed at extremely low abundance (Fig. 7A). Mass spectrometry analysis indicated that six protein bands are comprised of more than one protein family except the one with molecular mass of 60 kDa (Fig. 7 and Additional Table S4). According to the peptide intensity, one group of protein bands (8 and 11 kDa) were mainly comprised of relatively high abundance of LNX and SNX, and the other two groups (9, 16 and 60 kDa; 21 and 40 kDa) mainly contained PLA2 and CRISP, respectively (Additional Table S4). Protein bands with molecular masses higher than 21 kDa were strongly recognized by B. multicinctus and N. atra antivenoms but only one (9 kDa) of the three protein bands with the lowest molecular mass could be recognized by these two antivenoms (Fig. 7B and C). The cross-reaction was mainly located in the protein bands with molecular masses of 9, 21, and 40 kDa, with N. atra antivenom affording a higher cross-reaction intensity than that from B. multicinctus antivenom (Fig. 7B and C). Overall, these results indicated that N. atra antivenom presented a higher capacity to immunocapture H. curtus venom than B. multicinctus antivenom and, both antivenoms displayed relatively higher capacity to immunocapture CRISP and PLA2 than that to 3-FTx (LNX and SNX).

Fig. 6
figure 6

Cross-reaction between H. curtus venom and commercial antivenom determined by ELISA. Normal horse serum was used as negative control. Data are expressed as mean ± SD (n = 3)

Full size image
Fig. 7
figure 7

Cross-reaction between H. curtus venom and commercial antivenom determined by western blotting. A SDS-PAGE profiles of venom protein (identified proteins were listed in Table S4); B B. multicinctus antivenom; C N. atra antivenom (original images were listed in Figure S2). PLA2, phospholipase A2; CRISP, cysteine-rich secretory protein; LNX, long chain α-neurotoxin; SNX, short chain α-neurotoxin

Full size image

Discussion

Knowledge of venom-gland transcriptomic and venomic profiles is generally considered to be useful in elucidating the evolutionary process and potential function of snake venom, developing antivenom, and even exploring venom-based drugs [36,37,38,39,40,41,42]. Notably, the maturation and widespread application of high throughput sequencing and MS technologies have rendered it much easier to unravel the complexity of snake venoms worldwide at the -omics level. In the present study on H. curtus from the South China Sea, 22 and 3 toxin families were identified in the venom-gland transcriptome and venom proteome, respectively, using next-generation sequencing technology and a combined proteomic strategy.

Notably, the composition and abundance of H. curtus toxins identified in the present study apparently differed from those in previous investigations either at the mRNA or protein level [6, 8, 19]. For example, the toxin transcripts of the H. curtus venom gland was shown to be mainly comprised of 3-FTx, PLA2, CRISP, and PDGF families using Sanger sequencing, with the first three families accounting for 55.5% of the total venom-gland transcriptome and the 3-FTx family expressed more abundantly than PLA2 at a ratio of 4.4:1 [19]. Comparatively speaking, although 3-FTx, PLA2, and CRISP transcripts also constituted the most abundant components in the venom-gland transcriptome of H. curtus in the present study, accounting for 99.68% of total toxin transcripts, 3-FTx was less abundant than PLA2 at a ratio of 1:2.2. Moreover, 19 toxin families expressed with low abundance (0.32%) among total toxin transcripts were also detected in our sample. Next-generation sequencing technology therefore appeared to be more powerful than Sanger sequencing for detecting the snake venom-gland transcripts with low abundance. Moreover, the venom protein/mRNA sequence database was relatively insufficient 10 years ago, such that a high proportion (24.3%) of transcripts in a previous study could not be matched to the targeted toxins, potentially resulting in additional toxin families being omitted.

When compared with previously investigated H. curtus venoms from Australia and Malaysia [6, 8], the venom from the South China Sea population evaluated in the presented study diverged significantly with regard to the relative abundance of predominant toxins at the protein level. The South China Sea population exhibited the lowest abundance of PLA2 but the highest abundance of 3-FTx, especially SNX, which was approximately 4- and 1.45-fold higher than those from Australia and Malaysia, respectively. Notably, given that SNX was determined to be the most toxic component in H. curtus venom [6], the South China Sea population might possess the highest venomic lethality among the three localities.

Moreover, although the venom-gland transcriptome and venom proteome of H. curtus were dominated by 3-FTx, PLA2, and CRISP families, of which the total abundance expressed at the protein level was highly consistent with that at the mRNA level, an apparent discrepancy was still observed with regard to the composition and abundance of toxins at the family level. Given that comparisons at the family level may not necessarily reflect the mechanisms affecting specific toxin transcripts, a detailed analysis of the correlation between mRNA and protein abundances of 11 toxins for 3-FTx, PLA2, and CRISP families was conducted, with the results agreeing with those at the family level. Such discordance between both levels has been observed in many venomous snakes, being potentially attributed to transcriptional, post-transcriptional, translational, or post-translational regulatory mechanisms [11,12,13, 43,44,45,46] in addition to the time span between collection of the venom and venom-gland tissue [14]. In contrast, some studies revealed little differential expression of toxins at both -omics levels despite individual variation [18, 47, 48]. Notably, a strong correlation between both mRNA and protein levels was detected in the present study in the remaining toxin transcripts when SNX was excluded. Thus, the mechanisms regulating the mRNA and protein abundances of each toxin from H. curtus may differ. Nevertheless, as the majority of H. curtus toxins could only be detected at the mRNA but not the protein level, it remains to be explored whether this phenomenon is regulated by more complex mechanisms such as post-transcriptional silencing and protein buffering.

Considered as a highly effective and biochemically complex weapon involved in predation, digestion, and defense, snake venom is generally considered to have undergone positive selection during adaptive evolution, which has been successfully estimated at individual toxin and venom-gland transcriptomic levels [49,50,51,52,53,54,55,56]. To deal with diverse preys, a relatively high proportion of toxins appears to experience positive selection at transcriptomic level in many terrestrial venomous snakes, such as 24 of 27 full-length toxin-coding transcripts from venom-gland transcriptome were detected to be underwent positive selection in Crotalus adamanteus [52] and, nearly all of the venom-gland toxin ESTs had sites that clearly showed evidence of positive selection signals in B. multicinctus and N. atra [56]. Consistent with the results of these studies, the majority of the toxin transcripts (24/33) from H. curtus, accounting for 73% of the total analyzed transcripts, were deemed in the present study to have experienced positive selection according to the likelihood ratio from M1 and M2; and it was also verified by M7 and M8 (the results were not shown in this paper as they nearly identical to that from M1 and M2). That is to say the adaptive evolution of toxin genes from H. curtus can still be extensively driven by natural positive selection, even as the diet of H. curtus are simpler than those of its terrestrial relatives. Some studies proposed that toxin genes experience both stronger positive selection and lower selective constraint than non-toxin genes, and toxins possessed with the predominant physiological functions and expressed relatively high abundance evolve more rapidly due to the positive selection [53, 57, 58]. Here, it could be easily observed that the average dN/dS of 3-FTx, PLA2 and CRISP families are relatively larger than that of the other toxin families in H. curtus. It is thus confirmed that the force strength of positive selection is similarly stronger in those H. curtus venom toxins expressed important physiological functions and higher abundance, which might imply a potential interplay between abundance and adaptive evolution of toxins. Furthermore, the rapidly adaptive evolution driven by positive selection in toxins might be accompanied with mutations in vital nucleotide sites and the followed reduction (even loss) of toxin efficiency, while purifying selection would probably aid in preserving the venom potency by filtering out these mutations with deleterious effects [59]. It is thus believed that the different strength of purifying selection in H. curtus toxins will provide evolutionary implications in optimizing the trade-off between venom potency and adaption to a simplified diet.

The antivenomic analysis approach based on venomic profile, originally developed as an in vitro platform to qualitatively and quantitatively assess the efficacy of antivenom, has since been updated to the 3rd generation and displays a powerful potential as a substitute for many traditional strategies to evaluate antivenom efficacy [35]. Similar to several recently reported studies using 3rd generation antivenomics [60,61,62], the immunocapture capacity of antivenoms in the present study increased as the incubation amount of H. curtus venom increased but varied according to different venom components. Overall, N. atra antivenom exhibited higher capacity to immunocapture whole H. curtus venom than B. multicinctus antivenom, especially for SNX, PLA2, and CRISP. However, the former demonstrated a relatively lower capability to immunocapture the LNX. This suggests that N. atra antivenom possesses a greater abundance of F(ab′)2 to efficiently recognize H. curtus venom than B. multicinctus antivenom at the same dose, but the amounts of F(ab′)2 to recognize specific venom components are likely different between these two antivenoms. This divergence was also roughly verified by western blotting and ELISA analyses. Given that six of seven protein bands from H. curtus venom directly separated by SDS-PAGE were comprised of more than one protein family, western blotting analysis displayed weak ability to qualitatively and quantitatively evaluate the capacity of both antivenoms to recognize each venom component. Although it was found both antivenoms cannot efficiently recognized most of 3-FTx in H. curtus venom by 3rd generation antivenomic and western blotting analyses, it was still hard to associate the results from both two analyses precisely. Briefly, western blotting analysis only roughly indicated that the capacity of both antivenoms to recognize the venom components with low molecular mass (≤ 15 kDa, mainly 3-FTx and PLA2) was relatively weaker than that for recognizing high molecular mass components (mainly CRISP and PLA2), with no ability to recognize the components (mainly 3-FTx) with molecular masses of 8 and 11 kDa. This might be due to the lack of homologous components between B. multicinctus/N. atra and H. curtus venoms, or because 3-FTx generally induce limited or no antibodies because of their weak immune response.

Given the high similarity in venom composition among closely related snakes, appropriate heterologous antivenoms to treat envenomations caused by snakes with no commercial antivenom might be selected according to the phylogenetic relationship. For example, patients envenomed by sea snakes including H. curtus have been often been treated by B. multicinctus and N. atra antivenoms in some hospitals in China [31]. However, based on the antivenomic profiles generated in the present study, N. atra antivenom appears to be more suitable for treatment of patients envenomed by H. curtus than B. multicinctus antivenom, although the phylogenetic relationship between H. curtus and B. multicinctus is much closer than that between H. curtus and N. atra [63]. Thus, the phylogenetic relationship should not be a unique criterion for selection of alternative antivenoms in clinic treatment of snakebites.

Specifically, based on the maximal immunocapture capacity, one vial of commercial B. multicinctus (approximately 57 mg/mL, 10 mL/vial) and N. atra (124 mg/mL, 10 mL/vial) antivenoms theoretically could be inferred to absorb 3.4 and 10.7 mg H. curtus venom components, respectively, which appeared to approximate 28 and 89% of the average venom yield of 12.0 mg (lyophilized venom, calculated from six specimens in the present study). However, it should be noted that the maximal amount of venom components immunocaptured by the antivenoms may not be equal to that of whole venom; e.g., 20 mg N. atra antivenom could maximally absorb 172.9 μg venom components but only absorbed 69.6% of 150 μg whole H. curtus venom. Moreover, both antivenoms only exhibited high immunocapture efficiency when incubated with relatively low amounts of venom, with up to 53.5 and 90.6% of the total amount of 50 μg venom being absorbed by B. multicinctus and N. atra antivenoms, respectively. The abundance of unabsorbed venom components increased as the ratio of antivenom to venom decreased, especially for SNX, the most toxic component in H. curtus venom [6]. Consistent with this, the clinical efficacy of both B. multicinctus and N. atra antivenoms in the treatment of sea snake envenomations has also been relatively poor [31]. Accordingly, the development and application of species-specific antivenom against H. curtus venom in China is urgently needed.

Conclusion

An integrated omics-strategy was used to reveal the venom-gland transcriptomic, venomic, and antivenomic profiles of H. curtus from the South China Sea. Our findings verify discordance between the venom-gland transcriptome and venom proteome in H. curtus. Briefly, the venom evolved to as a ‘biochemically simple’ weapon at the protein level to adapt to the simplified diet available when the sea snake returned to the marine environment; conversely, the toxin transcripts with high diversity in the venom-gland transcriptome of H. curtus could be defined as a ‘genetically complicated’ mixture at the mRNA level. In addition, the majority of the full-length toxin-coding unigenes may have experienced positive selection in the evolutionary history. However, it should not be ignored that the transcripts of our sample have not been assembled based on a targeted reference genome, which might potentially hamper the identification of toxins that are not detected in venom-gland transcriptome and proteome and, exploration of more insights into toxin evolution including gene duplication and pseudogenization, etc. Moreover, N. atra antivenom has a stronger capability to immunocapture H. curtus venom components than does B. multicinctus antivenom, but both antivenoms only exhibit high immunocapture efficiency when they are incubated with relatively low amounts of venom. Thus, the clinical application of commercial N. atra and B. multicinctus antivenoms is not recommended in treatment of patients envenomed by H. curtus. Rather, the development and application of species-specific antivenom of H. curtus venom is urgently required in China. Taken together, our results provide a foundation for additional studies of population-specific venom molecular and evolutionary analyses and the development of effective H. curtus antivenom, and the antivenomic strategy based on venomic profile for qualitatively and quantitatively evaluating the efficiency of antivenom would present a well development future.

Methods

Venom and antivenoms

We obtained six adult specimens of H. curtus as by-catch in the South China Sea and transported them to our laboratory at Hainan Tropical Ocean University, and maintained them in a circulating sea water system for venom extraction. Fresh venom from individual snakes was extracted repeatedly at a 15-day interval using a 100-μl plastic pipette, centrifuged to remove impurities for 15 min at 10,000 g 4 °C, lyophilized, pooled equally, and stored at − 80 °C until use. Commercial monovalent B. multicinctus antivenom (batch 20,070,601, expiry date: 06/2010) and N. atra antivenom (batch 20,070,601, expiry date: 06/2010) consisting of purified F(ab′)2 fragments from the plasma of hyperimmunized horses were purchased from Shanghai Serum Biological Technology Co., Ltd., China. Both antivenoms were subdivided, lyophilized, and stored at − 80 °C immediately upon receipt in the laboratory. Snakes were collected under the permit (2017-09-12) issued by Hainan Tropical Ocean University, and our experimental procedures was approved by the Animal Research Ethics Committee of Hangzhou Normal University (AREC2019109).

Venom gland cDNA synthesis and sequencing

The snakes were sacrificed for venom-gland transcriptomic analysis when sufficient crude venom was accumulated. After the last round of venom extraction, the snakes were allowed to recover for 4 days to maximize the venom gland transcription and then anesthetized heavily with sodium pentobarbital (s.c. 30 mg/kg) until they showed no tail-retraction reflex. Two venom glands of each snake were removed and split into pieces with a diameter of < 2 mm, then pooled and permeated with RNAlater® (Qiagen) at 4 °C overnight before transferring to − 80 °C storage until use. The snakes then were euthanized by injection of sodium pentobarbital (i.p. 100 mg/kg), and the carcasses were preserved with 10% formalin and deposited in the animal collections of Hainan Tropical Ocean University. Total RNA of the pooled venom gland tissue from each sample was isolated using TRIzol (Life Technologies) according to the manufacturer’s instructions. RNA was purified, concentrated, and then resuspended in 100 μl THE Ambion® RNA storage solution (Life Technologies). The degree of contamination and degradation of RNA in each sample was monitored by 1% agarose gel electrophoresis. The purity and integrity of RNA was assessed using a Implen NanoPhotometer (Implen GmbH) and Agilent 2100 Bioanalyzer system (Agilent Technologies), respectively, whereas the concentration of RNA was measured using a Qubit 2.0 fluorometer (ThermoFisher Scientific). Subsequently, an RNA mixture for sequencing was prepared by mixing equal amounts of RNA from the above six samples.

The cDNA library for RNA sequencing was generated using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (New England Biolabs) according to the manufacturer’s instructions. In brief, mRNA was purified and enriched from total RNA with oligo (dT)-attached magnetic beads. Then, the mRNA was fragmented and used as a template to synthesize first-strand cDNA with reverse transcriptase (RNase H) and random hexamer primers. The following step for synthesizing second-strand cDNA was conducted with dNTPs, RNase H and DNA polymerase I based on first-strand cDNA. The resultant double-stranded cDNA was purified using AMPure XP beads (Beckman Coulter) and subjected to the processes of end repair and ligation of a poly (A) tail and adapters. The adapter-ligated fragments were preferentially screened for those at 250–300 bp in length, amplified by PCR, and purified using AMPure XP beads to create the final cDNA library. Quality of the library was assessed using the Agilent 2100 Bioanalyzer system (Agilent Technologies). High throughout sequencing of the cDNA library was performed using the Illumina HiSeq™2500 platform at Novogene Bioinformatics Technology Co., Ltd., China (www.novogene.cn).

Transcriptome assembly, annotation, and quantification

Prior to transcriptome assembly, raw reads (raw data) from Illumina sequencing were cleaned by eliminating the reads containing adapters and > 10% unknown sequences (‘N’s), as well as the reads with low quality (Q ≤ 20). The clean reads (remaining reads) were assembled into contigs using Trinity 2.4.0 based on Inchworm, Chrysalis, and Butterfly modules [64]. Initially, clean reads were grouped into different read sets and assembled into a k-mer dictionary (k = 25), which was then developed into a collection of linear contigs by greedily searching using Inchworm. If the contigs shared at least one k – 1-mer and the reads spanned the junction between contigs, they were pooled and built into de Bruijn graphs using Chrysalis. Finally, the de Bruijn graphs with clean reads and paired-end reads were reconciled and the full-length transcripts for spliced isoforms and paralogous sequences were reconstructed and arranged using Butterfly. Subsequently, the longest sequence with no redundancy was generated from the transcripts in each gene locus by Corset 1.05 and defined as a unigene, which was used as a reference for downstream analyses. Gene annotation was performed by BLAST 2.2.28 searching against the NCBI NT database and Diamond 0.8.22 searching against the NCBI NR and UniProt protein databases (strict to the taxa Serpentes). The E-value threshold was set to 1 × 10− 5. To estimate the abundance of unigenes, all clean reads were initially matched with the reference unigenes using RSEM 1.2.15 [65]. The number of reads matched to a given unigene was defined as the readcount, which was then calculated as FPKM [66] for evaluating the abundance.

Isolation and characterization of venom proteins

Venom powder (3 mg) collected from the six snakes was reconstituted in 0.1% TFA and centrifuged for 15 min at 10,000 g, 4 °C. The supernatant was collected and applied to a Kromasil 300 RP-C18 column (250 × 4.6 mm, 5 μm; AkzoNobel) using a Waters E600 HPLC system (Waters). The venom proteins were separated at a flow rate of 1 ml/min with a mobile phase system of 0.1% TFA in water (solution A) and 100% ACN (solution B) at the following linear gradient: 0–15% B for 30 min, followed by 15–45% B for 120 min and 45–70% B for 20 min. Protein detection was conducted at 215 nm. Fractions were collected and concentrated in a CentriVap® Centrifugal Concentrator (Labconco). Protein concentrations of venom fractions were determined according to the Bradford method [67]. The proteins in each fraction were separated by 18% SDS-PAGE. The gels were stained in 0.2% CBB R-250 and scanned using a UMAx2100 densitometer (Novax Technologies).

Protein bands in the gels were excised, destained with 0.1 M NH4HCO3 in 30% ACN, rinsed with 0.1 M NH4HCO3, and then digested with trypsin (Promega) for 20 h at 37 °C. The peptide mixture was lyophilized, and re-dissolved in 2 μl of 20% ACN. The solution was spotted onto a sample holder and air-dried, to which 0.5 μl of 5 mg/ml α-cyano-4-hydroxycinnamic acid (ABI) in 50% ACN and 0.1% TFA was added; then, the sample was dried completely and subjected to 4800 Plus MALDI-TOF/TOF-MS mass spectrometer (AB SCIEX) according to the instruction manual. For LC-MS/MS analysis, the peptide mixture was re-dissolved in 0.1% TFA and desalted using a Zorbax 300 SB-C18 column (150 × 0.3 mm, 5 μm; Agilent Technologies), then separated by reverse phase capillary HPLC using an RP-C18 column (150 × 0.15 mm, 5 μm; Column Technology Inc.) with a mobile phase system of 0.1% FA in water (solution A) and 84% ACN in 0.1% FA (solution B) as follows: 4–50% B for 30 min, followed by 50–100% B for 4 min and 100% B for 1 min. The eluted peptides were applied to a Q Exactive mass spectrometer (ThermoFisher Scientific) according to the manufacturer’s instructions. The raw MS/MS spectra were interpreted using FlexAnalysis 3.4 or Xcalibur 2.2 software and the assignment of amino acid sequence similarity was executed using PEAKS X against the UniProt Serpentes database or an in-house database constructed using toxin transcripts extracted from the H. curtus venom-gland transcriptome. The mass tolerance was set at 0.4 Da for MALDI-TOF/TOF-MS and 0.1 Da for LC-MS/MS. Carbamidomethyl (C) was set as a fixed modification and oxidation (M) was set as a variable modification.

Integration of the collected HPLC fractions and densitometry of the protein bands in the SDS-PAGE electropherogram were used to estimate the relative abundance of venom composition according to Calvete et al. [68] and Shan et al. [69]. Briefly, the relative abundance of fractions was calculated by peak area measurement using Empower 2.0 software. When the fractions presented only one protein band in SDS-PAGE, the relative abundance was directly obtained from the peak area measurement, whereas for the fractions presenting more than one protein band, the relative abundance of each band was estimated by densitometry using Tan4100 software.

Detection of adaptive molecular evolution

Tests of adaptive molecular evolution were performed on the transcripts with full-length CDS from 18 toxin families using positive selection analysis. Prior to analysis, sequences homologous to transcripts were obtained from the NCBI NT database with an approximate threshold of 10% nucleotide sequence divergence. Sequences were aligned using Geneious 4.8.3 on the basis of the amino-acid sequence. Gaps, stop codons, and signal peptides were excluded from all analyses. The best-fitting model for partitions was determined using PartitionFinder 1.1.1 [70] based on the Bayesian Information Criterion and greedy search. Then, a Bayesian phylogeny was constructed using BEAST 2.2 [71] and each analysis was conducted by performing four replicate searches with 100 million generations, retaining trees every 10,000 generations. The chains were ensured to converge and mix adequately using Tracer 2.2. The maximum clade credibility tree for the target tree was obtained using the TreeAnnotator 2.2 suite, with the first 10% of sampled generations being excluded.

A likelihood-ratio test for positive selection was performed using codeml in PAML 4.8, based on five models [72] as follows. The nearly neutral model (M1), which is defined as the null model, allows for a group of sites to have experienced neutral selection (dN/dS = 1) and another group to be constrained under purifying selection (dN/dS < 1) in evolutionary history. The positive model (M2), which is defined as the alternative model, advances an additional parameter that can indicate the group of sites that experienced positive selection (dN/dS > 1). To test whether a group of sites experienced positive selection, the difference in log likelihoods between these two models was compared to a χ2 distribution with two degrees of freedom. A similar test by comparing models M7 (beta) and M8 (beta with positive selection) was conducted to verify the initial results, using a χ2 distribution with two degrees of freedom. An overall dN/dS was indicated by the M0 model, which introduces a single ratio for all sites based on an average dN/dS across the entire phylogeny. If the transcript could only be aligned with one full-length homolog from the database, then dN and dS were directly analyzed by yn00 in PAML 4.8, and the dN/dS value for selection assessment was calculated manually.

Antivenomics analysis

The capability of commercial monovalent B. multicinctus antivenom and N. atra antivenom to recognize H. curtus venom was assessed using the third-generation antivenomics analysis [35]. The antivenom was dissolved and dialyzed against MilliQ water, then lyophilized and re-dissolved in coupling buffer (0.2 M NaHCO3, 0.5 M NaCl, pH 8.3). An immunoaffinity column coupled with B. multicinctus or N. atra antivenom was prepared simultaneously. CNBr-activated Sepharose 4B matrix (1 ml) (GE Healthcare) was packed in a 3 ml column and washed with 15 CV (matrix volumes) of 1 mM ice-cold HCl, followed by 3 CV of coupling buffer. The matrix was incubated overnight at 4 °C in 0.5 CV of coupling buffer containing 70 mg antivenom. To calculate the antivenom coupling yield, the concentrations of antivenom before and after incubation with the matrix were determined spectrophotometrically at 280 nm based on an extinction coefficient (ε) of 1.36 for 1 mg/ml protein in a 1 cm light path length. The matrices of the two columns were coupled with 20.1 mg B. multicinctus and 20.5 mg N. atra F(ab′)2 antivenom fragments. After the coupling, the non-reacting groups were blocked gently with 1 CV of blocking buffer (0.1 M Tris-HCl, pH 8.0) for 4 h at room temperature on an orbital shaker. Both columns containing matrix were alternately washed with six repetitions of 3 CV of low (0.1 M acetate buffer, 0.5 M NaCl, pH 4.0) and high (0.1 M Tris-HCl, pH 8.5) pH buffer and equilibrated with 3 CV of binding buffer (20 mM PBS, pH 7.2). Finally, the matrix coupled with 20 mg of antivenom was retained in each column. For immunoaffinity analysis, 0.5 ml binding buffer containing 50, 100, 150, 300, and 600 μg H. curtus venom was loaded onto the column and incubated gently for 1 h at room temperature on an orbital shaker. After recovering the non-retained venom fractions using 3 CV binding buffer, the retained fractions were eluted and collected with 3 CV of 0.1 M glycine-HCl, pH 2.0, and immediately brought to neutral pH with 1 M Tris-HCl (pH 9.0). Both non-retained and retained fractions were concentrated using a CentriVap® Centrifugal Concentrator and analyzed by RP-HPLC as described above.

Enzyme-linked immunosorbent assay (ELISA)

Microplates (96 wells) were coated with 100 μl H. curtus venom proteins (2 μg/ml in 0.1 M Na2CO3-NaHCO3, pH 9.6) per well overnight at 4 °C. The unbound proteins were washed off with PBST (0.05% Tween-20 in 10 mM PBS, pH 7.4) and the plate was blocked with 150 μl 2% BSA in PBST at 37 °C for 1 h. After washing three times, 100 μl of serially diluted horse serum/commercial B. multicinctus/N. atra antivenom (initial concentration 4 μg/μl) in PBST containing 1% BSA was added into each well and incubated at 37 °C for 1 h. Subsequently, the plate was washed again and 1:10000 diluted HRP-labelled anti-horse IgG (Sigma) was added into each well and incubated at 37 °C for 1 h. Finally, the unbound secondary antibodies were thoroughly rinsed from the plate with PBST. Aliquots (100 μl) of substrate solution (0.5 mg/ml OPD and 0.006% H2O2 in 0.15 M citrate buffer, pH 5.0) was added to each well and incubated at room temperature for 20 min. The reaction was stopped with 50 μl 2.5 M sulphuric acid and the absorbance was recorded at 490 nm using a SpectraMax 384 microplate reader (Molecular Devices).

Western blotting

Venom proteins of H. curtus were separated by 18% SDS-PAGE under reduced conditions according to Laemmli (1970). After electrophoresis, the proteins on the gel were either transferred to a PVDF membrane (0.45 μm, Millipore) in a semi-dry system (Bio-Rad) or stained with CBB R-250. The membrane was then blocked with 5% BSA in TBST (20 mM Tris-HCl pH 8.0, containing 150 mM NaCl and 0.05% Tween-20) overnight at 4 °C. After washing with TBST, the membrane was incubated with 1:500 diluted commercial antivenom at 37 °C for 1 h. Then, the membrane was washed and incubated with 1:3000 diluted HRP-labelled antihorse IgG at 37 °C for 1 h. After washing off the unbound secondary antibodies, the membrane was incubated with Pierce™ ECL western blotting substrate and exposed to X-ray film. The film was scanned using a UMax2100 densitometer (Novax Technologies) and analyzed with Tan4100 software. In addition, the protein bands in the CBB R-250 stained gel were excised and identified by LC-MS/MS according to the procedure above mentioned.

Statistical analyses

To detect the correlation between mRNA and protein level abundance of individual genes for each toxin family, the data were centred log-ratio (clr) transformed prior to analyses according to Rokyta et al. [73]. Two coefficients (Spearman’s rank correlation coefficient and Pearson’s correlation coefficient) used for assessment of correlation were calculated by non-parametric correlation and linear regression analyses using Statistica 8.0. The significance level was set at α = 0.05.

Availability of data and materials

The venom-gland transcriptome raw data were submitted to the NCBI Sequence Read Archive (SRA) under the accession number PRJNA670846.

Abbreviations

3-FTx:

Three finger toxin

5′NT:

5′ nucleotidase

AchE:

Acetylcholinesterase

AP:

Aminopeptidase

Bm:

B. multicinctus

CREGF:

Cysteine-rich EGF-like domain

CRISP:

Cysteine-rich secretory protein

CTL:

C-type lectin

HA:

Hyaluronidase

dN :

Nonsynoymous substitution rate

dS :

Synonymous substitution rate

LAAO:

l-amino acid oxidase

LNX:

Long chain α-neurotoxin

Na:

N. atra

NGF:

Nerve growth factor

PAML:

Phylogenetic analysis by maximum likelihood

PDE:

Phosphodiesterase

PLA2 :

Phospholipase A2

PLB:

Phospholipase B

QC:

Glutaminyl-peptide cyclotransferases

SNX:

Short chain α-neurotoxin

SVMP:

Snake venom metalloproteinase

SVSP:

Snake venom serine proteinase

VEGF:

Vascular endothelial growth factor

VF:

Venom factor

References

  1. Karthikeyan R, Balasubramamian T. Species diversity of sea snake (Hydrophiidae) distributed in the coramantal coast (East Coast of India). Int J Zool Res. 2007;3(3):107–31.

    Article  Google Scholar 

  2. Rasmussen AR, Murphy JC, Ompi M, Gibbons JW, Uetz P. Marine Reptiles. PLoS One. 2011;6(11):e27373. https://doi.org/10.1371/journal.pone.0027373.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Tan CH, Tan KY, Lim SE, Tan NH. Venomics of the beaked sea snake, Hydrophis schistosus: a minimalist toxin arsenal and its cross-neutralization by heterologous antivenoms. J Proteomics. 2015;126:121–30. https://doi.org/10.1016/j.jprot.2015.05.035.

    Article  CAS  PubMed  Google Scholar 

  4. Calvete JJ, Ghezellou P, Paiva O, Matainaho T, Ghassempour A, Goudarzi H, et al. Snake venomics of two poorly known Hydrophiinae: comparative proteomics of the venoms of terrestrial Toxicocalamus longissimus and marine Hydrophis cyanocinctus. J Proteomics. 2012;75(13):4091–101. https://doi.org/10.1016/j.jprot.2012.05.026.

    Article  CAS  PubMed  Google Scholar 

  5. Li M, Fry BG, Kini RM. Eggs-only diet: its implications for the toxin profile changes and ecology of the marbled sea snake (Aipysurus eydouxii). J Mol Evol. 2005;60(1):81–9. https://doi.org/10.1007/s00239-004-0138-0.

    Article  CAS  PubMed  Google Scholar 

  6. Tan CH, Tan KY, Ng TS, Sim SM, Tan NH. Venom proteome of spine-bellied sea snake (Hydrophis curtus) from Penang, Malaysia: toxicity correlation, immunoprofiling and cross-neutralization by sea snake antivenom. Toxins. 2019;11(1):3.

    Article  CAS  Google Scholar 

  7. Durban J, Sasa M, Calvete JJ. Venom gland transcriptomics and microRNA profiling of juvenile and adult yellow-bellied sea snake, Hydrophis platurus, from Playa del Coco (Guanacaste, Costa Rica). Toxicon. 2018;153:96–105. https://doi.org/10.1016/j.toxicon.2018.08.016.

    Article  CAS  PubMed  Google Scholar 

  8. Neale V, Sotillo J, Seymour J, Wilson D. The venom of the spine-bellied sea snake (Hydrophis curtus): proteome, toxin diversity and intraspecific variation. Int J Mol Sci. 2017;18(12):2695. https://doi.org/10.3390/ijms18122695.

    Article  CAS  PubMed Central  Google Scholar 

  9. Tan KY, Tan CH, Fung SY, Tan NH. Neutralization of the principal toxins from the venoms of Thai Naja kaouthia and Malaysian Hydrophis schistosus: insights into toxin-specific neutralization by two different antivenoms. Toxins. 2016;8(4):86. https://doi.org/10.3390/toxins8040086.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Lomonte B, Pla D, Sasa M, Tsai W-C, Solórzano A, Ureña-Díaz JM, et al. Two color morphs of the pelagic yellow-bellied sea snake, Pelamis platura, from different locations of Costa Rica: snake venomics, toxicity, and neutralization by antivenom. J Proteomics. 2014;103:137–52. https://doi.org/10.1016/j.jprot.2014.03.034.

    Article  CAS  PubMed  Google Scholar 

  11. Casewell NR, Wagstaff SC, Wüster W, Cook DAN, Bolton FMS, King SI, et al. Medically important differences in snake venom composition are dictated by distinct postgenomic mechanisms. Proc Natl Acad Sci U S A. 2014;111(25):9205–10. https://doi.org/10.1073/pnas.1405484111.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Rodrigues RS, Boldrini-França J, Fonseca FPP, de la Torre P, Henrique-Silva F, Sanz L, et al. Combined snake venomics and venom gland transcriptomic analysis of Bothropoides pauloensis. J Proteomics. 2012;75(9):2707–20. https://doi.org/10.1016/j.jprot.2012.03.028.

    Article  CAS  PubMed  Google Scholar 

  13. Margres MJ, McGivern JJ, Wray KP, Seavy M, Calvin K, Rokyta DR. Linking the transcriptome and proteome to characterize the venom of the eastern diamondback rattlesnake (Crotalus adamanteus). J Proteomics. 2014;96:145–58. https://doi.org/10.1016/j.jprot.2013.11.001.

    Article  CAS  PubMed  Google Scholar 

  14. Tan CH, Tan KY, Fung SY, Tan NH. Venom-gland transcriptome and venom proteome of the Malaysian king cobra (Ophiophagus hannah). BMC Genomics. 2015;16(1):1–21.

    Article  Google Scholar 

  15. Modahl CM, Frietze S, Mackessy SP. Transcriptome-facilitated proteomic characterization of rear-fanged snake venoms reveal abundant metalloproteinases with enhanced activity. J Proteomics. 2018;187:223–34. https://doi.org/10.1016/j.jprot.2018.08.004.

    Article  CAS  PubMed  Google Scholar 

  16. Xu N, Zhao H-Y, Yin Y, Shen S-S, Shan L-L, Chen C-X, et al. Combined venomics, antivenomics and venom gland transcriptome analysis of the monocoled cobra (Naja kaouthia) from China. J Proteomics. 2017;159:19–31. https://doi.org/10.1016/j.jprot.2017.02.018.

    Article  CAS  PubMed  Google Scholar 

  17. Aird SD, da Silva NJ Jr, Qiu L, Villar-Briones A, Saddi VA, de Campos Telles PM, et al. Coralsnake venomics: analyses of venom gland transcriptomes and proteomes of six Brazilian taxa. Toxins. 2017;9(6):187. https://doi.org/10.3390/toxins9060187.

    Article  CAS  PubMed Central  Google Scholar 

  18. Hofmann EP, Rautsaw RM, Strickland JL, Holding ML, Hogan MP, Mason AJ, et al. Comparative venom-gland transcriptomics and venom proteomics of four sidewinder rattlesnake (Crotalus cerastes) lineages reveal little differential expression despite individual variation. Sci Rep. 2018;8(1):15534. https://doi.org/10.1038/s41598-018-33943-5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Pahari S, Bickford D, Fry BG, Kini RM. Expression pattern of three-finger toxin and phospholipase A2 genes in the venom glands of two sea snakes, Lapemis curtus and Acalyptophis peronii: comparison of evolution of these toxins in land snakes, sea kraits and sea snakes. BMC Evol Biol. 2007;7(1):175. https://doi.org/10.1186/1471-2148-7-175.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Rasmussen AR. Sea snakes. In: Carpenter KE, Niem VH, editors. FAO species identification guide for fishery purposes//The living marine resources of the Western Central Pacific, vol. 6. Rome: Food and Agriculture Organization of the United Nations; 2001. p. 3988–4008.

    Google Scholar 

  21. Wang N-P, Li Q-B, Li B-Y, Li Z-Y, Li H-P, Tang S-X, et al. An epidemiological study on the snakebites in Guangxi proince, China in 1990 (in Chinese). J Snake. 1993;5(1):10–7.

    Google Scholar 

  22. Cao NV, Tao NT, Moore A, Montoya A, Rasmussen AR, Broad K, et al. Sea snake harvest in the gulf of Thailand. Conserv Biol. 2014;28(6):1677–87. https://doi.org/10.1111/cobi.12387.

    Article  PubMed  Google Scholar 

  23. Milton DA, Fry GC, Dell Q. Reducing impacts of trawling on protected sea snakes: by-catch reduction devices improve escapement and survival. Mar Freshw Res. 2009;60(8):824–32. https://doi.org/10.1071/MF08221.

    Article  Google Scholar 

  24. Phillips CM. Sea snake envenomation. Dermatol Ther. 2002;15(1):58–61. https://doi.org/10.1046/j.1529-8019.2002.01504.x.

    Article  Google Scholar 

  25. Tu AT, Fulde G. Sea snake bites. Clin Dermatol. 1987;5(3):118–26. https://doi.org/10.1016/S0738-081X(87)80018-4.

    Article  CAS  PubMed  Google Scholar 

  26. Ali SA, Alam JM, Abbasi A, Zaidi ZH, Stoeva S, Voelter W. Sea snake Hydrophis cyanocinctus venom. II. Histopathological changes, induced by a myotoxic phospholipase A2 (PLA2-H1). Toxicon. 2000;38(5):687–705. https://doi.org/10.1016/S0041-0101(99)00184-1.

    Article  CAS  PubMed  Google Scholar 

  27. Fry BG, Wüster W, Kini RM, Brusic V, Khan A, Venkataraman D, et al. Molecular evolution and phylogeny of elapid snake venom three-finger toxins. J Mol Evol. 2003;57(1):110–29. https://doi.org/10.1007/s00239-003-2461-2.

    Article  CAS  PubMed  Google Scholar 

  28. Yamazaki Y, Morita T. Structure and function of snake venom cysteine-rich secretory proteins. Toxicon. 2004;44(3):227–31. https://doi.org/10.1016/j.toxicon.2004.05.023.

    Article  CAS  PubMed  Google Scholar 

  29. World Health Organization. Guidelines for the production, control and regulation of snake antivenom immunoglobulins. In: WHO Expert Committee on Biological Standardization Sixty-seventh Report Annex 5 (WHO Technical Report Series, No 964). Geneva: World Health Organization of the United Nations; 2018.

    Google Scholar 

  30. White J. A clinician's guide to Australian venomous bites and stings: incorporating the updated CSL antivenom handbook. Parkville: Australia CSL Ltd.; 2013.

    Google Scholar 

  31. Zeng Z-Y, Wu Z-M, Wu X-Z. Observation of clinic efficacy of Sheduqing mixture on treatment of bites by sea snakes (in Chinese). J Snake. 2012;24(1):22–3.

    Google Scholar 

  32. Li Y, Tao H, Zhang L-M, Tuo Y. Progress in medical aid to the snakebites caused by sea snake (in Chinese). Nurs J Chin People’s Liberation Army. 2003;20(11):55–6.

    CAS  Google Scholar 

  33. Zhao E-M. Snakes of China (in Chinese). Heifei: Anhui Science and Technology Publishing House; 2006.

    Google Scholar 

  34. Anisimova M, Bielawski JP, Yang Z-H. Accuracy and power of the likelihood ratio test in detecting adaptive molecular evolution. Mol Biol Evol. 2001;18(8):1585–92. https://doi.org/10.1093/oxfordjournals.molbev.a003945.

    Article  CAS  PubMed  Google Scholar 

  35. Pla D, Rodríguez Y, Calvete JJ. Third generation antivenomics: pushing the limits of the in vitro preclinical assessment of antivenoms. Toxins. 2017;9(5):158.

    Article  PubMed Central  Google Scholar 

  36. Calvete JJ, Sanz L, Angulo Y, Lomonte B, Gutiérrez JM. Venoms, venomics, antivenomics. FEBS Lett. 2009;583(11):1736–43. https://doi.org/10.1016/j.febslet.2009.03.029.

    Article  CAS  PubMed  Google Scholar 

  37. Georgieva D, Öhler M, Seifert J, von Bergen M, Arni RK, Genov N, et al. Snake venomic of Crotalus durissus terrificus–correlation with pharmacological activities. J Proteome Res. 2010;9(5):2302–16. https://doi.org/10.1021/pr901042p.

    Article  CAS  PubMed  Google Scholar 

  38. Calvete JJ. Snake venomics: from the inventory of toxins to biology. Toxicon. 2013;75:44–62. https://doi.org/10.1016/j.toxicon.2013.03.020.

    Article  CAS  PubMed  Google Scholar 

  39. Lomonte B, Fernandez J, Sanz L, Angulo Y, Sasa M, Gutierrez JM, et al. Venomous snakes of Costa Rica: biological and medical implications of their venom proteomic profiles analyzed through the strategy of snake venomics. J Proteomics. 2014;105:323–39. https://doi.org/10.1016/j.jprot.2014.02.020.

    Article  CAS  PubMed  Google Scholar 

  40. Gao J-F, Wang J, He Y, Qu Y-F, Lin L-H, Ma X-M, et al. Proteomic and biochemical analyses of short-tailed pit viper (Gloydius brevicaudus) venom: age-related variation and composition–activity correlation. J Proteomics. 2014;105:307–22. https://doi.org/10.1016/j.jprot.2014.01.019.

    Article  CAS  PubMed  Google Scholar 

  41. Tan NH, Fung SY, Tan KY, Yap MKK, Gnanathasan CA, Tan CH. Functional venomics of the Sri Lankan Russell’s viper (Daboia russelii) and its toxinological correlations. J Proteomics. 2015;128:403–23. https://doi.org/10.1016/j.jprot.2015.08.017.

    Article  CAS  PubMed  Google Scholar 

  42. Xie B, Huang Y, Baumann K, Fry B, Shi Q. From marine venoms to drugs: efficiently supported by a combination of transcriptomics and proteomics. Mar Drugs. 2017;15(4):103. https://doi.org/10.3390/md15040103.

    Article  CAS  PubMed Central  Google Scholar 

  43. Valente RH, Guimarães PR, Junqueira M, Neves-Ferreira AGC, Soares MR, Chapeaurouge A, et al. Bothrops insularis venomics: a proteomic analysis supported by transcriptomic-generated sequence data. J Proteomics. 2009;72(2):241–55. https://doi.org/10.1016/j.jprot.2009.01.001.

    Article  CAS  PubMed  Google Scholar 

  44. Corrêa-Netto C, Junqueira-de-Azevedo ILM, Silva DA, Ho PL, Leitão-de-Araújo M, MLM A, et al. Snake venomics and venom gland transcriptomic analysis of Brazilian coral snakes, Micrurus altirostris and M. corallinus. J Proteomics. 2011;74(9):1795–809. https://doi.org/10.1016/j.jprot.2011.04.003.

    Article  CAS  PubMed  Google Scholar 

  45. Durban J, Juárez P, Angulo Y, Lomonte B, Flores-Diaz M, Alape-Girón A, et al. Profiling the venom gland transcriptomes of Costa Rican snakes by 454 pyrosequencing. BMC Genomics. 2011;12(1):259. https://doi.org/10.1186/1471-2164-12-259.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Durban J, Pérez A, Sanz L, Gómez A, Bonilla F, Rodríguez S, et al. Integrated “omics” profiling indicates that miRNAs are modulators of the ontogenetic venom composition shift in the central American rattlesnake, Crotalus simus simus. BMC Genomics. 2013;14(1):234. https://doi.org/10.1186/1471-2164-14-234.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Wagstaff SC, Sanz L, Juárez P, Harrison RA, Calvete JJ. Combined snake venomics and venom gland transcriptomic analysis of the ocellated carpet viper, Echis ocellatus. J Proteomics. 2009;71(6):609–23. https://doi.org/10.1016/j.jprot.2008.10.003.

    Article  CAS  PubMed  Google Scholar 

  48. Aird SD, Watanabe Y, Villar-Briones A, Roy MC, Terada K, Mikheyev AS. Quantitative high-throughput profiling of snake venom gland transcriptomes and proteomes (Ovophis okinavensis and Protobothrops flavoviridis). BMC Genomics. 2013;14(1):790. https://doi.org/10.1186/1471-2164-14-790.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Lynch VJ. Inventing an arsenal: adaptive evolution and neofunctionalization of snake venom phospholipase A2 genes. BMC Evol Biol. 2007;7(1):2. https://doi.org/10.1186/1471-2148-7-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Gibbs HL, Rossiter W. Rapid evolution by positive selection and gene gain and loss: PLA2 venom genes in closely related Sistrurus rattlesnakes with divergent diets. J Mol Evol. 2008;66(2):151–66. https://doi.org/10.1007/s00239-008-9067-7.

    Article  CAS  PubMed  Google Scholar 

  51. Casewell NR, Wagstaff SC, Harrison RA, Renjifo C, Wüster W. Domain loss facilitates accelerated evolution and neofunctionalization of duplicate snake venom metalloproteinase toxin genes. Mol Biol Evol. 2011;28(9):2637–49. https://doi.org/10.1093/molbev/msr091.

    Article  CAS  PubMed  Google Scholar 

  52. Rokyta DR, Wray KP, Lemmon AR, Lemmon EM, Caudle SB. A high-throughput venom-gland transcriptome for the eastern diamondback rattlesnake (Crotalus adamanteus) and evidence for pervasive positive selection across toxin classes. Toxicon. 2011;57(5):657–71. https://doi.org/10.1016/j.toxicon.2011.01.008.

    Article  CAS  PubMed  Google Scholar 

  53. Margres MJ, Aronow K, Loyacano J, Rokyta DR. The venom-gland transcriptome of the eastern coral snake (Micrurus fulvius) reveals high venom complexity in the intragenomic evolution of venoms. BMC Genomics. 2013;14(1):531. https://doi.org/10.1186/1471-2164-14-531.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Rokyta DR, Wray KP, Margres MJ. The genesis of an exceptionally lethal venom in the timber rattlesnake (Crotalus horridus) revealed through comparative venom-gland transcriptomics. BMC Genomics. 2013;14(1):394. https://doi.org/10.1186/1471-2164-14-394.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. McGivern JJ, Wray KP, Margres MJ, Couch ME, Mackessy SP, Rokyta DR. RNA-seq and high-definition mass spectrometry reveal the complex and divergent venoms of two rear-fanged colubrid snakes. BMC Genomics. 2014;15(1):1061. https://doi.org/10.1186/1471-2164-15-1061.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Jiang Y, Li Y, Lee WH, Xu X, Zhang Y, Zhao RP, et al. Venom gland transcriptomes of two elapid snakes (Bungarus multicinctus and Naja atra) and evolution of toxin genes. BMC Genomics. 2011;12(1):1. https://doi.org/10.1186/1471-2164-12-1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Aird SD, Aggarwal S, Villar-Briones A, Tin MM-Y, Terada K, Mikheyev AS. Snake venoms are integrated systems, but abundant venom proteins evolve more rapidly. BMC Genomics. 2015;16(1):647. https://doi.org/10.1186/s12864-015-1832-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Aird SD, Arora J, Barua A, Qiu L, Terada K, Mikheyev AS. Population genomic analysis of a pitviper reveals microevolutionary forces underlying venom chemistry. Genome Biol Evol. 2017;9(10):2640–9. https://doi.org/10.1093/gbe/evx199.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Sunagar K, Moran Y. The rise and fall of an evolutionary innovation: contrasting strategies of venom evolution in ancient and young animals. PLoS Genet. 2015;11(10):e1005596. https://doi.org/10.1371/journal.pgen.1005596.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Sanz L, Quesada-Bernat S, Chen P-Y, Lee C-D, Chiang J-R, Calvete JJ. Translational venomics: third-generation antivenomics of anti-siamese Russell’s viper, Daboia siamensis, antivenom manufactured in Taiwan CDC’s vaccine center. sTrop Med Infect Dis. 2018;3(2):66. https://doi.org/10.3390/tropicalmed3020066.

    Article  Google Scholar 

  61. Al-Shekhadat RI, Lopushanskaya KS, Segura Á, Gutiérrez JM, Calvete JJ, Pla D. Vipera berus berus venom from Russia: venomics, bioactivities and preclinical assessment of microgen antivenom. Toxins. 2019;11(2):90. https://doi.org/10.3390/toxins11020090.

    Article  CAS  PubMed Central  Google Scholar 

  62. Deka A, Gogoi A, Das D, Purkayastha J, Doley R. Proteomics of Naja kaouthia venom from north East India and assessment of Indian polyvalent antivenom by third generation antivenomics. J Proteomics. 2019;207:103463. https://doi.org/10.1016/j.jprot.2019.103463.

    Article  CAS  PubMed  Google Scholar 

  63. Pyron RA, Burbrink FT, Wiens JJ. A phylogeny and revised classification of Squamata, including 4161 species of lizards and snakes. BMC Evol Biol. 2013;13(1):93. https://doi.org/10.1186/1471-2148-13-93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52. https://doi.org/10.1038/nbt.1883.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323. https://doi.org/10.1186/1471-2105-12-323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Dutta S, Chanda A, Kalita B, Islam T, Patra A, Mukherjee AK. Proteomic analysis to unravel the complex venom proteome of eastern India Naja naja: correlation of venom composition with its biochemical and pharmacological properties. J Proteomics. 2017;156:29–39. https://doi.org/10.1016/j.jprot.2016.12.018.

    Article  CAS  PubMed  Google Scholar 

  67. Bradford MM. A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding. Anal Biochem. 1976;72(1–2):248–54. https://doi.org/10.1016/0003-2697(76)90527-3.

    Article  CAS  PubMed  Google Scholar 

  68. Calvete JJ, Juárez P, Sanz L. Snake venomics. Strategy and applications. J Mass Spectrom. 2007;42(11):1405–14. https://doi.org/10.1002/jms.1242.

    Article  CAS  PubMed  Google Scholar 

  69. Shan L-L, Gao J-F, Zhang Y-X, Shen S-S, He Y, Wang J, et al. Proteomic characterization and comparison of venoms from two elapid snakes (Bungarus multicinctus and Naja atra) from China. J Proteomics. 2016;138:83–94. https://doi.org/10.1016/j.jprot.2016.02.028.

    Article  CAS  PubMed  Google Scholar 

  70. Lanfear R, Calcott B, Ho SYW, Guindon S. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol Biol Evol. 2012;29(6):1695–701. https://doi.org/10.1093/molbev/mss020.

    Article  CAS  PubMed  Google Scholar 

  71. Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu C-H, Xie D, et al. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2014;10(4):e1003537. https://doi.org/10.1371/journal.pcbi.1003537.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Yang Z-H. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24(8):1586–91. https://doi.org/10.1093/molbev/msm088.

    Article  CAS  PubMed  Google Scholar 

  73. Rokyta DR, Margres MJ, Calvin K. Post-transcriptional mechanisms contribute little to phenotypic variation in snake venoms. G3-Genes|Genom|Genet. 2015;5(11):2375–82.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We thank Zhong-Yin Chen, Qi-Ze Liu, Jin-Geng Lv, Tao Jiang, Qiang Wang, Yi-Fan Zhao and Zi-Xuan Zhang for help in the laboratory.

Funding

This work was supported by grants from the Finance Science and Technology Project of Hainan Province (ZDYF2018120), Natural Science Foundation of China (31471995 and 31770428), Science and Technology Development Plan Project of Hangzhou (20180533B04). The funders had no role in the study design, data collection, analysis and interpretation, decision to publish, and preparation of the manuscript.

Author information

Authors and Affiliations

  1. Hangzhou Key Laboratory for Animal Adaptation and Evolution, College of Life and Environmental Sciences, Hangzhou Normal University, Hangzhou, 311121, Zhejiang, China

    Hong-Yan Zhao, Lin Wen, Yu-Feng Miao, Yan Sun, Yin Yin, Long-Hui Lin & Jian-Fang Gao

  2. Hainan Key Laboratory of Herpetological Research, College of Fisheries and Life Science, Hainan Tropical Ocean University, Sanya, 572022, Hainan, China

    Yu Du & Chi-Xian Lin

  3. MOE Key Laboratory of Utilization and Conservation for Tropical Marine Bioresources, Hainan Tropical Ocean University, Sanya, 572022, Hainan, China

    Yu Du, Chi-Xian Lin & Xiang Ji

  4. Jiangsu Key Laboratory for Biodiversity and Biotechnology, College of Life Sciences, Nanjing Normal University, Nanjing, 210023, Jiangsu, China

    Yu Du & Xiang Ji

  5. College of Life and Environmental Sciences, Wenzhou University, Wenzhou, 325035, Zhejiang, China

    Xiang Ji

Authors
  1. Hong-Yan Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Lin Wen
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Yu-Feng Miao
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yu Du
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Yan Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Yin Yin
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Chi-Xian Lin
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Long-Hui Lin
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Xiang Ji
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Jian-Fang Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

JFG and XJ conceived and designed the study. YD, CXL and LHL conducted fieldwork. HYZ, LW, YFM, YS, YY and JFG conducted laboratory and data analysis. HYZ, XJ and JFG wrote and edited the manuscript. All authors reviewed and contributed to editing of the manuscript and approved of its final publication.

Corresponding authors

Correspondence to Xiang Ji or Jian-Fang Gao.

Ethics declarations

Ethics approval and consent to participate

Snakes were collected under the permit (2017-09-12) issued by Hainan Tropical Ocean University, and our experimental procedures was approved by the Animal Research Ethics Committee of Hangzhou Normal University (AREC2019109).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

FPKM% of protein family in the venom gland transcriptome of the spine-bellied sea snake (H. curtus) from the South China Sea. Table S2. Amino acid sequences translated from toxin unigenes in the venom gland transcriptome of the spine-bellied sea snake (H. curtus) from the South China Sea. Sequences with partial CDS are indicated in italic. Table S3. Sequences used for positive selection detecting. Table S4. Identification of protein bands from whole venom separated by SDS-PAGE.

Additional file 2: Figure S1.

Original images of the SDS-PAGE profiles of venom fractions separated by RP-HPLC (1–11 peak number). Figure S2. Original images of SDS-PAGE profiles of whole venom protein (left panel) and cross-reaction between H. curtus venom and commercial antivenoms by western blotting (right panel) (Hcu: H. curtus).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhao, HY., Wen, L., Miao, YF. et al. Venom-gland transcriptomic, venomic, and antivenomic profiles of the spine-bellied sea snake (Hydrophis curtus) from the South China Sea. BMC Genomics 22, 520 (2021). https://doi.org/10.1186/s12864-021-07824-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-021-07824-7

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us