Population Genomic Analysis of a Pitviper Reveals Microevolutionary Forces Underlying Venom Chemistry.

Venoms are among the most biologically active secretions known, and are commonly believed to evolve under extreme positive selection. Many venom gene families, however, have undergone duplication, and are often deployed in doses vastly exceeding the LD50 for most prey species, which should reduce the strength of positive selection. Here, we contrast these selective regimes using snake venoms, which consist of rapidly evolving protein formulations. Though decades of extensive studies have found that snake venom proteins are subject to strong positive selection, the greater action of drift has been hypothesized, but never tested. Using a combination of de novo genome sequencing, population genomics, transcriptomics, and proteomics, we compare the two modes of evolution in the pitviper, Protobothrops mucrosquamatus. By partitioning selective constraints and adaptive evolution in a McDonald-Kreitman-type framework, we find support for both hypotheses: venom proteins indeed experience both stronger positive selection, and lower selective constraint than other genes in the genome. Furthermore, the strength of selection may be modulated by expression level, with more abundant proteins experiencing weaker selective constraint, leading to the accumulation of more deleterious mutations. These findings show that snake venoms evolve by a combination of adaptive and neutral mechanisms, both of which explain their extraordinarily high rates of molecular evolution. In addition to positive selection, which optimizes efficacy of the venom in the short term, relaxed selective constraints for deleterious mutations can lead to more rapid turnover of individual proteins, and potentially to exploration of a larger venom phenotypic space.

Introduction

Understanding how natural selection leads to adaptation requires a mechanistic explanation of how alternative genotypes lead to fitness differences (Storz and Wheat 2010). Venoms, particularly proteinaceous venoms of snakes, provide an excellent model system for the study of adaptation, because venom secretions represent quantifiable phenotypes, which are under control of easily characterizable genomic loci. Venom proteins often evolve rapidly (Casewell etal. 2013; Aird etal. 2015), and phenotypic studies have shown that they are more efficacious against prey types typically encountered and consumed (Mackessy 1996; Jorge da Silva and Aird 2001; Gibbs and Mackessy 2009). The study of snake venom protein evolution, like that of many adaptive traits, has emphasized the role of positive selection on protein structure (Sunagar etal. 2014, 2015). The role of drift for protein evolution has not been explicitly considered, except as a hypothesis to be rejected, a view that may fail to capture the complexity of selective regimes.

The strike of a snake initiates complex biochemical cascades within the prey organism that have immediate life and death consequences (Aird 2002). The venom cocktail must act quickly to immobilize the prey before it escapes or has a chance to injure the snake. Protein chemistry mediates the outcome of envenomation, and venom genes show extraordinarily high rates of evolution (Nakashima etal. 1995; Deshimaru etal. 1996). Yet, despite decades of interest in this phenomenon, still no consensus exists about the dominant forces driving venom protein sequence changes. On one hand, a large body of molecular evidence argues that venom components are subject to strong positive selection (i.e., “the adaptive hypothesis;” Nakashima etal. 1995; Kordis and Gubensek 2000; Zupunski etal. 2003; Zhu etal. 2004; Lynch 2007; St Pierre etal. 2008; Barghi etal. 2015; Pedroso etal. 2015). On the other hand, the complexity of venoms entails pharmacological redundancy, with more venom components injected than minimally necessary to kill or otherwise immobilize the prey (i.e., the “overkill hypothesis”), predicting relaxed functional constraints upon individual venom components (Mebs 2001).

