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 20070601, expiry date: 06/2010) and N. atra antivenom (batch 20070601, 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 four 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 (s.c. 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 and Agilent 2100 Bioanalyzer system, respectively, whereas the concentration of RNA was measured using a Qubit 2.0 fluorometer. 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® UltraTM 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. High throughout sequencing of the cDNA library was performed using the Illumina HiSeqTM2500 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 based on Inchworm, Chrysalis, and Butterfly modules [57]. 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 and defined as a unigene, which was used as a reference for downstream analyses. Gene annotation was performed by BLAST searching against the NCBI NT database and Diamond 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 [58]. The number of reads matched to a given unigene was defined as the readcount, which was then calculated as FPKM [59] 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. 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 [60]. 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.
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 a-cyano-4-hydroxycinnamic acid (ABI) in 50% ACN and 0.1% TFA was added; then, the sample was dried completely and subjected to ABI 4800 Plus MALDI-TOF/TOF-MS mass spectrometer 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 according to the manufacturer’s instructions. The raw MS/MS spectra were interpreted using FlexAnalysis or Xcalibur 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. [61]and Shan et al. [62]. Briefly, the relative abundance of fractions was calculated by peak area measurement using Empower 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 [63] based on the Bayesian Information Criterion and greedy search. Then, a Bayesian phylogeny was constructed using BEAST 2.2 [64] 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 [65] 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 [32]. 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, pH8.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.
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 PierceTM ECL western blotting substrate and exposed to X-ray film. The film was scanned using a UMax2100 densitometer and analyzed with Tan4100 software.
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. [66]. 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