Venom of the Annulated Sea Snake Hydrophis cyanocinctus: A Biochemically Simple but Genetically Complex Weapon
by
, and
Hong-Yan Zhao
1, 
Yan Sun
1,
Yu Du
2,3,4,
Jia-Qi Li
4,
Jin-Geng Lv
2,3,
Yan-Fu Qu
4,
Long-Hui Lin
1,
Chi-Xian Lin
2,3,*,
Xiang Ji
3,4,5,* and
Jian-Fang Gao
1,* 1
Hangzhou Key Laboratory for Animal Adaptation and Evolution, College of Life and Environmental Sciences, Hangzhou Normal University, Hangzhou 311121, China
2
Hainan Key Laboratory of Herpetological Research, College of Fisheries and Life Science, Hainan Tropical Ocean University, Sanya 572022, China
3
MOE Key Laboratory of Utilization and Conservation for Tropical Marine Bioresources, Hainan Tropical Ocean University, Sanya 572022, China
4
Jiangsu Key Laboratory for Biodiversity and Biotechnology, College of Life Sciences, Nanjing Normal University, Nanjing 210023, China
5
College of Life and Environmental Sciences, Wenzhou University, Wenzhou 325035, China
*
Authors to whom correspondence should be addressed.
Toxins 2021, 13(8), 548; https://0-doi-org.brum.beds.ac.uk/10.3390/toxins13080548
Submission received: 4 July 2021
/
Revised: 27 July 2021
/
Accepted: 30 July 2021
/
Published: 6 August 2021
(This article belongs to the Section Animal Venoms)
Abstract
:Given that the venom system in sea snakes has a role in enhancing their secondary adaption to the marine environment, it follows that elucidating the diversity and function of venom toxins will help to understand the adaptive radiation of sea snakes. We performed proteomic and de novo NGS analyses to explore the diversity of venom toxins in the annulated sea snake (Hydrophis cyanocinctus) and estimated the adaptive molecular evolution of the toxin-coding unigenes and the toxicity of the major components. We found three-finger toxins (3-FTxs), phospholipase A2 (PLA2) and cysteine-rich secretory protein (CRISP) in the venom proteome and 59 toxin-coding unigenes belonging to 24 protein families in the venom-gland transcriptome; 3-FTx and PLA2 were the most abundant families. Nearly half of the toxin-coding unigenes had undergone positive selection. The short- (i.p. 0.09 μg/g) and long-chain neurotoxin (i.p. 0.14 μg/g) presented fairly high toxicity, whereas both basic and acidic PLA2s expressed low toxicity. The toxicity of H. cyanocinctus venom was largely determined by the 3-FTxs. Our data show the venom is used by H. cyanocinctus as a biochemically simple but genetically complex weapon and venom evolution in H. cyanocinctus is presumably driven by natural selection to deal with fast-moving prey and enemies in the marine environment.
Key Contribution: We found an apparent discordance in venom composition between protein and mRNA and thus concluded that H. cyanocinctus uses venom as a biochemically simple but genetically complex weapon. The toxicity of H. cyanocinctus venom is largely determined by the 3-FTx family. Nearly half of the toxin-coding unigenes were found to have undergone positive selection; thus, the venom evolution in H. cyanocinctus is presumably driven by natural selection to deal with fast-moving prey and enemies in the marine environment.
1. Introduction
Snakebite envenomation, as a neglected tropical disease, affects at least 1.8–2.7 million people, with upper estimates of 81,000–138,000 deaths annually in the world [1]. Many people are impressed with snakebite envenomation occurring on land, where most of the venomous snakes are, but they know little about the envenomation that occurs in the sea, largely due to low encounter rates with sea snakes and their less aggressive performance [2,3,4]. However, the snakebite envenomation caused by sea snakes should not be underestimated. As the largest group of marine reptiles, sea snakes are widely distributed in many tropical and subtropical waters across the Indo-Pacific Ocean [5,6], with relatively large population sizes in several regions [7,8,9,10,11]. There is no consensus on sea snake conservation around the world. Harvesting and trading of sea snakes, driven by the use of their skin, meat and blood for food, medicines and leather, are common in some Asian countries and regions [7,12,13] and over half of sea snakebites occur when fishermen handle sea snakes in fishing nets [14,15]. The neurotoxicity of sea-snake venoms is strong; thus, although the incidence of snakebites caused by sea snakes is relatively low, the death rate can be as high as 50% if the victims are not urgently treated with the correct steps [16,17].
Compared with their terrestrial relatives, sea snakes tend to simplify venom components to the extreme at the protein level [18]. Sea-snake venoms are mainly comprised of three-finger toxins (3-FTxs) and phospholipase A2 (PLA2) [19,20,21,22,23,24], with the short- and long-chain neurotoxins of 3-FTxs being the major and strong neurotoxic components [25,26]. Most sea snakes have to deal with fast-moving prey and, as such, the lethal neurotoxic venom can help them to quickly paralyze or kill prey almost instantaneously without letting them escape and avoid hurting themselves in the prey struggle [18]. Thus, from an evolutionary perspective, the venom system in sea snakes can enhance their secondary adaption to the marine environment. Recent studies on sea snakes of the genus Hydrophis show that venom toxins are far more diverse at the mRNA level than at the protein level [27,28,29]. Further elucidation of the diversity of venom toxins at both protein and mRNA levels would therefore facilitate our understanding of adaptive radiation of sea snakes, especially the taxonomic groups, such as the genus Hydrophis, with high speciation rates [30].
Hydrophis cyanocinctus used to be a dominant species in the southeast costal area of China and the South China Sea [11,31,32,33]. However, wide populations of H. cyanocinctus have declined rapidly due to overexploitation. Fortunately, all sea snakes that can be found in Chinese waters, including H. cyanocinctus, are in the newly published List of Key Protected Wild Animals in China as protected animals of national class II, meaning that China has banned the capture of sea snakes. In order to better understand and protect H. cyanocinctus, we need to know more about the diversity, evolution and function of its venom. The profiles of H. cyanocinctus venom at the omics level remain a sparsely studied area, with only two previous studies on the proteomic diversity of venom samples from Haikou, China [21], and Hara, Iran [19,21]. Here, we investigated the venom composition profile of H. cyanocinctus from the Xisha (Paracel) Islands in the South China Sea (Hainan, China) using a classical strategy for venomic analysis [34] and the toxin-coding genes profile in the venom-gland transcriptome with a de novo NGS technique. We also tested the positive selection based on the assembled toxin-coding genes with full CDS and determined the median lethal doses of both crude venom and some major venom components using ICR mice. Specifically, we aimed to clarify the diversity and correlation in abundance of venom toxins at both protein and mRNA levels, the strength of natural selection on venom toxins and the contribution of major venom components to the toxicity of the whole venom.
2. Results
2.1. Venomic Profile
Twenty protein fractions (chromatographic peaks) were distinguished in venom by RP-HPLC and 26 protein bands from SDS-PAGE were identified by MS analysis (Figure 1A and Table 1). These protein bands could be classified into three toxin families: 3-FTx (58.09%), PLA2 (40.11%) and CRISP (1.80%) (Figure 1B and Table 1). Moreover, chromatographic peaks 1–9 were identified as 3-FTx and accounted for 54.94% of the total venom protein, peaks 10–19 were comprised of a high abundance of PLA2 (40.11%) and a relatively low abundance of 3-FTx (3.15%) and peak 20 contained only CRISP. Collectively, long-neurotoxin (LNX, 38.90%) was the predominant component in the 3-FTx family and short-neurotoxin (SNX, 19.19%) and acidic (17.53%) and basic (22.58%) PLA2s differed slightly from each other in expression abundance. Among four RP-HPLC factions of H. cyanocinctus venom, the SNX (i.p. 0.09 μg/g) and LNX (i.p. 0.14 μg/g) both presented fairly low LD50s for mice, whereas the basic and acidic PLA2s did not cause any death to the mice at the upper doses of 0.6 and 2.0 μg/g, respectively (Table 2). It is thus clear that the SNX is the most toxic component in H. cyanocinctus venom and that both SNX and LNX are far more toxic than PLA2s. Due to the high abundance of neurotoxins (SNX and LNX), H. cyanocinctus venom displayed relatively low LD50 for ICR mice (crude venom i.p. 0.26 μg/g and whole venom protein i.p. 0.16 μg/g). Another three major fractions (peaks 9, 16 and 17) were also determined to not cause death at a dose of 0.2 μg/g in mice; however, they were not further employed to determine the lethality with any other doses and thus are not displayed in the final results.
2.2. Transcriptome Assembly and Diversity of Toxin-Coding Unigenes
In the present study, an Illumina Hiseq 2500 de novo NGS platform generated 109,660,888 pairs of raw reads from the venom-gland transcriptome, and 106,472,772 pairs of these data filtered as clean reads were assembled into 172,458 unigenes (N50/N90 = 1939/314) using Trinity. A total of 93,350 unigenes (N50/N90 = 2330/638) were clustered from these transcripts using Corset and 49,458 unigenes passing the quality filter (FPKM > 1) were used for the further analysis. These unigenes were finally categorized as toxins (59), non-toxins (37,405) and unidentified components (11,994), based on the similarity alignment against the NCBI NT/NR and Uniprot databases (strictly limited to “Serpentes”) using either BLAST or Diamond. Alternatively, the toxins, non-toxins and unidentified components accounted for 85.6%, 13.5% and 0.9% of the total abundance in FPKM of the H. cyanocinctus transcriptome (Figure 2 and Figure 3A). Compared with the non-toxins (8.42 FPKM/unigene), the toxins were expressed at an extremely high redundancy (159,248.25 FPKM/unigene). According to the annotated genes, 59 toxin-coding unigenes could be classified into 24 protein families (Figure 2 and Figure 3A and Supplementary Materials, Tables S1 and S2). The PLA2 (28.97%) and 3-FTx (69.76%) families constituted the most abundant components in the toxin-coding unigenes. The remaining 22 toxin families were expressed in low abundances with a total FPKM of 1.27%, including C-type lectin (CTL), CRISP, kunitz, waprin, C-type natriuretic peptide (C-NP), snake venom metalloproteinase (SVMP), cysteine-type inhibitor (cystatin), hyaluronidase (HA), 5′-nucleotidase (5′NT), aminopeptidase (AP), phospholipase B (PLB), vascular endothelial growth factor (VEGF), dipeptidyl peptidase 4 (DPP IV), cysteine-rich with EGF-like domain (CREGF), glutaminyl-peptide cyclotransferases (QC), phosphodiesterase (PDE), PLA2 inhibitor, venom factor (VF), acid phosphomonoesterase (APME), nerve growth factor (NGF), metalloproteinase inhibitor (MP inhibitor) and l-amino acid oxidase (LAAO). Tables S1 and S2 in the Supplementary Materials show detailed information on all the toxin-coding unigenes arranged according to protein family. Moreover, an overview of the venom-gland transcriptomic profiles revealed by de novo NGS in sea snakes in the current study and three recent studies is further presented in Figure 3.
2.3. Correlation between Translational and Transcriptional Abundances of Toxins
To evaluate the correlation between translational and transcriptional abundances of toxins, the MS data of protein bands were further assigned to 11 genes from an in-house database constructed from toxin-coding transcripts (Table 1). Considering that each of the five protein bands, including three 3-FTxs, one PLA2 and one CRISP, could be matched with 2–3 toxin-coding transcripts due to having the same scores, the abundance of each toxin at the protein level was divided equally among these transcripts, although this may have arbitrarily reflected the protein abundance translated by these transcripts. The results indicated that the abundance of these 11 toxins at the protein level was strongly correlated with that at the mRNA level (ρ = 0.73, p = 0.01) according to Spearman’s rank correlation analysis, and this was validated by the linear regression analysis with a relatively high Pearson’s correlation coefficient (R = 0.84, p = 0.001). Furthermore, the transcript abundances could explain the majority of variation in protein abundance (R2 = 0.71, p = 0.001) (Figure 4).
2.4. Positive Selection in Evolutionary Adaption
Either the codeml or yn00 programs in PAML 4.8 were used to evaluate the positive selection in evolutionary adaption according to the full CDS homologs aligned to toxin-coding unigenes from H. cyanocinctus. Eighteen of the fifty-nine unigenes were excluded from the analysis due to a lack of full CDS homologs. Twenty-two codeml tests were performed on 27 full-length unigenes from 18 toxin families, with a significance level of 0.002 following Bonferroni correction (Table 3 and Table S3). Furthermore, the selection of 14 full-length unigenes from 10 toxin families together with their homologs was then directly analysed using yn00 (Table S4).
The 3-FTx (SNX), CRISP and NGF families were comprised of 2–3 unigenes and the sequence divergence in each family was lower than 10%. Thus, the unigenes from the same family were combined for codeml analysis with their homologs. The results rejected M1 (null hypothesis model) in favour of M2 (positive selection model), with all p < 0.05 after Bonferroni correction. For these tests, 28% (SNX), 22% (CRISP) and 15% (NGF) of codon sites exhibited positive selection of 4.74 ≤ ω ≤ 9.84. M0 showed that all sequence sites and branches in these unigenes exhibited an average strength of selection of 1.30 ≤ ω ≤ 2.79. The nucleotide sequence divergence between two PLA2 unigenes was >10%, and these unigenes were therefore analysed separately with their homologs. M1 could be easily rejected in favour of M2 with p < 0.001 following Bonferroni correction. Moreover, 32–38% of codon sites exhibited positive selection of 4.16 ≤ ω ≤ 8.00. M0 indicated that all the sequence sites and branches in these two PLA2 unigenes had an average strength of selection of 1.57 ≤ ω ≤ 1.95.
The nucleotide sequence divergence between two 5′NT unigenes was <10%, and these unigenes were therefore combined for positive selection analysis with their homologs. M1 could not be rejected with p = 0.12. M0 indicated that all the sequence sites and branches in these two unigenes had an average strength of selection of ω = 0.40. Three CTL unigenes presented relatively large divergences (>10%) in their nucleotide sequences; therefore, these unigenes were analysed separately with their homologs. However, M1 could not be rejected in these unigenes, with p > 0.05 in all cases. M0 indicated that all sequence sites and branches in three CTLs had an average strength of selection of 0.52 ≤ ω ≤ 0.62. Two kunitz unigenes were divided into two groups due to the large nucleotide sequence divergence. M1 could only be rejected in kunitz (2) in favour of M2, with p < 0.001 after Bonferroni correction. Thirty-eight percent of codon sites exhibited positive selection of ω = 7.95. M1 could not be rejected in kunitz (1), with p > 0.05 after Bonferroni correction. M0 showed that all the sequence sites and branches in these two unigenes had an average strength of selection of 0.33 ≤ ω ≤ 3.14.
For those unigenes subjected to condeml tests, six (cystatin, DPP IV, LAAO, PLB, SVMP and VF families) easily rejected M1, and they accepted M2 with p < 0.01 in all cases after Bonferroni correction. Additionally, 2–19% of codon sites were estimated to have positive selection of 3.82 ≤ ω ≤ 11.02. M0 indicated that all the sites and branches of the six unigenes had selection strengths of 0.43 ≤ ω ≤ 1.38. In contrast, tests for the AP, HA, PLA2 inhibitor, QC and VEGF unigenes could not reject M1, with all p > 0.05 after Bonferroni correction; M0 indicated a value for ω of 0.21–0.82.
Of the 14 remaining unigenes from 10 protein families, only one homolog with high similarity could be assigned to each sequence. In each pair of sequences, we found dN of 0–0.137 and dS of 0.009–0.347 for these unigenes from the venom-gland transcriptome. We calculated that dN/dS = 6.540, 3.673 and 4.251 for three LNX and 1.032 for one SVMP unigenes, from which we inferred that these sequences probably underwent positive selection. Moreover, 10 unigenes from eight protein families, including 5′NT, aminopeptidase, CREGF, cystatin, PDE, PLA2 inhibitor, MP inhibitor and waprin, were considered to exhibit purifying selection with dN/dS ranging from 0–0.529.
3. Discussion
As much attention has been paid to the composition of snake venom at the omics level, the strategy combining venomic analysis with de novo NGS analysis has been developed as an important means to unravel variation in snake venom and also provided an important reference for investigating the function of snake venom, the clinical symptoms of snakebites and the development of novel drugs. Here, we deployed a classic snake venomic workflow combined with RP-HPLC, SDA-PAGE and MS analyses [34] to decomplex the venomic profile of a pooled crude venom sample of H. cyanocinctus. With regard to the relative abundance of toxins at the protein level, 3-FTx and PLA2 could be defined as predominant toxin families, and this was similar to two previous venomic studies using samples from populations from Hara, Iran [19], and Haikou, Hainan, China [21]. However, the Xisha population exhibited the lowest abundance of 3-FTx but the highest abundance of PLA2 among the three populations and a comparatively low divergence between these two toxin families. Specifically, the expression of acidic and basic PLA2s was approximately 1.9- and 1.8-fold higher in the Xisha population than in the Haikou population, respectively. Furthermore, the ratio of SNX to LNX (1:2) in the Xisha population was opposite to that in the Hara (1.8:1) and Haikou (2.5:1) populations, and such an inconsistent divergence has also been found in H. curtus venoms [20,22,24]. This suggests that the potential symptoms of snakebites caused by H. cyanocinctus might vary among populations. Moreover, the lethality analysis indicated that the toxicity of H. cyanocinctus venom was largely determined by the 3-FTx family. The neurotoxins of H. cyanocinctus venom were less abundant in this study (58.09%) than in an earlier one using samples from the Hara population (81.1%; [19]).This finding explains why the venom from the Xisha population is less toxic than that from the Hara population (i.p. 0.172 μg/g). Importantly, the venom components with strong neurotoxicity in sea snakes have always been inferred to play an important role in paralyzing and killing fast-moving prey [18]; they may have evolved under the selective pressure imposed by such a diet in the marine environment. Thus, compared with the Hara population, the Xisha population might have experienced a regional diet with a relatively lower percentage of fast-moving prey and further evolved to express less abundant neurotoxins in its venom. However, this inference should be verified in further studies.
Of the studies carried out in the past decade to assess the diversity of toxin-coding genes at the venom-gland transcriptomic level in snakes, only four focused on sea snakes (Acalyptophis peronii, H. curtus and Hydrophis platurus) [27,28,29,35]. Compared with H. curtus (22 protein families, Figure 3B,C; [24,29]) and H. platurus (15 protein families, Figure 3D; [27]), H. cyanocinctus (24 protein families) displayed a relatively high diversity of toxin-coding unigenes from the venom-gland transcriptome. Similar to the venom-gland transcriptome reported previously for sea snakes [26,27,28,34], the PLA2 and 3-FTx families constituted the most abundant components in the toxin-coding unigenes of the H. cyanocinctus venom-gland transcriptome. However, the ratio of PLA2 to 3-FTx unigenes in abundance was 2.4:1 in H. cyanocinctus, which is significantly lower than that reported for H. curtus (4.1:1) [29] and H. platurus (7.6:1) [27]. Moreover, the ratio of SNX to LNX was 1:1.5 in H. cyanocinctus, which is extremely converse to the ratios reported for H. curtus (2.7:1) and H. platurus (8.7:1). However, whether this implies a lower neurotoxicity of venom in H. cyanocinctus (i.p. 0.26 μg/g crude venom in mice; current study) than in H. curtus (i.v. 0.20 μg/g; [20]) and H. platurus (i.p. 0.13 μg/g; [36]) needs to be further determined.
Concordance and discordance in venom compositions and abundance at the protein and mRNA levels always receive much scientific attention due to their role in uncovering the regulatory mechanisms at the protein/mRNA levels and the factors likely to affect the evolution and function of snake venom [37,38,39,40,41,42,43,44,45,46]. Our data show that H. cyanocinctus presents an apparent discordance in venom composition between protein (three major families) and mRNA (24 families). Due to the low abundance, some toxin families can easily be neglected at the protein level. Moreover, some toxin families with extremely low abundance might contribute little to the adaptive radiation of sea snakes and have even never been reported in closely related species. These toxin families can very likely be considered as pseudogenes that cannot be detected at the protein level. However, as our gene annotation was not based on a species-specific reference genome, whether the toxin-coding unigenes with low abundance and which are undetected at the protein level can be simply judged as pseudogenes or functional genes still needs to be verified. Concordance was verified in the abundance of 11 toxins at the protein and mRNA levels, suggesting that the post-transcriptional regulation might contribute little to the abundance of these toxins in H. cyanocinctus.
As a complicated biochemical arsenal, snake venom is believed to evolve to capture prey and attack potential enemies, thus comprising an ideal model to discuss the strength of natural selection exhibited by each toxin. Positive selection in the adaptive evolution of snake venom has been detected at both the individual toxin gene and venom-gland transcriptomic levels, with most studies conducted using site model analysis [28,47,48,49,50,51,52,53]. Here, considering that the power of the likelihood-ratio test might be reduced by the excessive sequence divergence [54], we divided the unigenes from the same toxin family into several groups at a threshold of 10% nucleotide sequence divergence and conducted codeml test separately with their homologs. Briefly, 17 of 27 toxin-coding unigenes from H. cyanocinctus were found to exhibit positive selection according to the likelihood ratio between M1 and M2, and this was verified by that between M7 and M8 (the results are not shown as they are nearly identical to those for M1 and M2). The percentage (63%, 17/27) of toxin-coding unigenes that underwent positive selection in H. cyanocictus was significantly lower than that for Crotalus adamanteus (89%, 24/27; [50]) and H. curtus (73%, 24/33; [24]), and it would even be much lower if the sequences analysed by yn00 were taken into account. Notably, all of the unigenes from the 3-FTx, PLA2 and CRISP families in H. cyanocictus, which could be detected in high abundance at the protein level and might have a practical function in predation and defence, were found to undergo positive selection. This finding suggests that the positive selection of toxin-coding unigenes in H. cyanocictus might be strongly driven by the fast-moving prey and enemies in the sea environment. Toxin-coding unigenes experiencing no positive selection might play no substantive role because H. cyanocictus has evolved to prefer a simplified diet consisting mainly of fish.
4. Conclusions
Here, we used an integrated omics strategy to investigate the diversity of venom toxins at the protein and mRNA levels in H. cyanocinctus from the South China Sea. We found an apparent discordance in venom composition between protein (three major families) and mRNA (24 families) and thus concluded that H. cyanocinctus uses venom as a biochemically simple but genetically complex weapon. The abundance of 11 toxins at the protein level was strongly correlated with that at the mRNA level, indicating that the post-transcriptional regulation contributes little to the abundance of these toxins. Nearly 51.2% (21 unigenes) of the toxin-coding unigenes with full lengths in H. cyanocinctus were found to have undergone positive selection, with the proportion being much lower than reported for other venomous snakes so far studied. Among four major venom components, the short-chain neurotoxin (SNX: i.p. 0.09 μg/g) and long-chain neurotoxin (LNX: i.p. 0.14 μg/g) of 3-FTx had fairly high toxicity (LD50) towards mice, whereas both basic and acidic PLA2s were not lethal to mice at the upper doses of 0.6 and 2.0 μg/g, respectively. Moreover, the crude venom and whole venom protein separately expressed toxicities of 0.26 and 0.16 μg/g toward mice. It is thus clear that the toxicity of H. cyanocinctus venom is largely determined by the 3-FTx family. For further recognition of the contribution of each venom component to the whole venom toxicity, the remaining RP-HPLC fractions should be employed to determine the lethality, especially for those with high abundance (e.g., peaks 9, 16 and 17).
5. Materials and Methods
5.1. Animals and Ethics
We collected six adults of the H. cyanocinctus species in the waters of the Xisha Islands in the South China Sea and transferred them to the Herpetological Research Center at Hainan Tropical Ocean University where they were maintained in a circulatory sea water system. The collection of snakes was approved by Hainan Tropical Ocean University (12 September 2017), and the experimental scheme was approved by the Animal Research Ethics Committee of Hangzhou Normal University (AREC2019109).
5.2. Isolation of Venom Proteins by RP-HPLC and SDS-PAGE
Fresh venom was extracted from each snake using a 100 μL plastic pipette, then centrifuged to remove impurities for 15 min at 10,000× g 4 °C, lyophilized, equally pooled and stored at −80 °C until use. Three milligrams of venom powder was re-dissolved in 0.1% TFA and centrifuged for 15 min at 10,000× g, 4 °C, and the supernatant was automatically loaded onto a Kromasil C18 column (250 × 4.6 mm, 5 μm particle size, 300 Å pore size; AkzoNobel, Bohus, Sweden) and separated at a flow rate of 1 mL/min using a Waters E600 HPLC system (Waters, Milford, MA, USA). The whole process was performed with a linear gradient of mobile phase A (0.1% TFA in water) and B (100% ACN): 0–15% B for 30 min, followed by 15–45% B for 120 min and 45–70% B for 20 min. Protein detection was monitored at 215 nm. The fractions were collected manually and concentrated in a Labconco CentriVap® Centrifugal Concentrator (Labconco, Kansas, MO, USA). Protein concentration was determined according to Bradford [55]. The proteins of each fraction were separated by 18% SDS-PAGE under reduced conditions, and the gels were stained in 0.2% Coomassie Brilliant Blue R-250 and imaged using a Tanon Imaging system (Tanon Science & Technology, Shanghai, China).
5.3. Identification and Relative Abundance Estimation of Venom Proteins
Protein bands in the gels were excised and split into pieces with sizes of 1×1 mm, destained with 50 mM NH4HCO3 in 50% ACN and rinsed with 100% ACN in a ThermoMixer (ThermoFisher Scientific, Waltham, MA, USA). The proteins in gels were further treated with 50 mM DTT in 50 mM NH4HCO3 and alkylated with 0.14 M IAA in 50 mM NH4HCO3, rinsed again with 100% ACN in a ThermoMixer, and then were digested with trypsin (Promega, Madison, WI, USA) for 16 h at 37 °C. The digests (peptide mixture) were collected, lyophilized and re-dissolved in 0.1% TFA, desalted and enriched using an Acclaim™ PepMap™ 100 C18 column (Trap Cartridge; 5 × 0.3 mm, 5 μm; ThermoFisher Scientific), then separated by capillary RP-HPLC using an Acclaim™ PepMap™ RSLC 100 C18 column (NanoViper; 75 μm × 15 cm, 2 μm; ThermoFisher Scientific) with mobile phase A of 0.1% FA in water and B of 20% ACN in 0.1% FA as follows: 4% B for 3 min, 4–50% B for 47 min, 50–99% B for 4 min and 99% B for 6 min. Peptide eluents were subjected to a Q Exactive Orbitrap platform (ThermoFisher Scientific) according to the manufacturer’s instructions. The original MS/MS spectra were processed using Xcalibur software, and the sequence similarity was conducted using PEAKS X based on the UniProt database (strictly limited to “Serpentes”) or an in-house database (toxin transcripts from H. cyanocinctus venom-gland transcriptome). The mass tolerance was set at 0.1 Da; carbamidomethyl (C) and oxidation (M) were set as fixed and variable modifications, respectively.
The integration of the fractions in a chromatogram and the densitometry of the protein bands in an electropherogram were used to estimate the relative abundance of the venom composition [34,56]. Briefly, the relative abundance of each fraction was defined as the peak area measured by Empower software. When the fraction contained only one protein band, the relative abundance was directly defined from the peak area measurement, whereas for fractions containing more than one protein band, the relative abundance of each band was estimated by densitometry using Tanon-3500R software (Tanon Science & Technology, Shanghai, China).
5.4. Venom Gland cDNA Synthesis and Sequencing
After extraction of the venom, the snakes were kept in reptile pet terrariums and allowed to recover for four days to maximize the transcription in the venom glands, then sacrificed by injection of sodium pentobarbital (i.p. 100 mg/kg). Venom glands from both sides in each specimen were removed and dissected into pieces with sizes of < 2 mm × 2 mm, pooled and permeated with RNAlater (Qiagen, Hilden, Germany) at 4 °C overnight, and then stored at −80 °C until further use. Venom gland RNA of each specimen was extracted using TRIzol (Life Technologies, Carslbad, CA, USA), purified, concentrated and resuspended in 100 μL THE Ambion RNA storage solution (Life Technologies). The purity, integrity and concentration of RNA were evaluated using an Implen NanoPhotometer (Implen, München, Germany), Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and Qubit 2.0 fluorometer (ThermoFisher Scientific), respectively. The mixture of total RNA for sequencing was prepared by equally mixing the RNA from the above six samples, in which the RNA integrity number (RIN) ranged from 7.4 to 8.0. Subsequently, mRNA was purified and enriched from the total RNA mixture using magnetic beads attached to oligo (dT), fragmented and then used as the template to synthesize a first-strand cDNA with random hexamer primers and reverse transcriptase (RNase H). After synthesis of the second-strand cDNA with dNTPs, RNase H and DNA polymerase I, the double-stranded cDNA was absorbed and purified using AMPure XP beads (Beckman Coulter), which was conducted with end repair and ligation of a poly (A) tail and adapters. The adapter-ligated fragments of 250–300 bp in length were preferentially screened, amplified by PCR and purified using AMPure XP beads to generate a final cDNA library. After assessment of the quality using an Agilent 2100 Bioanalyzer, the cDNA library was transferred to the Illumina HiSeqTM2500 platform (Illumina, San Diego, CA, USA) for high throughout sequencing.
5.5. Transcriptome Assembly, Annotation and Quantification
Prior to sequence assembly, the raw reads containing adapters and expressing low quality (Q ≤ 20) were eliminated. The clean reads were then assembled into contigs using Trinity [57] with the following steps: a k-mer dictionary (k = 25) was initially assembled from the different clean read sets and developed into a collection of linear contigs greedily searching using Inchworm; de Bruijn graphs were constructed using Chrysalis based on the contigs that shared at least one k–1-mer and the reads spanned the junction between contigs; the de Bruijn graphs with clean reads and paired-end reads were finally reconciled and the full-length transcripts for spliced isoforms and paralogs were reconstructed and arranged using Butterfly. The longest sequences with no redundancy derived from the transcripts in each gene locus by Corset [58] and CAP3 [59] were defined as unigenes and used as references for further analyses. Gene annotation was executed by searching against the NCBI NT/NR and UniProt protein databases (strictly limited to “Serpentes”). All clean reads were assigned to the reference unigenes using RSEM, and the number of reads matched to a given unigene was defined as the readcount and converted into FPKM for estimating the abundance of unigenes [60,61].
5.6. Detection of Positive Selection
Positive selection was detected on 42 transcripts with full-length CDS from 22 toxin families. Prior to analysis, the homologs of each transcript were gathered from the NCBI nucleotide database with a threshold of 10% divergence in sequences. The sequences were aligned by MAFFT 7.313 on the basis of the AA sequence. Signal peptides, gaps and stop codons were excluded from all analyses. Then, the best-fitting model for partitions was evaluated by PartitionFinder 2.1.1 [62] based on the Akaike information criterion (AIC) and a complete search. A maximum-likelihood phylogeny was constructed by IQ-TREE 1.6.8 in PhyloSuite 1.2.2 [63]; each operation was conducted by performing 5000 ultrafast bootstraps, and the SH-aLRT branch test was conducted by performing 1000 replicates. After confirmation that the chains were converged and mixed adequately, the maximum clade credibility tree was collected as the target tree.
Positive selection in nucleotide sites based on a likelihood-ratio test was evaluated using the codeml program in PAML 4.8 [64]. Generally, to evaluate whether the sites exhibited positive selection, the difference in log-likelihoods between models M1 (the null model refers to neutral selection with dN/dS = 1 and purifying selection with dN/dS < 1) and M2 (the positive model refers to positive selection with dN/dS > 1) was determined and compared to a χ2 distribution with two degrees of freedom. Then, the initial results were further verified by comparing models M7 and M8. A single ratio for all sites with average dN/dS was calculated by the M0 model. If only one homolog with complete CDS could be matched to our transcript, then dN and dS were directly calculated with the yn00 program in PAML 4.8, then the dN/dS value was calculated manually.
5.7. Lethality
Determination of the median lethal doses (LD50) was conducted by intraperitoneal injection of H. cyanocinctus venom and the major RP-HPLC fractions into ICR mice (22–26 g, N = 4; Laboratory Animal Center of Hangzhou Normal University) of either sex as previously described [65]. The mortality was recorded over 24 h, and LD50 was calculated using the Spearman–Karber method.
5.8. Statistical Analyses
The abundance of each toxin at the protein and mRNA levels was transformed by centred log-ratio (clr) as described by Rokyta et al. [66], and the correlation between the abundances of each toxin at both levels was assessed by three coefficients (Spearman’s rank correlation coefficient, Pearson’s correlation coefficient and determination coefficient), which were calculated through non-parametric correlation and linear regression analyses using Statistica 8.0 (StatSoft, Tulsa, OK, USA). The significance level was set at α = 0.05.
Supplementary Materials
The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/toxins13080548/s1, Table S1. Amino acid sequences translated from toxin unigenes in the venom-gland transcriptome of the annulated sea snake (Hydrophis cyanocinctus) from the South China Sea. Sequences with partial CDS are indicated in italics; Table S2. FPKM% of protein family in the venom-gland transcriptome of the annulated sea snake (Hydrophis cyanocinctus) from the South China Sea; Table S3. Sequences used for positive selection detection in the codeml test; Table S4. Sequences used for positive selection detection in the yn00 test.
Author Contributions
Conceptualization, J.-F.G., H.-Y.Z., C.-X.L. and X.J.; Investigation, H.-Y.Z., Y.S., J.-Q.L., J.-G.L., Y.-F.Q., L.-H.L. and J.-F.G.; Resources, Y.D. and C.-X.L.; Writing—original draft, H.-Y.Z., X.J. and J.-F.G.; Writing—review & editing, H.-Y.Z., X.J. and J.-F.G. All authors have read and agreed to the published version of the manuscript.
Funding
This work was supported by grants from the Natural Science Foundation of China (31770428, 31471995), the Science and Technology Development Plan Project of Hangzhou (20180533B04) and the Finance Science and Technology Project of Hainan Province (ZDYF2018120) and Special Funds from the Central Finance to Support the Development of Local Universities (Project for Construction of Sea Turtle Rescue and Conservation Center in Hainan).
Institutional Review Board Statement
The collection of snakes was approved by Hainan Tropical Ocean University (12 September 2017), and the experimental scheme was approved by the Animal Research Ethics Committee of Hangzhou Normal University (AREC2019109).
Informed Consent Statement
Not applicable.
Data Availability Statement
The venom-gland transcriptome raw data can be found in the NCBI Sequence Read Archive (SRA) under the accession number PRJNA736827.
Acknowledgments
We thank Lin Wen and Xiao-Man Zhang for technical assistance.
Conflicts of Interest
The authors declare no conflict of interest.
References
- Gutiérrez, J.M.; Calvete, J.J.; Habib, A.G.; Harrison, R.A.; Williams, D.J.; Warrell, D.A. Snakebite envenoming. Nat. Rev. Dis. Primers 2017, 3, 17063. [Google Scholar] [CrossRef]
- Udyawer, V.; Goiran, C.; Shine, R.; Buijs, A. Peaceful coexistence between people and deadly wildlife: Why are recreational users of the ocean so rarely bitten by sea snakes? People Nat. 2021, 3, 335–346. [Google Scholar] [CrossRef]
- Wang, Y.-G.; Mao, J.-F.; Jin, L.-T.; Cao, C.-L.; Huang, Y.-X.; Wu, L.-Q.; Xie, J.; Pan, Q.-L. Epidemiology and prevention and treatment strategies of snake bite in coastal areas of Zhejiang. Chin. J. Crit. Care Med. 2011, 31, 541–544. [Google Scholar]
- Wang, N.; Li, Q.; Li, B.; Li, Z.; Li, H.; Tang, S.; Sawai, Y.; Kawanura, Y.; Toriba, M.; Kobayashi, T. An epidemiological study on the snakebites in Guangxi Province, China in 1990. J. Snake 1993, 5, 10–17. [Google Scholar]
- 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, 107–131. [Google Scholar]
- Rasmussen, A.R.; Murphy, J.C.; Ompi, M.; Gibbons, J.W.; Uetz, P. Marine Reptiles. PLoS ONE 2011, 6, e27373. [Google Scholar] [CrossRef] [Green Version]
- Cao, N.V.; Tao, N.T.; Moore, A.; Montoya, A.; Rasmussen, A.R.; Broad, K.; Voris, H.K.; Takacs, Z. Sea snake harvest in the gulf of Thailand. Conserv. Biol. 2014, 28, 1677–1687. [Google Scholar]
- Udyawer, V.; Read, M.; Hamann, M.; Heupel, M.R.; Simpfendorfer, C.A. Importance of shallow tidal habitats as refugia from trawl fishing for sea snakes. J. Herpetol. 2016, 50, 527–533. [Google Scholar] [CrossRef]
- Karthikeyan, R.; Vijayalakshmi, S.; Balasubramamian, T. Feeding and parturition of female annulated sea snake Hydrophis cyanocinctus in captivity. Curr. Sci. 2008, 94, 660–664. [Google Scholar]
- Liu, L.; Chen, Z.; Liu, Y.; Shi, P. Studies on the sea snakes of the middle section of Fujian’s coastal waters and their feeding habits. Acta Herpetol. Sin. 1985, 4, 341–343. [Google Scholar]
- Zhang, L.-M.; Chen, Z.-L. Medicinal value of sea snakes. J. Tradit. Chin. Med. 2002, 30, 25–26. [Google Scholar]
- Suntrarachun, S.; Chanhome, L.; Sumontha, M. Identification of sea snake meat adulteration in meat products using PCR-RFLP of mitochondrial DNA. Food Sci. Hum. Wellness 2018, 7, 170–174. [Google Scholar] [CrossRef]
- Hao, E.-W.; Deng, J.-G.; Du, Z.-C.; Hou, X.-T. Study on conventional application and biological distribution of marine traditional Chinese medicine in Guangxi. J. Chin. Med. Mater. 2016, 39, 1262–1265. [Google Scholar]
- Phillips, C.M. Sea snake envenomation. Dermatol. Ther. 2002, 15, 58–61. [Google Scholar] [CrossRef]
- Milton, D.A.; Fry, G.C.; Dell, Q. Reducing impacts of trawling on protected sea snakes: By-catch reduction devices improve escapement and survival. Mar. Freshw. Res. 2009, 60, 824–832. [Google Scholar] [CrossRef]
- Tu, A.T.; Fulde, G. Sea snake bites. Clin. Dermatol. 1987, 5, 118–126. [Google Scholar] [CrossRef]
- Zhong, W.; Xie, G.-Q.; Zeng, G.-Q. Analysis of the secondary injury caused by marine organisms in drowning person in the South China Sea. Chin. J. Naut. Med. Hyperb. Med. 2008, 15, 175–177. [Google Scholar]
- Mackessy, S.P.; Tu, A.T. Biology of the Sea Snakes and Biochemistry of Their Venoms. In Toxin-Related Diseases: Poisons Originating from Plants, Animals and Spoilage; Tu, A.T., Ed.; Oxford & IBH Publishing Co.: New Delhi, India, 1993; pp. 305–351. [Google Scholar]
- Calvete, J.J.; Ghezellou, P.; Paiva, O.; Matainaho, T.; Ghassempour, A.; Goudarzi, H.; Kraus, F.; Sanz, L.; Williams, D.J. Snake venomics of two poorly known Hydrophiinae: Comparative proteomics of the venoms of terrestrial Toxicocalamus longissimus and marine Hydrophis cyanocinctus. J. Proteom. 2012, 75, 4091–4101. [Google Scholar] [CrossRef] [PubMed]
- Tan, C.H.; Tan, K.Y.; Ng, T.S.; Sim, S.M.; Tan, N.H. 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, 3. [Google Scholar] [CrossRef] [PubMed] [Green Version]
- Wang, B.; Wang, Q.-Q.; Wang, C.; Wang, B.-L.; Qiu, L.-L.; Zou, S.-J.; Zhang, F.-H.; Liu, G.-Y.; Zhang, L.-M. A comparative analysis of the proteomes and biological activities of the venoms from two sea snakes, Hydrophis curtus and Hydrophis cyanocinctus, from Hainan, China. Toxicon 2020, 187, 35–46. [Google Scholar] [CrossRef]
- 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, 2695. [Google Scholar] [CrossRef] [PubMed] [Green Version]
- Tan, C.H.; Wong, K.Y.; Tan, K.Y.; Tan, N.H. Venom proteome of the yellow-lipped sea krait, Laticauda colubrina from Bali: Insights into subvenomic diversity, venom antigenicity and cross-neutralization by antivenom. J. Proteom. 2017, 166, 48–58. [Google Scholar] [CrossRef]
- Zhao, H.-Y.; Wen, L.; Miao, Y.-F.; Du, Y.; Sun, Y.; Yin, Y.; Lin, C.-X.; Lin, L.-H.; Ji, X.; Gao, J.-F. Venom-gland transcriptomic, venomic, and antivenomic profiles of the spine-bellied sea snake (Hydrophis curtus) from the South China Sea. BMC Genom. 2021, 22, 520. [Google Scholar] [CrossRef]
- Tan, C.H.; Tan, K.Y.; Lim, S.E.; Tan, N.H. Venomics of the beaked sea snake, Hydrophis schistosus: A minimalist toxin arsenal and its cross-neutralization by heterologous antivenoms. J. Proteom. 2015, 126, 121–130. [Google Scholar] [CrossRef]
- Laustsen, A.H.; Gutiérrez, J.M.; Rasmussen, A.R.; Engmark, M.; Gravlund, P.; Sanders, K.L.; Lohse, B.; Lomonte, B. Danger in the reef: Proteome, toxicity, and neutralization of the venom of the olive sea snake, Aipysurus laevis. Toxicon 2015, 107 Pt B, 187–196. [Google Scholar] [CrossRef] [Green Version]
- Durban, J.; Sasa, M.; Calvete, J.J. 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. [Google Scholar] [CrossRef]
- Peng, C.; Ren, J.-L.; Deng, C.; Jiang, D.; Wang, J.; Qu, J.; Chang, J.; Yan, C.; Jiang, K.; Murphy, R.W.; et al. The genome of shaw’s sea snake (Hydrophis curtus) reveals secondary adaptation to its marine environment. Mol. Biol. Evol. 2020, 37, 1744–1760. [Google Scholar] [CrossRef]
- Tan, C.H.; Tan, K.Y. De novo venom-gland transcriptomics of spine-bellied sea snake (Hydrophis curtus) from Penang, Malaysia-next-generation sequencing, functional annotation and toxinological correlation. Toxins 2021, 13, 127. [Google Scholar] [CrossRef] [PubMed]
- Sanders, K.L.; Rasmussen, A.R.; Elmberg, J.; De Silva, A.; Guinea, M.L.; Lee, M.S.Y. Recent rapid speciation and ecomorph divergence in Indo-Australian sea snakes. Mol. Ecol. 2013, 22, 2742–2759. [Google Scholar] [CrossRef] [Green Version]
- Zhao, E.-M. Snakes of China; Anhui Science and Technology Publishing House: Hefei, China, 2006. [Google Scholar]
- Zheng, J.; Wu, S.; Pang, F.; Mo, Y.; Xu, L.; Su, B.; Chang, M.; Zhuo, Z.; Lao, Z. Studies on the components of Lapemis hardwickii venom. Acta Sci. Nat. Univ. Sunyatseni 1986, 3, 122–127. [Google Scholar]
- Shi, F.; Zheng, W.-R. Determination of the venom yield and toxicity (LD50) in Hydrophis cyanocinctus from Dongshan, Fujian. J. Fujian Med. Coll. 1987, 21, 110–111. [Google Scholar]
- Calvete, J.J.; Juárez, P.; Sanz, L. Snake venomics. Strategy and applications. J. Mass Spectrom. 2007, 42, 1405–1414. [Google Scholar] [CrossRef]
- Pahari, S.; Bickford, D.; Fry, B.G.; Kini, R.M. 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, 175. [Google Scholar] [CrossRef] [Green Version]
- Lomonte, B.; Pla, D.; Sasa, M.; Tsai, W.-C.; Solórzano, A.; Ureña-Díaz, J.M.; Fernández-Montes, M.L.; Mora-Obando, D.; Sanz, L.; Gutiérrez, J.M.; 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. Proteom. 2014, 103, 137–152. [Google Scholar] [CrossRef] [PubMed] [Green Version]
- Casewell, N.R.; Wagstaff, S.C.; Wüster, W.; Cook, D.A.N.; Bolton, F.M.S.; King, S.I.; Pla, D.; Sanz, L.; Calvete, J.J.; Harrison, R.A. Medically important differences in snake venom composition are dictated by distinct postgenomic mechanisms. Proc. Natl. Acad. Sci. USA 2014, 111, 9205–9210. [Google Scholar] [CrossRef] [PubMed] [Green Version]
- Valente, R.H.; Guimarães, P.R.; Junqueira, M.; Neves-Ferreira, A.G.C.; Soares, M.R.; Chapeaurouge, A.; Trugilho, M.R.O.; León, I.R.; Rocha, S.L.G.; Oliveira-Carvalho, A.L.; et al. Bothrops insularis venomics: A proteomic analysis supported by transcriptomic-generated sequence data. J. Proteom. 2009, 72, 241–255. [Google Scholar] [CrossRef]
- Corrêa-Netto, C.; Junqueira-de-Azevedo, I.d.L.M.; Silva, D.A.; Ho, P.L.; Leitão-de-Araújo, M.; Alves, M.L.M.; Sanz, L.; Foguel, D.; Zingali, R.B.; Calvete, J.J. Snake venomics and venom gland transcriptomic analysis of Brazilian coral snakes, Micrurus altirostris and M. corallinus. J. Proteom. 2011, 74, 1795–1809. [Google Scholar] [CrossRef]
- Rodrigues, R.S.; Boldrini-França, J.; Fonseca, F.P.P.; De la Torre, P.; Henrique-Silva, F.; Sanz, L.; Calvete, J.J.; Rodrigues, V.M. Combined snake venomics and venom gland transcriptomic analysis of Bothropoides pauloensis. J. Proteom. 2012, 75, 2707–2720. [Google Scholar] [CrossRef] [PubMed]
- Margres, M.J.; McGivern, J.J.; Wray, K.P.; Seavy, M.; Calvin, K.; Rokyta, D.R. Linking the transcriptome and proteome to characterize the venom of the eastern diamondback rattlesnake (Crotalus adamanteus). J. Proteom. 2014, 96, 145–158. [Google Scholar] [CrossRef]
- Durban, J.; Juárez, P.; Angulo, Y.; Lomonte, B.; Flores-Diaz, M.; Alape-Girón, A.; Sasa, M.; Sanz, L.; Gutiérrez, J.M.; Dopazo, J.; et al. Profiling the venom gland transcriptomes of Costa Rican snakes by 454 pyrosequencing. BMC Genom. 2011, 12, 259. [Google Scholar] [CrossRef] [PubMed] [Green Version]
- Durban, J.; Pérez, A.; Sanz, L.; Gómez, A.; Bonilla, F.; Rodríguez, S.; Chacón, D.; Sasa, M.; Angulo, Y.; Gutiérrez, J.M.; 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 Genom. 2013, 14, 234. [Google Scholar] [CrossRef] [Green Version]
- Tan, C.H.; Tan, K.Y.; Fung, S.Y.; Tan, N.H. Venom-gland transcriptome and venom proteome of the Malaysian king cobra (Ophiophagus hannah). BMC Genom. 2015, 16, 1–21. [Google Scholar] [CrossRef] [Green Version]
- Aird, S.D.; Watanabe, Y.; Villar-Briones, A.; Roy, M.C.; Terada, K.; Mikheyev, A.S. Quantitative high-throughput profiling of snake venom gland transcriptomes and proteomes (Ovophis okinavensis and Protobothrops flavoviridis). BMC Genom. 2013, 14, 790. [Google Scholar] [CrossRef] [PubMed] [Green Version]
- Hofmann, E.P.; Rautsaw, R.M.; Strickland, J.L.; Holding, M.L.; Hogan, M.P.; Mason, A.J.; Rokyta, D.R.; Parkinson, C.L. 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, 15534. [Google Scholar] [CrossRef] [PubMed]
- Lynch, V.J. Inventing an arsenal: Adaptive evolution and neofunctionalization of snake venom phospholipase A2 genes. BMC Evol. Biol. 2007, 7, 2. [Google Scholar] [CrossRef] [Green Version]
- Gibbs, H.L.; 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, 151–166. [Google Scholar] [CrossRef]
- Casewell, N.R.; Wagstaff, S.C.; Harrison, R.A.; Renjifo, C.; Wüster, W. Domain loss facilitates accelerated evolution and neofunctionalization of duplicate snake venom metalloproteinase toxin genes. Mol. Biol. Evol. 2011, 28, 2637–2649. [Google Scholar] [CrossRef] [Green Version]
- Rokyta, D.R.; Wray, K.P.; Lemmon, A.R.; Lemmon, E.M.; Caudle, S.B. 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, 657–671. [Google Scholar] [CrossRef]
- Margres, M.J.; Aronow, K.; Loyacano, J.; Rokyta, D.R. The venom-gland transcriptome of the eastern coral snake (Micrurus fulvius) reveals high venom complexity in the intragenomic evolution of venoms. BMC Genom. 2013, 14, 531. [Google Scholar] [CrossRef] [Green Version]
- Rokyta, D.R.; Wray, K.P.; Margres, M.J. The genesis of an exceptionally lethal venom in the timber rattlesnake (Crotalus horridus) revealed through comparative venom-gland transcriptomics. BMC Genom. 2013, 14, 394. [Google Scholar] [CrossRef] [Green Version]
- McGivern, J.J.; Wray, K.P.; Margres, M.J.; Couch, M.E.; Mackessy, S.P.; Rokyta, D.R. RNA-seq and high-definition mass spectrometry reveal the complex and divergent venoms of two rear-fanged colubrid snakes. BMC Genom. 2014, 15, 1061. [Google Scholar] [CrossRef] [Green Version]
- Anisimova, M.; Bielawski, J.P.; Yang, Z.-H. Accuracy and power of the likelihood ratio test in detecting adaptive molecular evolution. Mol. Biol. Evol. 2001, 18, 1585–1592. [Google Scholar] [CrossRef] [Green Version]
- Bradford, M.M. A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding. Anal. Biochem. 1976, 72, 248–254. [Google Scholar] [CrossRef]
- Shan, L.-L.; Gao, J.-F.; Zhang, Y.-X.; Shen, S.-S.; He, Y.; Wang, J.; Ma, X.-M.; Ji, X. Proteomic characterization and comparison of venoms from two elapid snakes (Bungarus multicinctus and Naja atra) from China. J. Proteom. 2016, 138, 83–94. [Google Scholar] [CrossRef]
- Grabherr, M.G.; Haas, B.J.; Yassour, M.; Levin, J.Z.; Thompson, D.A.; Amit, I.; Adiconis, X.; Fan, L.; Raychowdhury, R.; Zeng, Q.; et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 2011, 29, 644–652. [Google Scholar] [CrossRef] [Green Version]
- Davidson, N.M.; Oshlack, A. Corset: Enabling differential gene expression analysis for de novoassembled transcriptomes. Genome Biol. 2014, 15, 410. [Google Scholar] [CrossRef] [Green Version]
- Huang, X.-Q.; Madan, A. CAP3: A DNA sequence assembly program. Genome Res. 1999, 9, 868–877. [Google Scholar] [CrossRef] [Green Version]
- Li, B.; Dewey, C.N. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011, 12, 323. [Google Scholar] [CrossRef] [Green Version]
- Dutta, S.; Chanda, A.; Kalita, B.; Islam, T.; Patra, A.; Mukherjee, A.K. Proteomic analysis to unravel the complex venom proteome of eastern India Naja naja: Correlation of venom composition with its biochemical and pharmacological properties. J. Proteom. 2017, 156, 29–39. [Google Scholar] [CrossRef]
- Lanfear, R.; Frandsen, P.B.; Wright, A.M.; Senfeld, T.; Calcott, B. PartitionFinder 2: New methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol. Biol. Evol. 2017, 34, 772–773. [Google Scholar] [CrossRef] [Green Version]
- Zhang, D.; Gao, F.; Jakovlic, I.; Zou, H.; Zhang, J.; Li, W.X.; Wang, G.T. PhyloSuite: An integrated and scalable desktop platform for streamlined molecular sequence data management and evolutionary phylogenetics studies. Mol. Ecol. Resour. 2020, 20, 348–355. [Google Scholar] [CrossRef] [PubMed]
- Yang, Z.-H. PAML 4: Phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 2007, 24, 1586–1591. [Google Scholar] [CrossRef] [PubMed] [Green Version]
- Gao, J.-F.; Wang, J.; He, Y.; Qu, Y.-F.; Lin, L.-H.; Ma, X.-M.; Ji, X. Proteomic and biochemical analyses of short-tailed pit viper (Gloydius brevicaudus) venom: Age-related variation and composition–activity correlation. J. Proteom. 2014, 105, 307–322. [Google Scholar] [CrossRef]
- Rokyta, D.R.; Margres, M.J.; Calvin, K. Post-transcriptional mechanisms contribute little to phenotypic variation in snake venoms. G3 Genes Genomes Genet. 2015, 5, 2375–2382. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1.
The venom proteomic profile of H. cyanocinctus. (A) Protein elution profile of venom. The proteins were separated by a C18 column as described in the Materials and Methods section. Fractions were collected and analysed by SDS-PAGE under reduced conditions. Protein bands were excised, tryptic-digested and analysed by nESI-MS/MS for their assignment to known protein families. (B) Relative abundance of venom toxin families. 3-FTx, three-finger toxin; LNX, long-chain α-neurotoxin; SNX, short-chain α-neurotoxin; PLA2, phospholipase A2; CRISP, cysteine-rich secretory protein. The details are listed in Table 1.
Figure 1.
The venom proteomic profile of H. cyanocinctus. (A) Protein elution profile of venom. The proteins were separated by a C18 column as described in the Materials and Methods section. Fractions were collected and analysed by SDS-PAGE under reduced conditions. Protein bands were excised, tryptic-digested and analysed by nESI-MS/MS for their assignment to known protein families. (B) Relative abundance of venom toxin families. 3-FTx, three-finger toxin; LNX, long-chain α-neurotoxin; SNX, short-chain α-neurotoxin; PLA2, phospholipase A2; CRISP, cysteine-rich secretory protein. The details are listed in Table 1.
Figure 2.
The venom-gland transcriptomic profile of H. cyanocinctus. The details are listed in the Supplementary Materials, Tables S1 and S2.
Figure 2.
The venom-gland transcriptomic profile of H. cyanocinctus. The details are listed in the Supplementary Materials, Tables S1 and S2.
Figure 3.
Comparison of the overall venom-gland transcriptomic profiles of three sea snakes: (A) current study; (B–D) cited and rearranged from Zhao et al. [24], Tan et al. [29] and Durban et al. [27], respectively. NC, non-conventional neurotoxin; SVSP, snake venom serine proteinase; AchE, acetylcholinesterase. Relative abundances less than 0.01% are not indicated in the toxin families.
Figure 3.
Comparison of the overall venom-gland transcriptomic profiles of three sea snakes: (A) current study; (B–D) cited and rearranged from Zhao et al. [24], Tan et al. [29] and Durban et al. [27], respectively. NC, non-conventional neurotoxin; SVSP, snake venom serine proteinase; AchE, acetylcholinesterase. Relative abundances less than 0.01% are not indicated in the toxin families.
Figure 4.
Correlation between mRNA and protein abundances of individual transcripts for each toxin family. N, number of toxin transcripts; ρ, Spearman’s rank correlation coefficient; R, Pearson’s correlation coefficient; R2, determination coefficient.
Figure 4.
Correlation between mRNA and protein abundances of individual transcripts for each toxin family. N, number of toxin transcripts; ρ, Spearman’s rank correlation coefficient; R, Pearson’s correlation coefficient; R2, determination coefficient.
Table 1.
Assignment of the chromatographic fractions and electrophoretic bands from H. cyanocinctus venom to protein families.
Table 1.
Assignment of the chromatographic fractions and electrophoretic bands from H. cyanocinctus venom to protein families.
| Peak | % | MW (kDa) | Peptide Ion | MS/MS-Derived Sequence | Protein Family/Species/Accession/Transcript ID | |
|---|---|---|---|---|---|---|
| m/z | z | |||||
| 1 | 19.19 | 9.3 | 775.3 | 2 | TTTNCAESSCYKK | 3-FTx (SNX); Hydrophis cyanocinctus; P25494; Hcy|29338, Hcy|29549, Hcy|28274 |
| 2 | 3.74 | 7.3 | 537.3 727.3 | 3 2 | XEFGCAATCPTVBR XEFGCAATCPTVB | 3-FTx (LNX); Laticauda colubrina; P0C8R8; Hcy|29273 |
| 3 | 0.27 | 7.3 | 537.3 727.3 | 3 2 | XEFGCAATCPTVBR XEFGCAATCPTVB | 3-FTx (LNX); L. colubrina; P0C8R8; Hcy|29273 |
| 4 | 0.22 | 7.3 | 805.4 727.3 | 2 2 | XEFGCAATCPTVBR XEFGCAATCPTVB | 3-FTx (LNX); L. colubrina; P0C8R8; Hcy|29273 |
| 5 | 0.23 | 7.4 | 805.4 727.3 | 2 2 | XEFGCAATCPTVBR XEFGCAATCPTVB | 3-FTx (LNX); L. colubrina; P0C8R8; Hcy|29273 |
| 6 | 1.48 | 7.0 | 537.3 727.3 | 3 2 | XEFGCAATCPTVBR XEFGCAATCPTVB | 3-FTx (LNX); L. colubrina; P0C8R8; Hcy|29273 |
| 7 | 6.17 | 7.3 | 531.9 479.9 623.2 | 3 3 2 | RXEMGCAATCPTVB XEMGCAATCPTVB SWCDAFCGSR | 3-FTx (LNX); Hydrophis hardwickii; Q8UW28; Hcy|29273 |
| 8 | 5.85 | 7.3 | 588.3 580.3 495.8 | 2 4 2 | GBVXEXGCTAB THPYBPETCPPGBNXCYBB VXEXGCTAB | 3-FTx (LNX); H. hardwickii; A3FM53; Hcy|29273 |
| 9 | 17.79 | 7.6 | 588.3 638.2 788.4 580.2 | 2 2 2 4 | GBVXEXGCTAB SWCDAFCSSR VXEXGCTABCPTVB THPYBPETCPPGBNXCYBB | 3-FTx (LNX); H. hardwickii; A3FM53; Hcy|29273, Hcy|29140 |
| 10 | 1.22 | 17.2 | 591.3 646.0 822.3 | 3 3 2 | BVCDCDVAAAECFAR NXVBFSYVXTCANHNR VCDCDVAAAECFAR | Basic PLA2; Hydrophis schistosus; P00610; Hcy|29208 |
| 11.8 | 591.3 646.0 548.6 | 3 3 3 | BVCDCDVAAAECFAR NXVBFSYVXTCANHNR VCDCDVAAAECFAR | Basic PLA2; H. schistosus; P00610; Hcy|29208 | ||
| 0.77 | 7.3 | 588.3 587.6 495.8 638.2 | 2 3 2 2 | GBVXEXGCTAB GBVXEXGCTABCPTVB VXEXGCTAB SWCDAFCSSR | 3-FTx (LNX); H. hardwickii; Q8UW29; Hcy|29273 | |
| 11 | 6.19 | 11.7 | 618.8 821.9 | 2 2 | BVCDCDVAAAB VCDCDVAAABCFAR | Basic PLA2; H. hardwickii; Q8UW30; Hcy|29208 |
| 7.3 | 591.3 548.6 646.0 | 3 3 3 | BVCDCDVAAAECFAR VCDCDVAAAECFAR NXVBFSYVXTCANHNR | Basic PLA2; H. schistosus; P00610; Hcy|29208 | ||
| 12 | 3.35 | 11.7 | 591.3 548.6 757.8 646.0 446.2 | 3 3 2 3 3 | BVCDCDVAAAECFAR VCDCDVAAAECFAR NAYNNANYNXDTB NXVBFSYVXTCANHNR XHDDCYGEAEB | Basic PLA2; H. schistosus; P00610; Hcy|29208 |
| 0.49 | 7.4 | 588.3 495.8 | 2 2 | GBVXEXGCTAB VXEXGCTAB | 3-FTx (LNX); H. hardwickii; Q8UW29; Hcy|29273 | |
| 13 | 1.51 | 12.0 | 822.3 591.3 | 2 3 | VCDCDVAAAECFAR BVCDCDVAAAECFAR | Basic PLA2; H. schistosus; P00610; Hcy|29208 |
| 0.80 | 7.3 | 588.3 638.2 587.6 | 2 2 3 | GBVXEXGCTAB SWCDAFCSSR GBVXEXGCTABCPTVB | 3-FTx (LNX); H. hardwickii; A3FM53; Hcy|29273 | |
| 14 | 0.29 | 12.0 | 591.3 548.6 | 3 3 | BVCDCDVAAAECFAR VCDCDVAAAECFAR | Basic PLA2; H. schistosus; P00610; Hcy|29123 |
| 1.10 | 7.7 | 773.4 588.3 587.6 495.8 638.2 | 3 2 3 2 2 | THPYBPETCPPGBNXCYBB GBVXEXGCTAB GBVXEXGCTABCPTVB VXEXGCTAB SWCDAFCSSR | 3-FTx (LNX); H. hardwickii; A3FM53; Hcy|29273, Hcy|29296 | |
| 15 | 2.05 | 12.3 | 455.2 | 2 | TAAXCFAR | Acidic PLA2; Tropidechis carinatus; Q45Z26; Hcy|29123 |
| 16 | 7.15 | 11.4 | 455.2 644.3 | 2 3 | TAAXCFAR DNNDECBAFXCNCDR | Acidic PLA2; T. carinatus; Q45Z28; Hcy|29123 |
| 17 | 10.02 | 10.4 | 455.8 687.7 | 2 3 | XTXYSWB CFABAPYNNBNYNXDTB | Basic PLA2; Austrelaps superbus; Q9PUH5; Hcy|29057 |
| 18 | 3.21 | 10.7 | 644.3 455.2 | 3 2 | DNNDECBAFXCNCDR AFXCNCDR | Acidic PLA2; Notechis scutatus scutatus; Q9PSN5; Hcy|29123 |
| 19 | 5.11 | 11.4 | 455.3 601.8 528.2 | 2 2 2 | TAAXCFAR GGSGTPVDEXDR AFXCNCDR | Acidic PLA2; N. s. scutatus; Q9PSN5; Hcy|29057, Hcy|29123 |
| 20 | 1.80 | 23.6 | 448.2 | 3 | CTFAHSPEHTR | CRISP; H. hardwickii; Q8UW11; Hcy|29512, Hcy|29717 |
X: Leu/Ile; B: Lys/Gln. 3-FTx, three-finger toxin; PLA2, phospholipase A2; CRISP, cysteine-rich secretory protein; LNX, long-chain α-neurotoxin; SNX, short-chain α-neurotoxin. Transcripts are listed in the Supplementary Materials, Table S1.
Table 2.
The median lethal doses (LD50) of H. cyanocinctus venom and the major fractions separated by RP-HPLC.
Table 2.
The median lethal doses (LD50) of H. cyanocinctus venom and the major fractions separated by RP-HPLC.
| Fraction | Toxin | Intraperitoneal LD50 (μg/g) * |
|---|---|---|
| - | Crude venom | 0.26 (0.23–0.30) |
| - | Whole venom protein ** | 0.16 (0.12–0.20) |
| 1 | Short neurotoxin | 0.09 (0.07–0.12) |
| 7 | Long neurotoxin | 0.14 (0.09–0.21) |
| 11 | Basic PLA2 | >0.6 |
| 18 | Acidic PLA2 | >2.0 |
*: values in parentheses are 95% confidence limits. **: Whole venom protein was defined as the protein in crude venom and quantified by the crude venom after determination of the protein concentration.
Table 3.
Summary of codeml tests for positive selection of toxins from venom-gland transcriptome in H. cyanocinctus.
Table 3.
Summary of codeml tests for positive selection of toxins from venom-gland transcriptome in H. cyanocinctus.
| Toxins (No.) | M1: Nearly Neutral | −lnL | M2: Positive Selection | −lnL | M0: ω | Δ a | p-Value b | |||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 3-FTx (1, 2, 3) | p: | 0.56 | 0.44 | 608.01 | p: | 0.45 | 0.27 | 0.28 | 588.51 | 2.79 | 39.00 | 3.40 × 10−9 * |
| ω: | 0.00 | 1.00 | ω: | 0.00 | 1.00 | 9.84 | ||||||
| 5NT (1, 2) | p: | 0.60 | 0.40 | 3831.74 | p: | 0.90 | 0.00 | 0.10 | 3829.66 | 0.40 | 4.16 | 0.12 |
| ω: | 0.00 | 1.00 | ω: | 0.21 | 1.00 | 2.57 | ||||||
| AP (1) | p: | 0.63 | 0.37 | 5224.95 | p: | 0.93 | 0.00 | 0.07 | 5221.58 | 0.40 | 6.74 | 0.03 |
| ω: | 0.00 | 1.00 | ω: | 0.22 | 1.00 | 3.95 | ||||||
| CRISP (1, 2) | p: | 0.47 | 0.53 | 2499.04 | p: | 0.40 | 0.39 | 0.22 | 2458.11 | 1.38 | 81.87 | 0.00* |
| ω: | 0.00 | 1.00 | ω: | 0.00 | 1.00 | 4.74 | ||||||
| CTL (1) | p: | 0.56 | 0.44 | 1046.59 | p: | 0.82 | 0.00 | 0.18 | 1045.29 | 0.52 | 2.60 | 0.27 |
| ω: | 0.00 | 1.00 | ω: | 0.19 | 1.00 | 2.49 | ||||||
| CTL (2) | p: | 0.60 | 0.40 | 1177.31 | p: | 0.80 | 0.07 | 0.13 | 1172.58 | 0.56 | 9.46 | 8.83 × 10−3 |
| ω: | 0.06 | 1.00 | ω: | 0.22 | 1.00 | 3.07 | ||||||
| CTL (3) | p: | 0.55 | 0.45 | 1186.50 | p: | 0.60 | 0.00 | 0.40 | 1184.23 | 0.62 | 4.55 | 0.10 |
| ω: | 0.00 | 1.00 | ω: | 0.00 | 1.00 | 1.65 | ||||||
| Cystatin (1) | p: | 0.29 | 0.71 | 979.80 | p: | 0.96 | 0.00 | 0.04 | 969.73 | 1.15 | 20.13 | 4.25 × 10−5 * |
| ω: | 0.04 | 1.00 | ω: | 0.84 | 1.00 | 11.02 | ||||||
| DPP IV | p: | 0.64 | 0.36 | 4740.05 | p: | 0.83 | 0.11 | 0.06 | 4730.48 | 0.43 | 19.12 | 7.05 × 10−5 * |
| ω: | 0.00 | 1.00 | ω: | 0.15 | 1.00 | 3.82 | ||||||
| HA | p: | 0.45 | 0.55 | 2976.76 | p: | 0.52 | 0.00 | 0.48 | 2975.36 | 0.62 | 2.81 | 0.25 |
| ω: | 0.00 | 1.00 | ω: | 0.00 | 1.00 | 1.38 | ||||||
| Kunitz (1) | p: | 0.82 | 0.18 | 1502.79 | p: | 0.82 | 0.11 | 0.08 | 1502.79 | 0.33 | 0.00 | 1.00 |
| ω: | 0.21 | 1.00 | ω: | 0.21 | 1.00 | 1.00 | ||||||
| Kunitz (2) | p: | 0.44 | 0.56 | 944.01 | p: | 0.18 | 0.44 | 0.38 | 916.66 | 3.14 | 54.70 | 1.32 × 10−12 * |
| ω: | 0.05 | 1.00 | ω: | 0.00 | 1.00 | 7.95 | ||||||
| LAAO | p: | 0.44 | 0.56 | 6029.83 | p: | 0.39 | 0.43 | 0.17 | 5970.56 | 1.03 | 118.53 | 0.00* |
| ω: | 0.00 | 1.00 | ω: | 0.00 | 1.00 | 4.00 | ||||||
| NGF (1, 2) | p: | 0.51 | 0.49 | 1366.53 | p: | 0.39 | 0.46 | 0.15 | 1342.54 | 1.30 | 47.98 | 3.81 × 10−11 * |
| ω: | 0.07 | 1.00 | ω: | 0.13 | 1.00 | 5.87 | ||||||
| PLA2 (1) | p: | 0.47 | 0.53 | 1913.94 | p: | 0.41 | 0.21 | 0.38 | 1877.93 | 1.57 | 72.01 | 2.22 × 10−16 * |
| ω: | 0.01 | 1.00 | ω: | 0.00 | 1.00 | 4.16 | ||||||
| PLA2 (2) | p: | 0.45 | 0.55 | 808.89 | p: | 0.68 | 0.00 | 0.32 | 796.33 | 1.95 | 25.12 | 3.51 × 10−6 * |
| ω: | 0.00 | 1.00 | ω: | 0.15 | 1.00 | 8.00 | ||||||
| PLA2 inhibitor (1) | p: | 0.48 | 0.52 | 1900.98 | p: | 0.50 | 0.00 | 0.50 | 1900.96 | 0.50 | 0.04 | 0.98 |
| ω: | 0.00 | 1.00 | ω: | 0.00 | 1.00 | 1.07 | ||||||
| PLB | p: | 0.56 | 0.44 | 3750.75 | p: | 0.53 | 0.45 | 0.02 | 3742.03 | 0.54 | 17.44 | 1.63 × 10−4 * |
| ω: | 0.00 | 1.00 | ω: | 0.00 | 1.00 | 6.82 | ||||||
| QC | p: | 0.78 | 0.22 | 2194.10 | p: | 0.79 | 0.00 | 0.21 | 2194.09 | 0.21 | 0.01 | 1.00 |
| ω: | 0.00 | 1.00 | ω: | 0.00 | 1.00 | 1.05 | ||||||
| SVMP (1) | p: | 0.48 | 0.52 | 8066.12 | p: | 0.35 | 0.46 | 0.19 | 7910.53 | 1.38 | 311.17 | 0.00 * |
| ω: | 0.04 | 1.00 | ω: | 0.00 | 1.00 | 5.06 | ||||||
| VEGF | p: | 0.31 | 0.69 | 1101.80 | p: | 0.87 | 0.00 | 0.13 | 1101.01 | 0.82 | 1.58 | 0.45 |
| ω: | 0.00 | 1.00 | ω: | 0.49 | 1.00 | 3.83 | ||||||
| VF | p: | 0.53 | 0.47 | 13764.79 | p: | 0.51 | 0.46 | 0.04 | 13,726.91 | 0.61 | 75.76 | 0.00 * |
| ω: | 0.00 | 1.00 | ω: | 0.00 | 1.00 | 4.96 | ||||||
*, significance at the 5% level after Bonferroni correction; a, negative twice the difference in lnL between M1 and M2; b, p-value before correction.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. |
© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
Share and Cite
MDPI and ACS Style
Zhao, H.-Y.; Sun, Y.; Du, Y.; Li, J.-Q.; Lv, J.-G.; Qu, Y.-F.; Lin, L.-H.; Lin, C.-X.; Ji, X.; Gao, J.-F. Venom of the Annulated Sea Snake Hydrophis cyanocinctus: A Biochemically Simple but Genetically Complex Weapon. Toxins 2021, 13, 548. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins13080548
AMA Style
Zhao H-Y, Sun Y, Du Y, Li J-Q, Lv J-G, Qu Y-F, Lin L-H, Lin C-X, Ji X, Gao J-F. Venom of the Annulated Sea Snake Hydrophis cyanocinctus: A Biochemically Simple but Genetically Complex Weapon. Toxins. 2021; 13(8):548. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins13080548
Chicago/Turabian StyleZhao, Hong-Yan, Yan Sun, Yu Du, Jia-Qi Li, Jin-Geng Lv, Yan-Fu Qu, Long-Hui Lin, Chi-Xian Lin, Xiang Ji, and Jian-Fang Gao. 2021. "Venom of the Annulated Sea Snake Hydrophis cyanocinctus: A Biochemically Simple but Genetically Complex Weapon" Toxins 13, no. 8: 548. https://0-doi-org.brum.beds.ac.uk/10.3390/toxins13080548
Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.