The overkill hypothesis relies on the empirical observation that some individual snakes can carry enough venom to kill up to thousands of lab mice (Broad etal. 1979). However, it has generally been dismissed as an oversimplified model of the envenomation process (Casewell etal. 2013), and in the face of strong molecular data supporting positive selection on venom components (Blumstein 2006; Barlow etal. 2009). In addition, the overkill hypothesis may not hold for a number of species of specialized snake prey or predators that have evolved resistance to venom (Perez etal. 1978; Poran etal. 1987; Heatwole and Poran 1995; Weissenberg etal. 1997; Jansa and Voss 2011; Drabeck etal. 2015; Holding, Biardi, etal. 2016; Holding, Drabeck, etal. 2016). However, the overkill hypothesis has never been empirically tested, since merely detecting positive selection on a gene does not preclude its also being subject to more relaxed selection. In addition to the overkill hypothesis, there could be other reasons for high rates of drift in venom genes following from birth-and-death evolutionary models (Nei and Rooney 2005). Consequently, most discussions of venom structure focus on adaptive explanations, whereas the extent to which neutral forces actually shape it remains unknown.

Although both the adaptive and overkill hypotheses make clear-cut predictions, it has been difficult to test them comprehensively in the same analytical framework until now. Previous studies have typically employed ratios of nonsynonymous and synonymous substitutions (dn/ds) to study selection. This ratio results from the joint action of both neutral (ds) and selective forces (dn), making it difficult to partition and to examine them separately. The advent of population genomics introduced a wide array of other proxies for detecting selection acting on genes and genomic regions (Stinchcombe and Hoekstra 2008; Ekblom and Galindo 2011). However, few approaches quantify selection and drift within the same framework. Exceptionally, McDonald–Kreitman-type tests use population-level data to examine rates of fixation of synonymous and nonsynonymous mutations between two diverging lineages, allowing additional parameters to be estimated (McDonald and Kreitman 1991; Smith and Eyre-Walker 2002). Methods based on the MK test are the only ones that can, in principle, yield an unbiased quantitative assessment of the contribution of positive selection (Eyre-Walker 2006) to venom evolution.

In addition, this framework can simultaneously infer the proportion of substitutions driven by positive selection (α) and the selective constraint for each gene, that is, the rate at which deleterious mutations accumulate. These values are of key importance for assessing the merits of the adaptive and overkill hypotheses, which make contrasting predictions. Namely, under adaptive evolution, venom genes should accumulate adaptive substitutions at greater rates, compared with nonvenom genes. By contrast, if venom genes are tolerant of deleterious mutations, they would experience weaker selective constraint. In addition, the McDonald–Kreitman-type framework can also estimate other parameters, such as mutation rates (Θ = 4Neμ) for each gene. Higher mutation rates have been proposed as an alternative explanation for the more rapid evolutionary rates of venom proteins, but this hypothesis has never been comprehensively tested in a genome-wide analysis (Kini and Chan 1999; Kini and Chinnasamy 2010).

In the present study, we examined the relative contributions of adaptive selection and genetic drift to evolution of venom and housekeeping genes, the latter defined as genes expressed in the venom gland, but not secreted into the venom proteome. We focused on the elaborated venom system of a pitviper, the Taiwan habu (Protobothrops mucrosquamatus), comparing it with its sister species, the Sakishima habu (Protobothrops elegans), from which it split approximately 2.9 Ma (Guo etal. 2007; Hedges etal. 2015). We find that both positive selection and genetic drift play a role in the rapid sequence changes seen in venom genes. However, expression levels of individual venom proteins may modulate this relationship, with the most abundant proteins experiencing stronger drift.

Materials and Methods

Sample Collection

Thirty specimens of Protobothrops mucrosquamatus were collected at various localities throughout Okinawa, Japan: in Nago City (28) and Nakijin Village (1) in the Motobu Peninsula and in Onna Village (1), south of Nago. Venom was extracted from all specimens at Day 0 and venom glands were harvested at Days 1, 2, 4, and 8. Internal organs were also removed and frozen in liquid nitrogen, after which they were maintained at −80 °C until use. Venom gland transcriptomes were created for all specimens and a reference transcriptome was created from eight specimens, with two specimens selected randomly from each of the gland harvest days (Day 1: Pm_2, Pm_5; Day 2: Pm_8, Pm_11; Day 4: Pm_30, Pm_17; Day 8: Pm_23, Pm_27) (supplementary table S1, Supplementary Material online). The experimental protocols used for handling animals in this study have been approved by the Animal Resources Section of OIST.

