Tracking the recruitment and evolution of snake toxins using the evolutionary context provided by the Bothrops jararaca genome

Abstract

Venom is a key adaptive innovation in snakes, and how nonvenom genes were co-opted to become part of the toxin arsenal is a significant evolutionary question. While this process has been investigated through the phylogenetic reconstruction of toxin sequences, evidence provided by the genomic context of toxin genes remains less explored. To investigate the process of toxin recruitment, we sequenced the genome of Bothrops jararaca, a clinically relevant pitviper. In addition to producing a road map with canonical structures of genes encoding 12 toxin families, we inferred most of the ancestral genes for their loci. We found evidence that 1) snake venom metalloproteinases (SVMPs) and phospholipases A₂ (PLA2) have expanded in genomic proximity to their nonvenomous ancestors; 2) serine proteinases arose by co-opting a local gene that also gave rise to lizard gilatoxins and then expanded; 3) the bradykinin-potentiating peptides originated from a C-type natriuretic peptide gene backbone; and 4) VEGF-F was co-opted from a PGF-like gene and not from VEGF-A. We evaluated two scenarios for the original recruitment of nontoxin genes for snake venom: 1) in locus ancestral gene duplication and 2) in locus ancestral gene direct co-option. The first explains the origins of two important toxins (SVMP and PLA2), while the second explains the emergence of a greater number of venom components. Overall, our results support the idea of a locally assembled venom arsenal in which the most clinically relevant toxin families expanded through posterior gene duplications, regardless of whether they originated by duplication or gene co-option.

Keywords: co-option; gene recruitment; genome; snake venom; toxin evolution.

Fig. 1.

Schematic diagram of the genomic sequencing strategies used to obtain toxin genes and their flanking regions in B. jararaca.

Fig. 2.

Schematic architecture of venom gene loci of different toxins showing syntenic blocks among different species. Red pentagon, toxin gene; yellow pentagon, nontoxin ortholog of a toxin gene (or paralog if in the same species); and white pentagon, flanking nontoxin gene. In each box, the name of the ortholog representing the putative ancestral gene for the toxin family is noted bellow the toxin family name, followed in parenthesis by the gene ID of a reference gene (which is noted in blue and outlined in blue in the scheme) from an organism that do not contain the toxic character for this family. Gene names are indicated over the array of orthologs or within pentagons. Some paralogous genes are represented by one pentagon internally marked with the number of paralogs occurring in the species. Relevant pseudogenes are indicated with Ψ. Species were classified according to the following color code: green box, venomous snake; blue box, nonvenomous snake; orange box, nonsnake Squamata; and gray box, none of the above. Species codes and GenBank Genome ID or segment accession number are as follows: Bo.jara, Bothrops jararaca (this study); An.caro, Anolis carolinensis, ID: 708; Cr.adam, Crotalus adamanteus, PLA2 scaffold KX211996; Cr.atro, Crotalus atrox, PLA2 scaffold KX211994; Cr.scut, Crotalus scutulatus, ADAM28 scaffold MT032003.1; Cr.viri, Crotalus viridis, ID: 71654; Ga.gall, Gallus gallus, ID: 111; Ge.japo, Gekko japonicus, ID: 40475; La.chal, Latimeria chalumnae, ID: 3262; Mu.musc, Mus musculus, ID: 52; Op.hann, Ophiophagus hannah, ID: 10842; Po.mura, Podarcis muralis, ID: 8765; Po.vitt, Pogona vitticeps, ID: 7589; Pr.mucr, Protobothrops mucrosquamatus, ID: 18192; Ps.text, Pseudonaja textilis, ID: 72610; and Py.bivi, Python bivittatus, ID: 17893. For B. jararaca, the scheme is based on the combination of data gathered from the sequences obtained by the different strategies used in this work. CTLs and LAAO were not included in the figure since the B. jararaca scaffolds did not provide enough information to define the architecture of their loci.

Fig. 3.

The SVMP gene structure and arrangement at the locus. (A) Architecture of the ADAM28 genomic locus in different vertebrates (not in scale). The putative SVMP ancestral gene ADAM28 (yellow arrow) and flanking genes (STC1, NEFM, and NEFL: white arrows) form a syntenic block among vertebrates. Orange and beige arrows represent ADAM family genes in humans (ADAMDEC1 and ADAM7). Red and pink arrows represent genes from the SVMP classes P-III and P-II, respectively. Solid lines are contiguous sequences, and dotted lines indicate uncertain order or no contiguity. Blue bars represent regions covered by BAC. (B) Schematic alignment of SVMP gene structures showing the conservation of exons (squares) between SVMP P-III and P-II. A short and a long deletion at exon 14 of SVMP P-II are marked. These deletions result in the loss of the Cys-rich domain and the shortening of the disintegrin-like sequence through the acquisition of a new stop codon (red star) preceding the original one (black star). (C) Details of the nucleotide alignment with the encoded amino acid residues between the two neighboring SVMP genes belonging to the P-III and P-II classes (BJARBC_SVMP3_g20 and BJARBC_SVMP2_g07, respectively) in the region between exons 14 and 17.

Fig. 4.

BPP/CNP gene structure. (A) Schematic alignment of BPP/CNP and CNP genes from different organisms emphasizing the correspondence of introns and exons, the conservation of domain structures, and the extension of exon 1 harboring the BPPs in Viperidae. Species are classified according to the following color code: green box, venomous snake; orange box, nonsnake Squamata; and gray box: none of the above. (B) Part of the BPP/CNP gene sequence from B. jararaca and its translation, showing that BPPs are restricted to exon 1.

Fig. 5.

(A) Architecture of the venom VEGF-F gene and nonvenom VEGF-A gene loci in synteny among different organisms. Red pentagon, VEGF-F toxin gene; yellow pentagon, PGF-like or VEGF-F–like genes; green pentagon, VEGF-A gene; and the white pentagon represents adjacent nonrelated genes. Dots at the end of solid lines indicate scaffold ends. Blue bar represents region covered by BAC. Gene representations are not to scale. Species were classified according to the following color code: green box, venomous snake; blue box, nonvenomous snake; orange box, nonsnake Squamata; and gray box, none of the above. (B) Summarized phylogenetic tree of the VEGF family of growth factor focusing the origin of VEGF-F (snake venom VEGFs) from the PGF-like/VEGF-F–like ortholog. The complete phylogenetic analysis is shown in SI Appendix, Fig. S4. (C) Schematic comparison of VEGF-F and VEGF-A genes in three Squamata, pointing out the levels of conservation throughout these genes. Percentage values on the right represent pairwise identity of CDS regions.