Comparative analyses of putative toxin gene homologs from an Old World viper, Daboia russelii

Russell’s viper shed skin and DNA isolation

Freshly shed skin of Russell’s viper from Bangalore, India was a gift from Mr. Gerry Martin. The shed skin for the entire snake was obtained, cleaned thoroughly with 70% ethanol and with nuclease-free water three times each, dried thoroughly and frozen until the time of extraction of DNA. Genomic DNA was extracted following the protocol of Fetzner Jr (1999) with modifications.

Sequencing, read processing and assembly

Illumina paired-end read libraries (100 base paired-end reads with insert size of 350 bases) were prepared following the manufacturer instructions using amplification free genomic DNA library preparation kit and sequenced using Illumina HiSeq2500 instrument. Archaeal, bacterial and human sequence contamination were removed from the Russell’s viper sequence by DeConSeq (Schmieder & Edwards, 2011) using curated and representative genomes (https://www.ncbi.nlm.nih.gov/genome/browse/reference/).

Furthermore, the sequenced reads were post-processed to remove unpaired reads and quality analysis was performed using FastQC v0.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The rd_len_cutoff option was exercised during the read assembly step to trim off the low-quality bases, since the per-base quality was found to drop below 28 after the initial 50–70 bases of the read. The Russell’s viper read libraries with 26× coverage were assembled using SOAPdenovo2 (r240) (Luo et al., 2012).

Identifying toxin gene homologs, coding regions, and predicted gene structures

The DNA sequences for 51 out of 54 venom-associated genes (Fry, 2005) from Protobothrops mucrosquamatus  were downloaded (Table 1). These were used to fish genomic scaffolds bearing highly similar sequences in Russell’s viper genome assembly, using BLAST with an E-value threshold of 10⁻³. The fished scaffolds were then anchored to the respective coding sequences from Protobothrops mucrosquamatus, using a discontiguous megaBlast, to determine the correct frame of translation and extract the complete amino acid coding sequence (CDS) corresponding to putative toxin homologs in Russell’s viper. We obtained the exon-intron structures for all the putative toxin gene homologs in Russell’s viper by aligning the CDS with gene sequences using discontiguous megaBlast and plotted using the tool GSDS2.0 (Hu et al., 2015). The sequences for the Russell’s viper putative venom gene homologs were deposited in GenBank and their accession numbers are provided in Table S1.

Comparative analysis between venom proteins and their putative homologs

We obtained the accession IDs for the major toxin families from Russell’s viper of the Indian sub-continent (Fig. S1 in Sharma et al., 2015). Their corresponding protein sequences were matched using blastp with the amino acid sequences from the putative skin homologs. For the genes covered under each family, a percent identity metric, indicative of the extent of sequence similarity between the venom proteins and their skin homologs, was estimated. Similar comparative analyses were performed for king cobra (Ophiophagus hannah) using accession IDs provided in Additional File 4 of Tan et al. (2015), and the predicted toxin protein homologs from blood of king cobra (PRJNA201683 from Vonk et al., 2013). Comparative analyses were performed using blastp, with the venom protein sequence as the query, against PRJNA201683.

Comparative analyses of putative venom protein homolog domains

The amino acid sequences of all the Russell’s viper’s putative toxin homologs were subjected to domain search using Pfam (Finn et al., 2016) (Table S2). All domain sequences were aligned using blastp to non-redundant protein sequences from 18 snake species (Table S3). We wanted to compare the gene structures between the venomous and the non-venomous animals, hence included sequence information from the members of the later group. Five domains (NGF, PDGF, Kunitz BPTI, CAP and CRISP) from four genes (NGF, VEGF, CRISP/Serotriflin, and Kunitoxin) with variability across different snake groups were used for expansive comparative analyses (Table S4). The sequences used were from viperids (taxid: 8689), elapids (taxid: 8602), colubrids (taxid: 8578), boids (taxid: 8572), acrochordids (taxid: 42164), pythonids (taxid: 34894), lizards (squamates (taxid: 8509) minus snakes (taxid: 8570), crocodiles (taxid: 51964) and testudines (taxid: 8459).

3D structure prediction of the chosen domains

Consensus sequences were determined from NGF, PDGF, Kunitz BPTI, CAP and CRISP domain alignments using Simple Consensus Maker (https://www.hiv.lanl.gov/content/sequence/CONSENSUS/SimpCon.html) for crotalines (CR), viperines (VP) and elapids. The consensus sequences were submitted to the protein fold recognition server (Kelley et al., 2015) using standard mode (http://www.sbg.bio.ic.ac.uk/phyre2/html/page.cgi?id=index). The best 3D model was further investigated by Phyre2 to analyze the structural model using various open source tools.

Shed skin yielded fairly good quality and high-molecular weight genomic DNA

Genomic DNA isolated from the shed skin of Russell’s viper was fairly intact with most of the DNA in the size range of more than 5 kbp (Fig. S1). Sequenced short reads were assembled and then used to fish the sequences for the 51 putative toxin genes in Russell’s viper (see ‘Materials and Methods’). Next, we obtained the exon-intron structures for all putative homologs in Russell’s viper by aligning the CDS with gene sequences (Fig. S2). We found the average length of the exons in those sequences to be around 190 nucleotides (nt), matching well with the lengths of other vertebrate exons (Gelfman et al., 2012).

Similarity between venom proteins and their putative skin homologs

For the Russell’s viper, we found 45–100% sequence similarity between the major venom proteins and their predicted putative skin homologs (Fig. 1). The sequences for venom nerve growth factor (VNGF) and its putative skin homolog were identical. Similarly, VEGF and CRISPs from venom gland were highly similar to their putative skin homologs (99% and 92% sequence similarity, respectively). Other proteins like, KSPI, SVSPs and PLA2 showed 79%, 74% and 61% sequence identity, respectively (Fig. 1 and Fig. S3). In order to find out whether the sequence divergence between some of the venom gland proteins and their predicted putative skin homologs was specific to Russell’s viper, we performed similar analysis using venom proteins and their blood homologs from king cobra, Ophiophagus hannah (Vonk et al., 2013). In the case of Ophiophagus hannah, the differences between toxin proteins and their blood homologs were minor for most families studied (similarity ≥75%), except for PLA2, which had a low similarity of 23% (Figs. S4 and S5).

Comparative domain analyses

Among the genes, a larger pool of sequences were available only for NGF, PDGF domain of VEGF, Kunitz_BPTI domain of Kunitoxin, CRISP and CAP domains in CRISP and Serotriflin proteins, from various snake groups (Colubridae, Boidae, Pythonidae and Acrochordidae), non-snake reptilian groups (lizards, crocodiles and Testudines), venomous invertebrates (wasps, spiders and scorpions) and venomous vertebrates (fishes and mammals). Therefore, these domains were compared with their putative homologs from Russell’s viper. Comparative domain analysis was performed for all putative toxin gene homologs (Fig. S6) across 18 snake species for those where sequence information was available (Table S3). In the case of five domains: CAP and CRISP domains of CRISP and serotriflin genes (L and AL), Kunitz BPTI of kunitoxin (S), NGF (T) and PDGF of VEGFA (AP-AR) and VEGFF (AU), we found that the maximum number of species aligned to their domain sequences. Some protein domains, the CRISP, Kunitz BPTI, guanylate CYC, PDGF of VEGFF and WAP, showed long stretches of mismatches (Fig. S6) compared with Russell’s viper sequence. Out of these, only NGF and PDGF domains of VEGF had amino acid changes specific to the members of the group Crotalinae, that were completely absent in any other group used for comparison, including in lizards, crocodiles, and turtles (Fig. S7). Specific changes in these proteins and their implications are discussed below.

The putative skin-derived NGF gene homolog in Russell’s viper is a single exon gene with a 745nt transcript coding for a 244 amino acid protein consisting of a single NGF domain (Fig. 2A). The NGF domain bears 28% sequence conservation across all the five vertebrate phyla, namely, fishes, amphibians, reptiles, aves and mammals distributed along the length of the domain (Fig. 2B). Thirty-six percent out of these residues are conserved across other venomous vertebrates (fishes, squamates and mammals) and venomous invertebrates (scorpions and wasps) (Fig. 2C). Thirteen percent of the skin-derived putative NGF domain residues are variable with respect to the domain sequence in at least one among the NGF sequences in the groups of viperids and elapids (Fig. 2D). Although several amino acids in the NGF domain in crotalines seem to have changed from the putative domain in Russell’s viper and other vipers of the group Viperinae, their function probably remains unchanged. For example, phenylalanine (F) to isoleucine (I) at position 12 and serine (S) to asparagine (N) at position 19 between the crotalines and viperines does not change the function of the amino acids (from one hydrophobic amino acid to another and from one polar amino acid to another). However, it is also true that F changing to I removes the bulky aromatic ring, whereas S could be a phosphorylated site as opposite to the N in the same position. There are others, for example, threonine (T) and glutamine (Q), at positions 67 and 68, respectively, in the NGF domain of the crotalines, which were only there in that specific group. One of those, a polar amino acid glutamine at position 68, is a very important residue as its corresponding amino acid in any of the other snakes, except in colubrids, is a hydrophobic proline.

In Russell’s viper, the putative skin-derived VEGFA gene homolog comprises five exons coding for a 652 nt long transcript and a translated protein with two domains: PDGF and VEGF-C (Fig. 3A). The PDGF domain sequence exhibits conservation in 65% of its residues across the three vertebrate phyla (reptiles, birds and mammals) (Fig. 3B). Since sequence information from fishes and amphibians were not available, we could not include those in the comparison study. Out of the conserved residues in the above said domain, 21% were also conserved in venomous vertebrates (squamates and mammals) and venomous invertebrates (wasps). Fifteen percent of the PDGF domain residues were variable in at least one of the two snake groups: viperids and elapids (Figs. 3C and 3D). Like the NGF domain, the evolution of the putative skin-derived PDGF domain in crotalines at certain amino acids is striking. For example, in the crotalines, the position 67 is a polar amino acid tyrosine (Y) while in all other reptiles, venomous invertebrates and mammals; this is primarily a hydrophobic amino acid phenylalanine (F).

The putative skin-derived Kunitoxin gene homolog in Russell’s viper is a 3.1 kb gene comprising two exons, with a transcript length of 270 nt that codes for a 44 amino acids long single Kunitz BPTI domain (Fig. 4A). About 29% of the protein domain residues are conserved across the four vertebrate phyla (amphibians, reptiles, aves and mammals) (Fig. 4B). Since sequence information from the Kunitz BPTI for fishes was not available, we could not include those in the comparison. Out of these conserved residues, 76% are conserved in venomous vertebrates (squamates and mammals) and venomous invertebrates (scorpions and wasps) (Fig. 4C) and 56% of the domain residues are variable in at least one of two snake groups (viperids and elapids) (Fig. 4D). Of the residues that are evolved in the members of Crotalinae, the second residue, a positively charged one, arginine (R) is present only in the members of Viperinae, which is replaced by a hydrophobic residue, proline (P), in the crotalines and elapids. Residues 14–18 are very polymorphic in the crotalines and elapids, but not so in the viperines.

The skin-derived putative CRISP gene homolog in Russell’s viper is a 25 kb long gene, comprises of eight exons coding for a 787 nt transcript and two translated protein domains, CAP and CRISP (Fig. 5A). The CAP domain exhibits conservation in 7% of its residues across all the five vertebrate phyla (Fig. 5B). Forty-two percent of those residues are conserved across venomous vertebrates (amphibians, squamates and mammals) and venomous invertebrates (scorpions and wasps) (Fig. 5C). In addition, there are five residues conserved across all the venomous animals (Fig. 5C). Twenty-seven percent of the CAP domain residues are variable in at least one of the three snake groups (Fig. 5D). There are several extra residues for the CAP domain in the crotalines and elapids, but not in the viperines. The conserved residues comprised mostly of Cystines and to a lesser extent Asparagines (Fig. 5E) across venomous vertebrates (squamates and mammals) (Fig. 5F).Sixty percent of the residues in the putative homolog of the CRISP domain are variable in at least one viperine or elapid member with respect to the domain sequence of Russell’s viper (Fig. 5G).

Next, we explored the role of consensus domain sequences in the putative protein homologs and the possible role of conserved amino acids in those domains across viperids and elapids. We constructed the 3D structure models using Phyre2, followed by Phyre2 investigation, for further analyses on the structural model. As evident from the analyses, amino acid residues 18–19 and 117 of the NGF domain reflected a difference in mutation sensitivity as detected by SusPect algorithm (Yates et al., 2014), especially in the elapids compared to the viperids (Fig. 6). There were differences in certain residues across these two groups. Residue 18 is Valine in the viperines and Isoleucine in the elapids; residue 19 is Serine in the viperines and Asparagine in the crotalines; and residue 117 is Threonine in the elapids and Serine in the crotalines (Fig. 6A). This might have implications in the structure of the protein as the largest pockets detected by fpocket algorithm appear to be vastly different among the crotalines, viperines and elapids for the NGF, PDGF, CAP and CRISP domains (Fig. 6). The pockets appeared smallest in all cases for the elapids (quantitatively substantial for PDGF, CAP and CRISP: one-fourth that of viperines for PDGF, and two-thirds that of viperines for CAP and CRISP), and largest in the case of viperines (Fig. 6). Minor differences in clashes were observed at residues 10, 11 and 20 of the Kunitz domain and residue 38 of this domain showed a rotamer conflict in the case of the crotalines (Fig. 6C). Similarly, residue 46 of the CAP domain and residues 4 and 31 of the CRISP domain showed rotamer conflict for the viperines (Figs. 6D and 6E). The other protein quality and functional parameters were not affected across the 3D structure models for the three snake groups (Fig. S8).