Genome Assembly

The reference genome was assembled from liver tissue from sample Pm_15. We prepared PCR-free Illumina TruSeq libraries size-selected to 450 ± 25 bp, which were sequenced on a HiSeq 2500 instrument to a mean depth of 69×. Scaffolds were then assembled using DISCOVAR (Weisenfeld etal. 2014). We also prepared 5 and 15 kb mate-pair libraries using Illumina’s Nextera protocol and sequenced them on an Illumina HiSeq 2000 to a depth of approximately 100× each. Reads from mate-paired libraries were trimmed and filtered using NextClip (Leggett etal. 2014) and used to scaffold contigs using SSPACE (Boetzer etal. 2011).

Transcriptomic and Proteomic Analyses

Venom gland RNA seq libraries were prepared as described previously (Aird etal. 2013), except that ERCC92 synthetic spike-ins were added to the RNA extracts as described by Aird et al (2015) for quality control. In addition, pooled libraries were normalized using the Evrogen Trimmer-2 cDNA normalization kit and sequenced to improve genome annotation for nonvenom transcripts. Reads were mapped to predicted coding sequences using Bowtie2 (Langmead and Salzberg 2012) within the RSEM package (Li and Dewey 2011). Though expression data doesn’t accurately always reflect protein levels in the venom (Casewell etal. 2014), it does in Protobothrops (Aird etal. 2013, 2015). Therefore, we used the number of fragments per kilobase mapped (FPKM) as a measure of protein abundance. Mapped read counts for each sample, as estimated by RSEM, can be found in supplementary table S1, Supplementary Material online. Proteomic analysis was conducted as described previously (Aird etal. 2013, 2015).

Variant Calling

Variants in the coding sequences were called using four separate algorithms: GATK, freebayes, platypus, and samtools (Li etal. 2009; McKenna etal. 2010; Garrison and Marth 2012; Rimmer etal. 2014). Variants were converted to allelic primitives using GATK, and a consensus set was called by BAYSIC using a 0.80 posterior probability (Cantarel etal. 2014). Finally, indels, variants with >10% missing data in the ingroup samples, and variants with more than two alleles were filtered out. We then classified variants as synonymous or nonsynonymous using SnpEff (Cingolani etal. 2012), and noted which were fixed versus polymorphic with respect to the outgroup.

Resequencing

Twenty P. mucrosquamatus samples used for RNA-seq and one P. elegans were selected for whole-genome resequencing. Library construction and sequencing was performed by BGI (Shenzhen, China) on a HiSeq 2000 instrument in PE150 mode, using Illumina TruSeq libraries. BGI then trimmed sequencing adaptors using in-house scripts. Raw reads were then mapped to the reference assembly using NextGenMap (Sedlazeck etal. 2013). Sample coverages computed by GATK can be found in supplementary table S1, Supplementary Material online. Copy number variation was assessed with CNV-seq using default parameters (Xie and Tammi 2009). A gene was considered to have copy number variation if a CNV-seq window with a significant signature of copy number variation intersected an exon of the gene model.

Population Genomics

A table of fixed or polymorphic replacement and silent substitutions was used to infer population genomic parameters using two approaches: a Bayesian approach (SnIPRE) (Eilertson etal. 2012), and a maximum likelihood approach (MKTest) (Welch 2006). The two packages compute parallel measures of selection for every gene. Both of them model mutational constraint, although it is defined differently. In SnIPRE it is a regression model coefficient, with more negative values indicating that polymorphic mutations are either being fixed or eliminated at a higher rate than synonymous mutations, whereas the MKtest estimates a fraction where 1-f of the mutants are under strong purifying selection. In either case, larger estimated values indicate relaxed selective pressure. The overall results were qualitatively similar for the two estimates of selective constraint, but gave radically different results for serine proteases, with SnIPRE estimating relatively relaxed constraint for this gene family, and MKTest estimating stronger selective constraint than experienced by the average gene in the genome. Almost all of the MKTest’s estimates of f lay on boundary values (0, 1), suggesting a relatively poor fit to the model. We include them for completeness, but are inclined to trust them less than those of SnIPRE.

