The venom-gland transcriptome of the eastern diamondback rattlesnake (Crotalus adamanteus)

Abstract

Background: Snake venoms have significant impacts on human populations through the morbidity and mortality associated with snakebites and as sources of drugs, drug leads, and physiological research tools. Genes expressed by venom-gland tissue, including those encoding toxic proteins, have therefore been sequenced but only with relatively sparse coverage resulting from the low-throughput sequencing approaches available. High-throughput approaches based on 454 pyrosequencing have recently been applied to the study of snake venoms to give the most complete characterizations to date of the genes expressed in active venom glands, but such approaches are costly and still provide a far-from-complete characterization of the genes expressed during venom production.

Results: We describe the de novo assembly and analysis of the venom-gland transcriptome of an eastern diamondback rattlesnake (Crotalus adamanteus) based on 95,643,958 pairs of quality-filtered, 100-base-pair Illumina reads. We identified 123 unique, full-length toxin-coding sequences, which cluster into 78 groups with less than 1% nucleotide divergence, and 2,879 unique, full-length nontoxin coding sequences. The toxin sequences accounted for 35.4% of the total reads, and the nontoxin sequences for an additional 27.5%. The most highly expressed toxin was a small myotoxin related to crotamine, which accounted for 5.9% of the total reads. Snake-venom metalloproteinases accounted for the highest percentage of reads mapping to a toxin class (24.4%), followed by C-type lectins (22.2%) and serine proteinases (20.0%). The most diverse toxin classes were the C-type lectins (21 clusters), the snake-venom metalloproteinases (16 clusters), and the serine proteinases (14 clusters). The high-abundance nontoxin transcripts were predominantly those involved in protein folding and translation, consistent with the protein-secretory function of the tissue.

Conclusions: We have provided the most complete characterization of the genes expressed in an active snake venom gland to date, producing insights into snakebite pathology and guidance for snakebite treatment for the largest rattlesnake species and arguably the most dangerous snake native to the United States of America, C. adamanteus. We have more than doubled the number of sequenced toxins for this species and created extensive genomic resources for snakes based entirely on de novo assembly of Illumina sequence data.

Figure 1

Merging overlapping reads. (A) Reads are slid along each other until the number of matches exceeds the significance threshold. In the example shown, the optimal overlap is 74 nucleotides (nt). (B) The quality of reads declines dramatically toward their 3’ ends, where overlap occurs if the fragment length is less than twice the read length, allowing the actual quality to be much higher than the nominal values. The example shown is the average of pairs that overlap by exactly 50 nt.

Figure 2

Domination of the C. adamanteus venom-gland transcriptome by toxin transcripts. The 123 unique toxin sequences were clustered into 78 groups with less than 1% nucleotide divergence for estimation of abundances. (A) The vast majority of the extremely highly expressed genes were toxins. The inset shows a magnification of the top 200 transcripts. (B) Expression levels of individual toxin clusters are shown with toxin classes coded by color. The toxin clusters are in the same order as in Table 3.

Figure 3

Expression levels of major classes of toxins and nontoxins. More than 60% of the total reads have been accounted for with full-length annotated transcripts. (A) The major toxin classes were the CTLs, SVSPs, MYO, and SVMPs (types II and III). (B) As expected for a protein-secreting tissue, the venom gland expresses an abundance of proteins involved in proteostasis.

Figure 4

Comparison of gene ontology (GO) results for our annotated full-length nontoxin sequences with those of the contigs from ade novo assembly with NGen. Only level 2 GO terms are shown. The distributions of GO terms are similar across data sets, suggesting that the annotated transcripts provided a comprehensive characterization of the genes expressed in the venom gland. (A) The distributions of sequences reaching various stages of identification and annotation are shown. The level 2 GO terms are shown for molecular function (B), biological process (C), and cellular component (D).

Figure 5

The biological-process GO terms identified for the 2,879 annotated full-length nontoxin sequences. Terms specific for the production, processing, and export of proteins are highlighted in black. The inset shows the low-abundance portion of the full distribution.

Figure 6

The molecular-function GO terms identified for the 2,879 annotated full-length nontoxin sequences. Terms specific for the production, processing, and export of proteins are highlighted in black. The inset shows the low-abundance portion of the full distribution.

Figure 7

The cellular-components GO terms identified for the 2,879 annotated full-length nontoxin sequences. Terms specific for the production, processing, and export of proteins are highlighted in black. The inset shows the low-abundance portion of the full distribution.