Profiling the venom gland transcriptomes of Costa Rican snakes by 454 pyrosequencing

Abstract

Background: A long term research goal of venomics, of applied importance for improving current antivenom therapy, but also for drug discovery, is to understand the pharmacological potential of venoms. Individually or combined, proteomic and transcriptomic studies have demonstrated their feasibility to explore in depth the molecular diversity of venoms. In the absence of genome sequence, transcriptomes represent also valuable searchable databases for proteomic projects.

Results: The venom gland transcriptomes of 8 Costa Rican taxa from 5 genera (Crotalus, Bothrops, Atropoides, Cerrophidion, and Bothriechis) of pitvipers were investigated using high-throughput 454 pyrosequencing. 100,394 out of 330,010 masked reads produced significant hits in the available databases. 5.165,220 nucleotides (8.27%) were masked by RepeatMasker, the vast majority of which corresponding to class I (retroelements) and class II (DNA transposons) mobile elements. BLAST hits included 79,991 matches to entries of the taxonomic suborder Serpentes, of which 62,433 displayed similarity to documented venom proteins. Strong discrepancies between the transcriptome-computed and the proteome-gathered toxin compositions were obvious at first sight. Although the reasons underlaying this discrepancy are elusive, since no clear trend within or between species is apparent, the data indicate that individual mRNA species may be translationally controlled in a species-dependent manner. The minimum number of genes from each toxin family transcribed into the venom gland transcriptome of each species was calculated from multiple alignments of reads matched to a full-length reference sequence of each toxin family. Reads encoding ORF regions of Kazal-type inhibitor-like proteins were uniquely found in Bothriechis schlegelii and B. lateralis transcriptomes, suggesting a genus-specific recruitment event during the early-Middle Miocene. A transcriptome-based cladogram supports the large divergence between A. mexicanus and A. picadoi, and a closer kinship between A. mexicanus and C. godmani.

Conclusions: Our comparative next-generation sequencing (NGS) analysis reveals taxon-specific trends governing the formulation of the venom arsenal. Knowledge of the venom proteome provides hints on the translation efficiency of toxin-coding transcripts, contributing thereby to a more accurate interpretation of the transcriptome. The application of NGS to the analysis of snake venom transcriptomes, may represent the tool for opening the door to systems venomics.

Figure 1

Gene Ontology annotation of the non-toxin Serpentes reads according to their presumed biological process (panel A) and molecular function (panel B). The figure represents the combination of the reads from all eight species. However, each species transcriptome exhibited similar relative expression levels of GO-annotated non-toxin transcript classes.

Figure 2

Calculation of the minimum number of gene copies. Multiple alignment of the top six SVMP transcripts of B. asper (Car) (Additional file 1: Table S6) using the sequence (top) of the most similar database-annotated toxin sequence as template. Each line of the multiple sequence alignment displays a distinct set of contig(s), comprised by a unique set of reads indicated in parentheses (see also Additional file 1: Table S5). Since the short average length of the reads and the low coverage of reads per contig prevented the assemblage of reliable gene sequences, each line of the alignment corresponds to at least a distinct gene of the SVMP multigene family translated into the venom gland transcriptome of B. asper (car).

Figure 3

Cartesian graph and corresponding chart pie displaying the uneven distribution of the number of reads per contig among the 29 SVMP genes identified in the venom gland transcriptome of B. asper (car).

Figure 4

Transcriptomes versus proteomes. Comparison of the protein composition of the venom of Costa Rican snakes reported from proteomic analysis (chart pies labelled "a") (Additional file 1: Table S4)^49-52 and predicted from their venom gland transcriptomes (this work). Chart pies "b" display the relative occurrence of ORF-coding reads listed in Additional file 1: Table S4 normalized for the full-length DNA sequence of a canonical member of the protein family (%mol). Panels c show the relative abundance (mol%) of toxin families as in panels "b" but computing only toxins previously identified in the venom proteome^49-52. Chart pies depicted in panels d show the relative composition (reads%) of all venom protein family hits in each of the Costa Rican snake venom gland transcriptome (Table 2). Protein family names are abbreviated as in Table 1.

Figure 5

Principal Component Analysis (PCA) of the Costa Rican snake venom gland transcriptomes (A, B) and the corresponding proteomes (C, D). Panels A and C show, respectively, the contributions to PC1 and PC2 of the different toxin families of the transcriptomes and the proteomes. Panels B and D, score plots displaying the segregation of the transcriptomes (B) and the proteomes (D) into different categories. In B, PC1 and PC2 contribute equally and together explain 65% of the observed transcriptome variability; in D, PC1 and PC2 explain 70% of the variability among proteomes.

Figure 6

Cladogram of phylogenetic alliances among Central American snakes inferred from comparison of concatenated consensus sequences for SVMP, SP, BPP, LAO and PLA₂ gathered from analysis of their venom gland transcriptomes. Numbers at branching points indicate the degree of sequence divergence (0.1 = 10% sequence divergence).