Both packages also estimate mutation rates as Θ (= Ne*μ). In addition, MKTest estimates the proportion of nonsynonymous divergence driven by positive natural selection, a measure of true adaptive evolution (α) (Smith and Eyre-Walker 2002). In addition, both approaches model divergence times (SnIPRE: a coefficient called τ; MKTest: expected neutral divergence per site: μ*t) and, and weakly selected beneficial mutations resulting from mildly deleterious substitutions (γ = 4Nes). These are not presented here because they lack testable predictions and straightforward biological interpretations in this context.

We then used the R statistical package (R core team 2015) to conduct a range of statistical analyses on the population genomic coefficient estimates. Given the relatively small number of venom genes, and the frequently nonnormal distribution of the coefficient estimates, we used nonparametric statistics, namely Spearman’s rank correlation, and Kruskal–Wallis test adjusted for multiple comparisons using the method of Siegel and Castellan (Siegel and Castellan 1988). When comparing the fraction of adaptive substitutions acting on gene classes, we eliminated any genes with negative values of α, which is defined only for the range between 0 and 1, and negative values are produced by sampling error or violations of the model (Eyre-Walker 2002).

Mapping Evolutionary Changes in Structure of Phospholipases A2 and Serine Proteases

For each analysis, sequences of venom components from the reference genome were codon-aligned using MAFFT (Katoh etal. 2002), and sites under selection were analyzed using the REL algorithm implemented in the HyPhy package on the DataMonkey server (Kosakovsky Pond and Frost 2005; Pond and Frost 2005).

A consensus sequence of amino acids was formed from the codon alignment and used as input for developing a PLA2 model in SWISS-MODEL (Arnold etal. 2006; Kiefer etal. 2009). The suitable model was selected and edited using PyMOL (Delano 2002). Sites on the protein structure were color-coded based on the REL results, where a Bayes Factor of 50 or more for dN > dS, representing negative selection, was colored blue, whereas a Bayes Factor of 50 or more for dN < dS, representing positive selection, was colored red.

Data Accessibility

Genome annotation is available on NCBI (Protobothrops mucrosquamatus Annotation Release 100). Raw genomic reads and RNA-seq are available under NCBI accessions PRJNA313429 and PRJDB4386, respectively. Snakemake files (Köster and Rahmann 2012) for genome assembly and variant calling, as well as a script for R-based analysis and intermediate data are in a MySQL database available at https://github.com/mikheyev/mucrosquamatus-selection, last accessed September 29, 2017. A virtual machine necessary to reproduce the R-based analyses is provided at http://mybinder.org/repo/mikheyev/mucrosquamatus-selection, last accessed September 29, 2017. Note: the virtual machine relies on a third-party service, and its availability is not guaranteed. Therefore, we include a Dockerfile necessary to build such a virtual machine locally using the Docker platform. Finally, a human-readable version of the analysis leading to the principal results and generating draft versions of all the figures is available at https://mikheyev.github.io/mucrosquamatus-selection/, last accessed September 29, 2017.

Results

We produced the first genome assembly of a pitviper, P. mucrosquamatus. It was 1.63 gigabases long (contig and scaffold N50, 22 and 424 kb, respectively). The genome was annotated using RNA-seq data from 32 individual snakes, producing an annotation with 20,122 protein-coding genes. We also resequenced 20 additional individual P. mucrosquamatus genomes and a P. elegans genome as an outgroup (supplementary table S1, Supplementary Material online). These data were used to calculate the number of fixed and replacement synonymous and nonsynonymous substitutions between the two species, which were used to compute estimates of selective constraint and positive selection for every polymorphic gene.

Venom Composition

