Sequence Divergence in Venom Genes Within and Between Montane Pitviper (Viperidae: Crotalinae: Cerrophidion) Species is Driven by Mutation-Drift Equilibrium

Abstract

Snake venom can vary both among and within species. While some groups of New World pitvipers-such as rattlesnakes-have been well studied, very little is known about the venom of montane pitvipers (Cerrophidion) found across the Mesoamerican highlands. Compared to most well-studied rattlesnakes, which are widely distributed, the isolated montane populations of Cerrophidion may facilitate unique evolutionary trajectories and venom differentiation. Here, we describe the venom gland transcriptomes for populations of C. petlalcalensis, C. tzotzilorum, and C. godmani from Mexico, and a single individual of C. sasai from Costa Rica. We explore gene expression variation in Cerrophidion and sequence evolution of toxins within C. godmani specifically. Cerrophidion venom gland transcriptomes are composed primarily of snake venom metalloproteinases, phospholipase A[Formula: see text]s (PLA[Formula: see text]s), and snake venom serine proteases. Cerrophidion petlalcalensis shows little intraspecific variation; however, C. godmani and C. tzotzilorum differ significantly between geographically isolated populations. Interestingly, intraspecific variation was mostly attributed to expression variation as we did not detect signals of selection within C. godmani toxins. Additionally, we found PLA[Formula: see text]-like myotoxins in all species except C. petlalcalensis, and crotoxin-like PLA[Formula: see text]s in the southern population of C. godmani. Our results demonstrate significant intraspecific venom variation within C. godmani and C. tzotzilorum. The toxins of C. godmani show little evidence of directional selection where variation in toxin sequence is consistent with evolution under a model of mutation-drift equilibrium. Cerrophidion godmani individuals from the southern population may exhibit neurotoxic venom activity given the presence of crotoxin-like PLA[Formula: see text]s; however, further research is required to confirm this hypothesis.

Keywords: Gene family evolution; Mutation–drift equilibrium; Selection; Transcriptomics.

Fig. 1

Cerrophidion distribution in Mesoamerica. Map modified from VenomMaps (Rautsaw et al. 2022), with the localities of samples of each species used herein. Species are represented by different shapes, different outline colors correspond to Northern populations of that species.*Cerrophidion wilsoni was not included in this work. Species tree scaled with IQtree from the inferred Astral tree; support values correspond to the Astral species tree. Node shapes correspond with the populations in the map. Pie charts in the tips show the percentage of expression of the five most abundant toxin families; the remaining toxin families are included in the category “Other.” Photo Credit: Jason M. Jones (C. tzotzilorum)

Fig. 2

RSEM results for the consensus transcriptomes of A C. sasai, B. C. petlalcalensis, and C. C. tzotzilorum. In A (I) barplot of the log ranked expression of toxin genes, (II) pie charts of the percent expression of each toxin family average of all individuals. In B and C, (I) barplot of the log ranked expression of toxin genes, (II) stacked barplots with the percent expression of each toxin family by sampled individual, (III) pie charts of the percent expression of each toxin family average of all individuals. Photo Credit: A R. Wayne VanDevender (C. sasai), B Carlos E. Montaño-Ruvalcaba (C. petlalcalensis), and C Ramses A. Rosales-García (C. tzotzilorum)

Fig. 3

RSEM results for the A average transcriptome of all individual of C. godmani; B average of the northern population; and C average of the southern population. In A (I) barplot of the log ranked expression of toxin genes, (II) stacked barplots with the percent expression of each toxin family by sampled individual, (III) pie chart of the percent expression of each toxin family for individual populations and for all the individuals. In B and C, (I) barplot of the log ranked expression of toxin genes; (II) pie chart of the percent expression of each toxin family for individual populations and for all the individuals. Photo Credit: Carlos E. Montaño-Ruvalcaba (C. godmani)

Fig. 4

A Consensus maximum likelihood tree of the PLA₂s in Cerrophidion including PLA₂s used in Whittington et al. (2018), Mason et al. (2020), Neri-Castro et al. (2020b), and from Genbank (accession numbers in online resource 1, Table S1). The Cerrophidion PLA₂s (names highlighted in blue) are numbered by the toxin’s average expression for each species. Acidic and basic PLA₂s are identified by red and blue branches, respectively, based on the hypothetical isoelectric point of the amino acid sequences. Nodes with a black dot have >75 bootstrap support. Cerrophidion crotoxin subunit homologs are identified by a star (*) next to the name. B Species tree scaled with IQtree from the inferred Astral tree; lineages with crotoxin-like subunit homologs are purple; support values correspond to the Astral species tree. C Amino acid alignment of the gA2 clade and the hypothetical homolog from Cerrophidion, dots represent no change from the reference sequence (C._godmani_11). Sites represented with bars match cleavage sites identified in Whittington et al. (2018); a black star (*) at site 5 is the key substitution known in Bothriechis, Crotalus, and Gloydius; a red star (*) is at the alternative cleavage site in C. godmani based on the protein cutter tool from ExPASy server (http://web.expasy.org/peptide_cutter/; Gasteiger et al. 2005) (Color figure online)

Fig. 5

Heatmap showing the log TPM expression of toxins identified as differentially expressed in C. godmani ordered by the average expression. In the left columns (Pop & SVL) the darker colors represent significant differential expression agreement by both DESeq2 and edgeR (FDR< 0.05)

Fig. 6

Selection plots. Top: estimates of selection using A Tajima’s D, B F_ST, and C Likelihood Ratio Test (LRT) for the BUSTED model for Toxins and Nontoxins, each with the Nontoxin 95th percentile (dotted lines) to identify outlier toxins. The toxin family and the rank based on highest-to-lowest average expression in the transcriptome is displayed for toxins which fall outside the 95th percentile. Bottom: Linear regressions of the Toxin’s mean expression (Average TPM) and estimates of selection including D Tajima’s D, E F_ST, and F LRT of the BUSTED model. For Tajima’s D, dotted lines are regressions considering all the transcripts (center), just positive values (top) and just negative values (bottom)