We detected 76 proteins in the venom of P. mucrosquamatus using mass spectrometry. However, these venoms were dominated by a few phospholipase A2, serine protease, and metalloprotease genes (fig. 1), with the rest of the proteins, including other members of these three families, playing minor roles. Both snakes had strongly correlated patterns of gene expression, suggesting conservation of their regulatory machinery (rs = 0.73, N = 55, P < 2*10−16, supplementary, fig. S1, Supplementary Material online). Thus, both species should experience similar effects of gene expression on selection. To the best of our knowledge, there have been no formal dietary studies of either species; however, anecdotal reports and our own field observations suggest that both are dietary generalists, preying upon frogs, lizards, birds, rodents, and shrews. That being said, percentages of each prey type may vary locally with availability, making it is impossible to draw any firm conclusions about how venom composition in either species may be related to diet.

. 1.

—A small number of genes, belonging to just three families, dominates the venom of P. mucrosquamatus. Transcript proportions are correlated with proteomic abundance (Aird etal. 2013, 2015). Abbreviations: 5′NT 5′-nucleotidase, AChE acetylcholinesterase, CRISP cysteine-rich secretory proteins, CTL C-type lectin-like proteins, HYAL hyaluronidase, LAO L-amino acid oxidase, MP metalloprotease, NATR C-type natriuretic peptide, NGF nerve growth factor, PDE phosphodiesterase, PLA2 phospholipase A2, PLB Phospholipase B, PLCl phospholipase A2 inhibitor, QC glutaminyl cyclase, SP serine protease, VEGF vascular endothelial growth factor.

Structural Organization of Venom Gene Clusters

Protobothrops mucrosquamatus venom includes four highly diversified gene families: C-type lectin-like proteins (CTLs), metalloproteases, serine proteases, and phospholipases A2. The latter two families are located in tandem arrays of duplicated genes, each on its own scaffold (NW_015386730.1 and NW_015387341.1). Whether CTLs and metalloproteases have similar organization is not certain, given the fragmented nature of the assembly, but half of the genes in each family reside on just a few scaffolds. Previous work has shown that phospholipases A2 in the Okinawa habu (Protobothrops flavoviridis) also occur in a single gene cluster (Ikeda etal. 2010), though the copy number and orientation of genes differs from that reconstructed by the present study for P. mucrosquamatus. Two phospholipase A2 genes found in P. mucrosquamatus (LOC107291353 and LOC107291356) appear to be missing from the closely related P. elegans genome. There also appears to be copy number variation within many P. mucrosquamatus venom genes (supplementary fig. S2, Supplementary Material online). These findings suggest that venom gene clusters may be structurally polymorphic, with genes being duplicated and lost both within a species and between closely related species.

Do Venom Genes Mutate More Rapidly?

Differences in evolutionary rates can also potentially result from unequal rates of background mutation (Kini and Chinnasamy 2010). However, we found no difference in mutation rates (Θ) between venom and nonvenom genes (SnIPRE P = 0.059, MKTest P = 0.73), suggesting that drift and selection are the major forces behind venom protein sequence change.

Effects of Gene Expression on Venom Protein Evolution

There was no effect of average expression level of venom proteins on the proportion of substitutions driven by positive selection (rs = 0.017, P = 0.92). However, there was evidence that selective constraints on the most abundant toxins are diminished (SnIPRE rs = 0.016, P = 0.016; MKTest rs = 0.28, P = 0.075; fig. 3). Relaxed constraints lead to accumulation of deleterious mutations, and an observed increase in the dn/ds ratio.

. 3.

—More abundant venom components experience increasingly relaxed selective constraint. As in figure 2, higher values of the selective constraint coefficient, indicate less effective elimination of deleterious mutations. MKTest results were consistent with those of SnIPRE, but marginally nonsignificant (P = 0.075). Because estimates of mutational constraint differed widely between SnIPRE and MKTest, they are excluded from this analysis, though including them does not qualitatively change results. Abbreviations as in figure 1. These findings suggest that selective constraint can be modulated by protein expression level.

Discussion

Selective Forces Acting on Venom

Adaptive Evolution

The role of adaptive selection in snake venom evolution is uncontroversial, and we also find that venom genes have a greater proportion of substitutions driven by positive selection (α), than the rest of the genome (fig. 2). Overall, there was weaker support for selection acting on individual venom protein classes (supplementary fig. S2, Supplementary Material online), largely because per-class sample sizes are significantly smaller. Sequence-based analysis of the dominant venom components, the phospholipases A2 and serine proteases, identified multiple sites under positive selection, particularly on the exteriors of these proteins, which suggest accelerated evolution of functional components (Kini and Chan 1999) (supplementary fig. S3, Supplementary Material online). Strong adaptive selection is consistent with a wide suite of ecological observations, further supporting the adaptive value of snake venoms. For example, snake venoms tend to be more effective for subduing their preferred prey (Silva Jr and Aird 2001; Gibbs and Mackessy 2009) and even closely related taxa with divergent diets employ strongly divergent venoms (Sanz etal. 2006). Perhaps most dramatically, the adoption of a fish egg-only diet by the sea snake, Aipysurus eydouxii, has led to relaxed selection on its venom, and a variety of dysfunctional mutations in its phospholipase A2 genes (Li etal. 2005). It is therefore likely that positive selection is the dominant driver of snake venom evolution.

. 2.

—Venom proteins show higher rates of adaptive evolution (A), and relaxed selective constraint (B, C) compared with housekeeping genes. Violin plots summarize the distributions of the parameter estimates, with means given by red dots. Statistical significance is computed using Kruskal–Wallis tests. The proportion of adaptive substitutions fixed by natural selection refers to those in relatively recent evolutionary time, since the divergence of the P. mucrosquamatus and P. elegans, without distinguishing between the possible modes of natural selection (e.g., episodic vs. continuous). Selective constraint was estimated with two software packages that use different values for their coefficients, but in both cases higher values mean lower selective constraint. These data sets support both adaptive and overkill hypotheses, respectively.

The Overkill Hypothesis and the Role of Drift in Venom Proteome Evolution

We likewise find that venom components experience less selective constraint than housekeeping genes (fig. 2 and supplementary fig. S4, Supplementary Material online). Although this finding may, at first, seem to contradict the higher levels of adaptive evolution, it is important to remember that selective pressures may vary across the secondary and tertiary structure of a given protein (supplementary fig. S5, Supplementary Material online), and that, in principle, some parts of the protein can experience positive selection, whereas other parts are free to evolve neutrally. Consequently, relaxed selective constraints do not imply the absence of positive selection, but rather its relative strength. Indeed, numerous studies, and our data from P. mucrosquamatus, show that major venom components in snakes are under strong selection (Kordis etal. 1998; Kordis and Gubensek 2000; Zupunski etal. 2003; Juárez etal. 2008; Barlow etal. 2009; Casewell etal. 2011; Vonk etal. 2013; Malhotra etal. 2015). That being said, reduced selective constraint on venom proteins supports predictions of the overkill hypothesis, and in turn, has three interesting and underexplored consequences for venom proteins: 1) increased rates of pseudogenization, 2) higher levels of standing variation, and 3) venom evolution may happen neutrally during vicariant speciation.

First, accumulation of deleterious mutations could lead to pseudogenization. Indeed, there is evidence that this has occurred in the phospholipase A2 cluster of P. flavoviridis, which contains two pseudogenes and three active genes (Ikeda etal. 2010). Furthermore, allopatrically isolated populations of this species differ in the presence or absence of pseudogenes (Chijiwa etal. 2000; Ikeda etal. 2010). Although the authors attribute such interisland variability to selection, based on our findings, it also seems possible that at least some of the variation is due to neutral mutation and drift, which would be particularly strong in small isolated populations. Furthermore, as pseudogenized genes no longer produce functional proteins, their expression should be strongly selected against, ultimately eliminating them from the venom, whereas other proteins take their place. Consequently, drift, possibly followed by negative selection on pseudogenized genes, can lead to accelerated turnover of individual venom components. Indeed, gene loss was found to be a key factor in the evolution of phospholipase A2 genes in rattlesnakes (Dowell etal. 2016).

Second, higher levels of variability associated with relaxed selective constraint may expand the phenotypic space available to snake venoms. The concept of phenotypic space is analogous to Hutchinson’s concept of the niche, as an N-dimensional hypervolume (Hutchinson 1957). In venomic phenotypic space, each venom constituent multiplied by its concentration could be conceptualized as a dimension in N-dimensional phenotypic space, such that the incorporation or deletion of venom constituents, or changes in their abundance, would add, delete, lengthen, or shorten dimensions. It seems highly likely that for any given snake species, there could be multiple biochemical strategies that might be equally efficacious, particularly if it feeds on multiple prey types. Nonetheless, members of a venomous snake species, in essence, represent elements of a living array that performs what a pharmaceutical company might call “pharmaceutical lead optimization.” That is, many toxin variants and their combinations are screened continuously and simultaneously across populations.

Although positive selection promotes local optima, genetic drift may allow populations to cross fitness valleys, and potentially leads to new venom compositions (Wright 1932). For example, previously dominant toxins could accumulate deleterious mutations and become less abundant, whereas more minor components could acquire beneficial mutations and increase in concentration. This could result from either changes in their regulatory sequences, or from changes in copy number (Aird etal. 1991; Malhotra 2015). Although such evolutionary dynamics would be difficult to conclusively demonstrate, relative abundance of venom components can indeed change relatively rapidly. For instance, a phospholipase A2 comprising 73.4% of the P. elegans transcriptome, is barely present in venom of its closest relatives, P. flavoviridis (<0.05%) (Aird etal. 2015) and also in P. mucrosquamatus (present study, supplementary fig. S1, Supplementary Material online). Even within a species, there can be significant geographic variation in protein abundance (Mebs and Kornalik 1984; Chijiwa etal. 2003), as demonstrated by studies of venom composition among littermates (Chippaux etal. 1982; Pintor etal. 2011). Future work focused on phylogenetic reconstruction of ancestral states of proteins, and their expression levels, should provide important insights into the interplay between expression, selection, and venom evolution, as has been done for a range of other animals (Morandin etal. 2016).

Third, the observation that venoms are more susceptible to the effects of genetic drift than the rest of the genome, has important implications for venom differences between species. Vicariance effects, and in particular founder effects, will change effective population sizes and further exacerbate the effects of drift on venom components, possibly leading to significant amounts of differentiation between populations or nascent species. Such regional variation is found in Okinawan Protobothrops, which differs in venom composition between islands (Chijiwa etal. 2003). Many of these islands are small and isolated, likely harboring small populations of snakes. Even more recent work on the eastern diamondback rattlesnake (Crotalus adamanteus) has shown that populations vary greatly in copy number at toxin-encoding loci, and that selection for increased expression, rather than sequence diversity has driven duplications (Margres etal. 2017). Therefore, neutral divergence should be considered when comparing possible adaptive differences between species.

Effect of Expression on Venom Evolution

Such an increase has, in fact, recently been documented by us (Aird etal. 2015). At that time, we interpreted the higher dn/ds rate as evidence of stronger adaptive selection on the most abundant proteins, whereas the present analysis suggests that the opposite is true, highlighting the advantages of the McDonald–Kreitman framework over traditional dn/ds statistics.

At first glance, the relatively lower strength of selection on more abundant venom components seems counterintuitive. However, we propose a model to explain this pattern, based on several additional observations. First, in order to assure rapid immobilization of prey, snakes inject an optimal amount of venom, which is several times more than the minimally necessary dose (Allon and Kochva 1974; Hayes 1995; Herbert and Hayes 2008). Second, studies of P. flavoviridis venom found that the venom fraction containing hemorrhagic metalloproteases, second only to phospholipases in abundance (Aird etal. 2013), showed roughly twice the lethality of crude venom (Ohsaka 1960; Takahashi and Ohsaka 1970). However, Aird (2002) suggested that immobilization via hypotension and paralysis, rather than lethality, is the objective of envenomation, and by extension, the trait upon which selection acts. In envenomated mice, immobilization time is roughly 3–5× shorter than the time to death (Herbert and Hayes 2008). Therefore, if excess venom is injected to rapidly immobilize prey, concentrations of the most abundant venom components may indeed reach ‘overkill’ levels. Thus, an excess of the most abundant components could lead to relaxed selective constraints upon them, since deleterious mutations could be compensated by higher concentration.

However, can higher concentrations rescue mutated proteins? Indirect insights into this question come from site-directed mutagenesis studies that look at the activity of venom components. Because such studies typically target sites of likely biochemical importance, and make nonrandom alterations of the protein, they explore the more extreme end of the protein fitness landscape. However, site-directed mutagenesis studies generally find that, except for mutations in the active site, which tend to be strongly deleterious, most other mutations have either no effect, or relatively mild effects on enzymatic properties of the venom (e.g., Trémeau etal. 1995; Dennis etal. 1993; Zhang etal. 1997). Such decreases in catalytic or noncatalytic activity can, in principle, be rescued by increased protein concentration. Conversely, if the levels of the protein were high in the first place, mutations may be masked as long as the decrease in activity is relatively modest.

Analytical Caveats

Whole-genome inference of population-level parameters is difficult, and subject to a number of important assumptions. First, changes in demography may affect parameter estimates. Although this could affect absolute values of the parameters, our major conclusions should not be affected, since we are making comparisons within a genome. Another important assumption is the independence of genes. Since venom genes are found in large clusters, they may be effectively linked, possibly biasing results. However, the different venom classes are most likely unlinked, and the omnibus tests presented in figure 2, should be relatively robust. Unfortunately, linkage at the family level affects all previous studies of selection acting on venom genes as well. Finally, copy number variation may affect some parameter estimates by possibly inflating the extent of polymorphism existing in some genes. Filtering out any genes that display evidence of copy number variation still gives a higher overall level of adaptive substitution in venom genes (P = 0.013), and a significantly lower selective constraint by one of the two software packages (SnIPRE P = 0.022; MKTest P = 0.76). Thus, while copy number variation can potentially skew some parameter estimates, it is unlikely to be the driving force behind the statistical signal.

Another point that bears mention, is that the specimens used in this study came from populations introduced to Okinawa, and may represent only a subset of genetic diversity found in the native range populations. This will reduce number of polymorphic loci in the genome, and introduce a level of noise due to stochastic sampling from the source population. However, while these factors may change the absolute estimates of selective constraint and adaptive evolution, they should not affect the relative estimates within a genome, that is between venom and nonvenom components (fig. 2), or between different venom components (fig. 3).

Conclusions

An important characteristic of snake venoms is functional redundancy among components, resulting in multiple impacts upon critical prey systems, such as blood pressure regulation (Aird 2002; Aird etal. 2016). Redundancy is achieved both by employment of highly specific toxins that target particular receptors, ion channels, etc., and also through delivery and release of secondary messengers and modulators, such as purine nucleosides and polyamines, that affect multiple systems (Aird 2002). Individual venom components experience a wide variety of selective pressures, ranging from positive selection at certain sites, to purifying selection at others, to neutrality (supplementary fig. S5, Supplementary Material online). Our data suggest that selective regimes of venom components result from a complex adaptive landscape dominated by positive and negative selection, but also affected by varying degrees of genetic drift. In particular, changes in population size or population substructure, which affect genetic drift, may also affect trajectories of venom evolution. Although we have focused on understanding evolutionary forces acting on venom protein composition since the separation of two snake species, the increasing availability of genomes and genomic tools for other species will soon show whether these patterns may be more broadly applicable to other snakes, or perhaps even to all venomous organisms that employ peptidyl/proteinaceous secretions.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.