Skip to main content
Download PDF
  • Research Article
  • Open access
  • Published:

Trait differentiation and modular toxin expression in palm-pitvipers

BMC Genomics volume 21, Article number: 147 (2020) Cite this article

Abstract

Background

Modularity is the tendency for systems to organize into semi-independent units and can be a key to the evolution and diversification of complex biological systems. Snake venoms are highly variable modular systems that exhibit extreme diversification even across very short time scales. One well-studied venom phenotype dichotomy is a trade-off between neurotoxicity versus hemotoxicity that occurs through the high expression of a heterodimeric neurotoxic phospholipase A2 (PLA2) or snake venom metalloproteinases (SVMPs). We tested whether the variation in these venom phenotypes could occur via variation in regulatory sub-modules through comparative venom gland transcriptomics of representative Black-Speckled Palm-Pitvipers (Bothriechis nigroviridis) and Talamancan Palm-Pitvipers (B. nubestris).

Results

We assembled 1517 coding sequences, including 43 toxins for B. nigroviridis and 1787 coding sequences including 42 toxins for B. nubestris. The venom gland transcriptomes were extremely divergent between these two species with one B. nigroviridis exhibiting a primarily neurotoxic pattern of expression, both B. nubestris expressing primarily hemorrhagic toxins, and a second B. nigroviridis exhibiting a mixed expression phenotype. Weighted gene coexpression analyses identified six submodules of transcript expression variation, one of which was highly associated with SVMPs and a second which contained both subunits of the neurotoxic PLA2 complex. The sub-module association of these toxins suggest common regulatory pathways underlie the variation in their expression and is consistent with known patterns of inheritance of similar haplotypes in other species. We also find evidence that module associated toxin families show fewer gene duplications and transcript losses between species, but module association did not appear to affect sequence diversification.

Conclusion

Sub-modular regulation of expression likely contributes to the diversification of venom phenotypes within and among species and underscores the role of modularity in facilitating rapid evolution of complex traits.

Background

Modularity, the tendency for systems to organize into semi-independent discrete units, is a central theme in the evolution of biological systems and complex traits [1]. Modularity creates evolvability and the potential to adapt to novel environments rapidly by eliminating or reducing antagonistic pleiotropy while simultaneously permitting advantageous phenotypic changes through the use of conserved genetic machinery [2, 3]. Gene regulatory networks are an especially common mechanism for modular evolution within and among lineages [4]. Inducing, increasing, reducing, or eliminating expression of specific sub-modules can create or replicate advantageous phenotypes through the recombination of sub-modular features [5]. As such, modularity is a common characteristic of many adaptive traits because sub-module associated features can be rapidly modified without evolving ‘from scratch’ [2]. Heliconius butterflies provide a classic example where a variety of predator-deterring wing patterns have evolved and diversified through variation in modular elements (e.g., color and spot-pattern) controlled by just a few conserved genes (e.g., the optix transcription factor and the wntA signaling pathway) [5–7]. Identifying modules and their sub-modules underlying variation in highly variable modular traits can therefore provide valuable insight on the genetic basis of diversification across micro and macroscales.

Snake venoms are highly variable adaptive traits composed of 10–100 secreted proteins that collectively work to subdue prey or deter predation [8, 9]. Despite the perceived complexity of the venom system, venoms appear to evolve rapidly and respond to local selection pressures over short timescales [10, 11]. The exceptional degree of phenotypic variation observed in venoms can partially be contributed to the modularity of the venom system. Because toxin expression and production is localized to a specialized gland [12–15] (but see [16, 17]), the venom system is a functional module that is inherently more free to vary with limited pleiotropic effects. Moreover, venom functionality is, at least in part, dependent on the coordinated expression of specific toxins or toxin classes which may covary geographically or among species [18–20]. In many cases, recurrent patterns of variation in venom compositions suggest that expression of associated toxins represent sub-modules of variation, though empirical tests of sub-modularity of toxins are lacking.

One example of venom variation likely mediated by sub-modular regulation is an apparent phenotypic trade-off between neurotoxicity and hemotoxicity. In crotalid vipers (Viperidae: Crotalinae), hemorrhagic venoms are most common and are a function of high proportions of several toxin families, especially snake venom metalloproteinases (SVMPs) [21, 22]. However, in some lineages neurotoxicity has emerged as a principal phenotype [22]. An extremely well-documented manifestation of neurotoxicity in crotalid venoms is based on high expression of a heterodimeric β-neurotoxic phospholipase A2 (PLA2) complex [23, 24]. These phenotypes can manifest as interspecific, intraspecific, and/or ontogenetic variation [18–20, 22, 25–28], prompting the establishment of a “Type A/Type B" nomenclature to describe the variation in rattlesnakes. Type A venoms refer to those dominated by the neurotoxic PLA2s, and Type B venoms refer to those with high proportions of SVMPs. Notably, there are also descriptions of Type A+B venoms which have high proportions of neurotoxic PLA2s and hemorrhagic SVMPs, but these phenotypes are rare even in Type A - Type B contact zones [11, 19, 29]. Here, recurring phenotypic patterns, the lack of apparent phylogenetic signal (even over ecological time scales), and the usage of common genetic building blocks (i.e., toxin families) is suggestive of modularity mediating the evolution of these phenotypes.

An opportunity to test this exists in the arboreal pitvipers of the genus Bothriechis. One species, B. nigroviridis, exhibits a neurotoxic venom phenotype driven by the high abundance of a neurotoxic heterodimeric PLA2 named nigroviriditoxin [30, 31]. Bothriechis nigroviridis is unique among species with neurotoxic venom because of its ecological differentiation; B. nigroviridis is an arboreal high-elevation specialist while most others are mid-low elevation terrestrial species. The sister species to B. nigroviridis, B. nubestris, appears to occupy an extremely similar ecological niche based on its documented range and conserved morphology [32]. Although empirical studies of B.nubestris’ venom have yet to be conducted, its divergence from B. nigroviridis 6–10 mya would provide sufficient temporal opportunity for venom diversification [33]. Bothriechis nigroviridis and B. nubestris can therefore provide a test case for examining mechanisms of phenotypic diversification in a modular framework.

We sought to describe and compare the venom gland transcriptomes of B. nigroviridis and B. nubestris to understand toxin evolution in a modular framework. We characterize the venom gland transcriptomes of representatives of each species and identify key dimensions of variation within and between species. We identified conserved and unique toxins and used weighted-gene co-expression network analysis (WGCNA) to test for sub-modules of variation among distinct venom types. Based on the observation that neurotoxic and hemotoxic phenotypes occur independently, in combination, or as ontogenetic changes, we hypothesized that toxins associated with neurotoxic and hemorrhagic phenotypes (i.e., neurotoxic PLA2s and SVMPs) would segregate into distinct sub-modules of correlated expression variation. Additionally, we examine instances of intraspecific transcript duplication and loss and comparative sequence divergence. We hypothesized that if modular expression is a primary driver of variation, gene duplications and sequence diversification would be reduced in sub-module associated toxin families whose function has been selectively optimized and is primarily regulated by expression.

Results

Transcriptome characterization

To examine the evolutionary mechanisms underlying venom divergence we sequenced, assembled, and characterized the venom gland transcriptomes of two Bothriechis nigroviridis (CLP1856 and CLP1864) and two B. nubestris (CLP1859 and CLP1865) (Fig. 1, Table 1). The number of recovered toxins and recovered families were generally consistent with those of other viperid transcriptomes [25, 34–37] and with estimates of toxin family size in early high-throughput transcriptomes of B. schlegelii and B. lateralis [38] (Table 2, Table 3).

Fig. 1
figure 1

Phylogeny of Bothriechis based on [33] and a distribution map for B. nigroviridis and B. nubestris made in R v.3.5.3 (https://www.R-project.org/) based on ranges described in [74] and [33] and publicly available specimen localities in [32]. Sampled localities are shown as dots with specimen labels. Animal images were modified and used with permission from credit holder Alexander Robertson

Full size image
Table 1 Specimen information for Bothriechis individuals used in this work
Full size table
Table 2 Toxin transcripts recovered for Bothriechis nigroviridis and associated classifications as orthologs or paralogs, expected transcripts per million reads (TPM) estimated by RSEM, likely over expression classification as detected in intraspecific variation comparisons (i.e., above the 99th percentile of expected variance in expression based on a nontoxin null distribution), and coverage-based assessment of likely presence or absence
Full size table
Table 3 Toxin transcripts recovered for Bothriechis nubestris and associated classifications as orthologs or paralogs, expected transcripts per million reads (TPM) estimated by RSEM, over expression classification as detected in intraspecific variation comparisons (i.e., above the 99th percentile of expected variance in expression based on a nontoxin null distribution), and coverage-based assessment of likely presence or absence
Full size table

We recovered 1517 total transcripts for B. nigroviridis, which included 43 toxins from 13 toxin families. The venom transcriptome of B. nigroviridis was largely dominated by the expression of the heterodimeric neurotoxic PLA2, nigroviriditoxin [31], especially in the northern individual where it accounted for 60.3% of toxin expression (Fig. 2, Table 2). BPPs and SVSPs were also abundant in B. nigroviridis venoms, accounting for 7.6% and 14.6% of toxin expression, respectively (Fig. 2, Table 2). The high expression of the neurotoxic PLA2 complex observed in the northern individual is consistent with the neurotoxic phenotype previously described in individuals from a similar locality (∼50 km north of CLP1864’s locality, though from a different cordillera) [30] (Type A based on the rattlesnake nomenclature). Consistent with the Type A phenotype, there was low expression of CTL and SVMP variants which, in a previous proteomic study of B. nigroviridis, were not detected in the venom [30].

Fig. 2
figure 2

Venom characterization for Bothriechis nigroviridis. a Venom transcriptome compositions for B. nigroviridis based on average expression between two individuals. b Venom transcriptome compositions of each individual used. The venom of B. nigroviridis CLP1864 is largely consistent with the published proteome for this species. The high proportion of snake venom metalloproteinases (SVMPs) observed in the venom gland transcriptome of B. nigroviridis CLP1856 has not been described previously. c Intraspecific variation in transcript expression for B. nigroviridis. Data have been centered log-ratio transformed to account for their compositional nature. Dashed lines denote the 99% confidence interval of nontoxin expression and red lines are lines of best fit based on orthogonal residuals. B. nigroviridis displays substantially more variation in toxin expression, primarily in C-type lectins (CTLs), SVMPs, and snake venom serine proteinases (SVSP)s

Full size image

Unlike the northern B. nigroviridis, the southern B. nigroviridis showed substantial expression of the nigroviriditoxin subunits as well as SVMPs (Fig. 2, Table 2). Both subunits of nigroviriditoxin and seven of the nine SVMPS were identified as outliers in expression comparisons between the two individuals; nigroviriditoxin and one SVMP were found to be expressed outside of the 99th percentile of the null distribution in the northern B. nigroviridis while six SVMPs were expressed outside of the 99th percentile of the null distribution in the southern B. nigroviridis (Table 2). In addition to the toxin family differences, four CTL and 11 SVSP variants fell outside of the 99th percentile of the null distribution of expression divergence between individuals (Table 2). Of the 43 total toxins assembled for B. nigroviridis, 27 were expressed outside of the 99th percentile of the nontoxin null distribution. In many cases, expression differences could be explained by toxin absence. Overall, 14 toxins were found to be absent in one individual with six absences in the southern B. nigroviridis and eight absences in the northern B. nigroviridis. The overall pattern of toxin expression is more characteristic of a Type A+B phenotype than Type A [39].

For B. nubestris we recovered 1787 transcripts which included 42 toxins from 14 toxin families (Table 3). In contrast to B. nigroviridis, toxin expression and presences/absences were generally similar between the two sequenced individuals of B. nubestris (Fig. 3, Table 3). In total, 14 toxins were expressed outside of the 99th percentile of the nontoxin null distribution. Toxins whose expression was outside the 99th percentile spanned all major families including BPP, CTLs, PLA2s, SVMPs, and SVSPs. However, only two toxins, Bnube-BPP-1 and Bnube-SVMPIII-1, were found to be absent in one individual. The overall expression pattern for both individuals was broadly consistent with observed Type B venoms [18]. SVMPs and CTLs were highly abundant components in the venom making up, on average 34.9% and 40.4% of toxin expression, respectively. In addition to SVMPs and CTLs, B. nubestris also expressed three PLA2s at lower levels. Two of these PLA2s were orthologous to the alpha and beta subunits of nigroviriditoxin on average accounting for 0.2% and 0.5% of toxin expression, respectively. The third PLA2, Bnube-PLA2-3, made up 15.7% of toxin expression in one B. nubestris individual (CLP1865) and appears homologous to a non-enzymatic, myotoxic PLA2 in B. schlegelii [40, 41].

Fig. 3
figure 3

Venom characterization for Bothriechis nubestris. a Venom transcriptome compositions for B. nubestris based on average expression between two individuals for each species. b Venom transcriptome compositions of each individual used. The venom of B. nubestris is dominated by SVMPs and CTLs. c Intraspecific variation in transcript expression for B. nubestris. Data have been centered log-ratio transformed to account for their compositional nature. Dashed lines denote the 99% confidence interval of nontoxin expression and red lines are lines of best fit based on orthogonal residuals. The venoms of B. nubestris CLP1859 and CLP1865 are largely similar, though CLP1865 displays elevated expression of a basic PLA2 and BPPs

Full size image

Interspecific variation and submodule identification

OrthoFinder [42] identified 1282 one-to-one orthologs, which included 32 orthologous toxins. Due to the high variability in toxin expression observed between individuals of B. nigroviridis, we compared toxin expression of each individual to the average expression of B. nubestris (Fig. 4). High variation in ortholog expression was observed between the northern B. nigroviridis and B. nubestris, with 14 toxins detected as differentially expressed by DESeq2 (Fig. 4, Table 4). The most prominent pattern was the variation in expression of nigroviriditoxin subunits and SVMPs (Fig. 4); a pattern which supports the classification of the northern B. nigroviridis’ venom as Type A and B. nubestris’ venom as Type B. In contrast, only 8 orthologous toxins were detected as differentially expressed between the southern B. nigroviridis and B. nubestris (Fig. 4, Table 5). Moreover, the variance in orthologous expression between the southern B.nigroviridis and B. nubestris was substantially lower than observed in the previous comparison, due largely to increased expression of several SVMPs.

Fig. 4
figure 4

Interspecific comparisons of toxin expression between average Bothriechis nubestris toxin expression and a Type A B. nigroviridis and b Type A+B B. nigroviridis. TPM values have been centered log-ratio transformed to account for the compositional nature of the data. Dashed lines denote the 99% confidence interval of nontoxin expression and red lines are lines of best fit based on orthogonal residuals. Paralogs are shown near axes for each species

Full size image
Table 4 DESeq2 expression analyses for B. nigroviridis A versus B. nubestris toxins comparison
Full size table
Table 5 DESeq2 expression analyses for B. nigroviridis A + B versus B. nubestris toxins comparison
Full size table

We implemented WGCNA assigning three venom phenotypes as "treatments": Type A (B. nigroviridis CLP1864), Type A+B (B. nigroviridis CLP1856), and Type B (B. nubestris CLP1859 and CLP1865). After transcript filtering, 83 transcripts, including 22 toxin transcripts, were segregated into six modules (Fig. 5, in Additional file 1: Table S1). Most of the toxins associated with the Type A/Type B phenotypes segregated into two distinct modules. Module 2 contained five of the seven orthologous SVMPs while module 3 contained both nigroviriditoxin subunits. SVSPs were distributed across three modules, including module 2 and module 3. Similarly, BPPs were the only toxin assigned to module 1 which appeared to primarily capture intraspecific variation in B. nubestris. Of the three orthologous CTLs, one was removed during filtering and the remaining two were assigned to modules 2 and 6. Finally, two VEGFs were assigned to two separate modules as well. We did not identify any transcription factors associated with the putatively Type A or Type B modules. However, we did identify a translation initation factor, TIF-4E1, associated with module 2.

Fig. 5
figure 5

Expression profiles for the six expression modules identified by CEMiTool. Each line represents a transcript and its change in expression across treatments. Toxins assigned to each module are colored by class and labelled. Nontoxins associated with a module are shown as grey lines. Toxins generally associated with the Type A and Type B venom phenotypes (neurotoxic PLA2 subunits and SVMPs, respectively) largely separated into two modules: M2 and M3. B. nigroviridis with Type A+B venom showed generally intermediate expression of A-B associated toxins

Full size image

Gene family analyses

To better understand the dynamics of transcript turnover (i.e., gene duplications and transcript losses through either gene loss or gene silencing) in relation to families associated with specific modules, we inferred toxin family phylogenies for four highly expressed and diverse toxin families and identified species-specific gene duplication and transcript loss events. As expected, our results suggest that the majority of toxin genes in B. nigroviridis and B. nubestris were likely present in their common ancestor. In three of the four toxin families, OrthoFinder identified one-to-one orthologs for the majority of toxins, although expression levels were not necessarily conserved (Fig. 5). However, each toxin family exhibited at least one species-specific toxin loss and three of the families showed evidence of both losses and duplications.

Transcript turnover was lower in families with a higher proportion of toxins sorted into a specific submodule. The two CTLs were split between two expression submodules (M2 and M6) and had four deletions and one duplication. Similarly, five SVSPs were split between three modules with three SVSPs assigned to module 2. SVMPs were inferred to have a single duplication and loss and were similarly assigned to three modules (M2, M4, and M6), though the five consistently highly expressed SVMPs were assigned to M2. PLA2s were the only family to experience a single species-specific toxin transcript loss, and the two orthologous toxins were assigned to M3.

In both SVMPs and SVSPs, we observed sequence divergence that occurred in one or more toxin copies following a duplication event (Fig. 6). In the case of SVSPs, nucleotide sequence divergence was sufficient to give conflicting phylogenetic signal when compared to the amino acid-based phylogeny inferred by OrthoFinder (Fig. 6, in Additional File 1: Figure S1). Although we did not find a significant difference in expression of one-to-one toxin orthologs compared to duplicated or conserved toxins (p = 0.28), we did find a marginally significant interaction between species and expression of one-to-one orthologs versus duplicated or conserved toxins (p = 0.08, Fig. 7). More specifically, B. nubestris appeared to exhibit proportionately higher expression of toxins, but also disproportionately higher expression of duplicated and conserved toxins (Fig. 7).

Fig. 6
figure 6

Toxin family phylogenies and expression plots of a C-type lectins (CTLs), b phospholipase A2s (PLA2s), c snake venom metalloproteinases (SVMPs), and (d) snake venom serine proteases (SVSPs). Single copy toxin orthologs identified by OrthoFinder are marked by brackets in the phylogeny. Toxin transcript gains and losses were inferred based on a simple parsimony model and are shown on phylogenies as grey circles and rectangles, respectively. Expression plots are based on average expression of each toxin for each species and dashed lines denote 99% confidence interval established by nontoxin expression. Identified orthologs are shown as colored circles and losses as colored inverted triangles. Duplicated toxins are shown as colored diamonds and expression of each copy is plotted against expression of their orthologous counter part in the other species (identified with bracketing on plots)

Full size image
Fig. 7
figure 7

Violin plots comparing expression of orthologous and paralogous toxins for Bothriechis nigroviridis and B. nubestris. Orthologous and paralogous toxins were not differentially expressed between the species

Full size image

Sequence based selection analyses

To determine the extent and role of sequence diversification in differentiating venoms, we compared pairwise values of ω, dS, and dN between toxin and nontoxin orthologs. Toxin sequences exhibited significantly higher values of ω (p <0.001) with three toxins, CTL-2, SVMPII-1, and SVMPIII-5, having ω values >1 indicating positive selection (Fig. 7). Despite having a higher ω ratio than the background nontoxins, the overall mean ω for toxin sequences was 0.56. Additionally, we tested for differences in synonymous and nonsynonymous substitution rates between toxins and nontoxins under the expectation that toxins and nontoxins should display similar background synonymous substitution rates but differ in nonsynonomous substitutions resulting in diversifying selection. As expected, we found no differences in synonymous substitution rates between toxins and nontoxins (p = 0.252) but significantly higher nonsynonymous substitution rates (p <0.001). Moreover, nine toxins had nonsynonymous substitution above the 95th percentile of nontoxin sequences; nearly double the number of toxins above the 95th percentile of ω. However, four of these toxins were found to have synonymous substitution above the 95th percentile of nontoxin sequences.

Discussion

We tested the hypothesis that dimensions of the neurotoxic-hemorrhagic venom phenotype were associated with specific submodules of toxin expression. We identified six submodules of expression variation, which included a primarily Type A submodule containing both nigroviriditoxin homolog subunits and a primarily Type B submodule containing the majority of orthologous SVMPs. The findings supported our hypothesis and implicate submodular regulation as a mechanism for rapid venom diversification. Modular expression regimes would allow rapid transitions between phenotypes while avoiding or minimizing occurrence of low-fitness intermediates [2] and facilitate ontogentic shifts observed in many groups [27, 28, 43, 44]. In the Bothriechis system, modularity effectively explains many of the toxin expression differences between B. nigroviridis and B. nubestris. The patterns of modularity observed here are also consistent with on-going genomic research to elucidate the genomic architecture mediating venom phenotype evolution [15, 45, 46]. Taken together these findings provide strong support for a role of sub-modular variation mediating changes in snake venom phenotypes.

Modularity underlying the neurotoxic-hemorrhagic dichotomy

The patterns of modularity and submodular organization inferred by WGCNA analyses explained much of the inter- and intraspecific variation in toxin expression we observed for B. nigroviridis and B. nubestris. We recovered a venom gland transcriptome for the northern B. nigroviridis consistent with the published proteomic venom phenotype and Type A venom expression. The increase in expression of nigroviriditoxin/nigroviriditoxin homologs is accomplished primarily through modification of regulatory patterns in module 3. Similarly, modifications to regulatory elements in module 2 can mediate expression regime shifts of many toxins, especially SVMPs. The strong association of these modules with species-specific patterns of inheritance demonstrate how modularity can promote rapid phenotypic transition among recently diverged and/or eco-morphically conserved species.

Of note was the Type A+B expression pattern in the southern B. nigroviridis which suggested intermediate or combined expression of the Type A and Type B submodules. Although Type A+B venoms have been documented in multiple species [19, 39] they are primarily associated with species exhibiting population level neurotoxic-hemorrhagic dichotomies and often occur at lower frequency than either the Type A or Type B phenotypes [11]. If this pattern holds true in B. nigroviridis, it would suggest the existence of individuals or populations of B. nigroviridis that have primarily Type B venom. Population level sampling has been difficult to attain due to the inherent rarity of this species and the logistical challenges of sampling many of the undisturbed, high-elevation regions of its distribution. However, population level sampling will be key for understanding the ecological and evolutionary dynamics of venom variation in this species. More importantly, the occurrence of the Type A+B phenotype in B. nigroviridis and other species suggests that the Type A and B submodules are not mutually exclusive. Rather, each module likely has independent genetic architectures which can occur in variable combinations among populations and species.

Modular expression effectively explains Type A/Type B toxin variation among these two species, but several toxin families such as CTLs, SVSPs, and VEGFs did not fit this framework. The variation observed in these families underscores the diversity of expression patterns in venom toxins and presents an ongoing challenge for the future. Although a great deal of work has been devoted to dissecting broad patterns of venom variation (e.g., neurotoxic-hemorrhagic dichotomy), the mechanisms influencing variation in other diverse toxin families such as SVSPs and CTLs deserves further investigation.

While our findings present evidence for submodularity of toxin expression, it is important to note their limitations as well. WGCNA identifies submodular clusters based on positive and negative correlations in transcript expression across assigned treatments with the expectation that these transcripts may be influenced by common regulatory elements. Because coexpression network analyses are based on observed patterns of expression rather than experimental validation, they are better regarded as hypotheses of submodular association rather than empirical findings. Moreover, WGCNA are ideally implemented using thousands of candidate transcripts derived from thoroughly assembled and annotated genomes with tens of replicates across treatments for robust inference. Unfortunately, genomic resources remain limited for snakes and such large sample sizes are difficult to attain for many species. Here, we have implemented WGCNA with a much reduced sample size and far fewer candidate genes than is typically ideal, which may make module assignment less powerful and robust, especially for lowly expressed transcripts. Nevertheless, our analyses assigned many highly expressed toxins to biologically plausible submodules corresponding to known axes of phenotypic variation in snake venom. Thus, we believe that WGCNA as implemented here represent an important proof-of-concept in the relevance and potential of these methods and the conceptual framework of modularity for evolutionary study of venom differentiation.

Mechanisms promoting modularity

Although our WGCNA and similar approaches identify submodules of variation based on phenomenological rather than mechanistic models, observed patterns of expression and recent genomic work implicate several general mechanisms contributing to modularity of the system. For instance, one of the primary advantages of co-expression network approaches is the ability to identify regulatory components such as transcription factors that potentially mediate the identified expression differences. In sub-module 2, we identified one translation initiation factor that showed increased expression with progression towards the Type B phenotype. Translation initiation factors enhance translation by stabilizing mRNA and facilitating assembly of ribosomal complexes [47]. In mammals, TIF-4E is required for efficient translation and acts as a translational regulatory mechanism [47]. Here, its association with module 2 may reflect an effort to promote rapid translation of the relatively large and highly expressed SVMPs. Though concordant expression of TIF-4E and module 2 toxins does not necessarily imply a causative link, it does present a hypothesis to test through functional validation.

The identification of primarily neurotoxic and hemorrhagic submodules are also consistent with recent genomic evidence which show that Type A and Type B toxins are inherited as independent haplotypes [15, 45, 46]. In some cases, presence and absence differences in these genes have been implicated as the primary drivers of variance in Type A/Type B phenotypes. In the case of the northern B. nigroviridis, absence of the SVMP tandem array could account for both the low expression of SVMPs and their inferred absence from the transcriptome (Table 2). In contrast, both B. nubestris individuals express low levels of a nigroviriditoxin homolog. Despite patterns of low expression, the sequences of the B. nubestris PLA2s were highly conserved with respect to nigroviriditoxin; both subunits had over 99% nucleotide sequence similarity with three nonsynonomous substitutions occurring in the beta subunit and one synonymous substitution occurring in the alpha subunit. The conservation of these sequences suggests that the B. nubestris variants of nigroviriditoxin have likely retained their neurotoxic function and that convergence on a "low neurotoxicity" phenotype therefore occurs through regulatory evolution in Bothriechis rather than through gene loss/gain as is observed in other species [15, 45, 46].

If expression patterns of the Type A and Type B submodules are inherited as independent haplotypes with additive effects, we can hypothesize that combined phenotypes are possible and should exhibit intermediate expression of of each module. The expression patterns apparent in the southern B. nigroviridis support these predictions as it displayed intermediate expression between the Type A B. nigroviridis and the Type B B. nubestris for the majority of Type A and Type B associated toxins. Additive expression of species-specific toxins has also been observed in interspecific hybrids where the putatively heterozygous offspring exhibit lower expression levels than presumably homozygous parents [35]. In the case of B. nigroviridis, intermediate expression observed in the southern B. nigroviridis could feasibly be the result of heterozygosity at Type A and Type B loci, though such a hypothesis is largely postulation without genomic evidence. As such, comparative genomics approaches that test architectural mechanisms promoting and mediating modularity are a promising avenue for future work.

Transcript turnover and diversification in a modular system

We expected selective optimization for modularity of toxin expression to affect toxin transcript turnover and sequence diversification. We tested for these effects in four toxin families and found that although all four toxin families had experienced some turnover, rates of duplication and loss were higher in toxins less associated with specific modules. Many snake toxin families have experienced dramatic expansions since their common ancestor [9] though the frequency of toxin duplications and losses within species is not clear. The marginal decrease in transcript-turnover with increased association with a specific submodule suggests selection for maintaining these toxins. Duplications are often implicated as having a primary role in toxin neofunctionalization by creating functional redundancy that allows toxins to ‘explore’ the phenotype space [9, 48, 49], but can also occur as a mechanism to increase expression of beneficial toxins [50]. We observed both increased sequence divergence following duplication and a marginal increase in expression of duplicated or conserved (i.e., not deleted or silenced) toxins specific to the B. nubestris lineage. Whether the possible emphasis on expression of paralogous versus orthologous toxins reflect phenomena unique to the B. nubestris lineage or a broader trend in the evolution of the more complex, hemorrhagic venom types is not clear, especially given our limited sample size. However, increased sampling of lineages and their toxin compositions will provide improved resolution to test the extent and role of gene duplication and loss in venom diversification.

We expected sequence diversification to be lowest in module associated toxins, but we did not find evidence to support this. Two of the three toxins with ω above one were SVMPs associated with Module 2, suggesting that although regulation may be conserved/coordinated, functionality is not. Many of the toxins with elevated rates of nonsynonymous substitution had similarly high rates of synonymous substitutions, which may indicate an overall higher substitution rate than the genomic background. Notably, SVSPs, which were generally less associated with a specific module, displayed some of the highest values of both dN and dS. The overall elevated substitution rates of these toxins and the lack of correspondence to clear expression regimes may reflect higher rates of substitution and recombination in these gene regions, though patterns of gene expression and the organization of the genetic architecture of SVSP regions is not well understood. Overall, toxin ω values were generally below what is expected under positive selection with just a few toxins displaying ω values greater than 1. Instead, toxin evolution between species appears to function under a model of relaxed purifying selection, which has been similarly noted in other interspecific comparisons of toxin sequence evolution [20].

Conclusions

Snake venoms are key innovations that have allowed the diversification of species across the globe. Unfortunately, many of the genomic mechanisms governing rapid variation of phenotypes remain uncertain. Through comparative transcriptomics and coexpression network analyses, we demonstrated how rapid transition between a common phenotypic venom dichotomy can occur through submodular regulation of the associated toxins. Modularity of the venom system and submodular variation of venom classes likely contribute to broader patterns of variation observed across taxonomic levels [51]. As genomic and transcriptomic resources become more available for venomous snakes, systems-based approaches such as the co-expression network analyses used here will yield more comprehensive understanding of the evolution of venoms and other complex, modular traits. Although our work presents these findings in the limited context of a single species pair, it highlights the importance of considering how complex traits function and evolve as a modular system. Our understanding of the selective forces that generate modularity and how modularity in turn mediates and facilitates the evolution of complex traits remains incomplete. However, as we have shown here, on-going efforts to address these questions in dynamic adaptive systems can provide key insights that lead to a more integrated understanding of the genomics of rapid adaptation in complex traits.

Methods

Sample collection

We collected two individuals of Bothriechis nigroviridis and two B. nubestris in May-June of 2016 for venom gland extraction and sequencing. Due to the smaller range of B. nubestris, both individuals were collected from the same locality (∼1 km apart), San Gerardo de Dota, San Jose province, Costa Rica. Bothriechis nigroviridis occupies a wider range than B. nubestris and we collected two individuals from distant populations. One of these individuals (CLP1864), was collected from outside of the La Esperanza sector of Parque Tapanati, Cartago province, Costa Rica, a locality that is ∼50 km south of specimens collected and used in previous proteomic studies characterizing the venom of this species [30]. The second individual (CLP1856) came from the southern most portion of the species’ range in Costa Rica, Las Tablas, Puntarenas province, Costa Rica (Fig. 8) ∼200 km southeast of specimens used in [30].

Fig. 8
figure 8

Distribution of a pairwise dN/dS ratios, b synonymous substitution rates, and c nonsynonymous substitution rates of orthologous transcripts. Dashed red lines denote 95 percentiles based on distribution of nontoxins. Lines beneath plots indicate toxins, and toxins with values greater than the 95 percentile are marked with blue arrows. In c, toxins above the 95th percentile with elevated synonymous mutation rates (i.e., above the 95th percentile in b are colored yellow. Toxins had statistically higher dN/dS ratios and nonsynonymous substitution rates based on a Wilcoxon signed rank test. Toxin and nontoxin synonymous mutation rates were not significantly different

Full size image

Following collection, each individual had its venom collected via manual extraction. Collected venoms were lyophilized and stored at -20 C for later use. Each animal was sacrificed four days later when transcription of venom proteins is at its maximum [52], via injection of sodium pentobarbitol (100mg/kg). Venom glands were dissected and stored separately in approximately 2 mL of RNAlater preservative. Animal carcasses were preserved as museum specimens with 10% buffered formalin and deposited in the Universidad de Costa Rica. The above methods were approved by University of Central Florida Institutional Animal Care and Use Committee (IACUC) protocol 16-17W, Clemson University IACUC protocol number 2017-067, and Universidad de Costa Rica Comimté Institucional para el Cuidado y Uso de los Animales (CICUA) permit number CICUA-082-17.

Venom gland transcriptome sequencing

Total RNA was extracted from left and right glands independently using a standard, Trizol reagent extraction as described in [53]. Briefly, diced venom gland tissues were submerged in 500 μL of Trizol, homogenized with a sterile 20-gauge needle, and treated with an additional 500 μL of Trizol and 200 μL chloroform. RNA was then separated from tissue, cellular components, and DNA by centrifuging the total mix in a 5Prime phase lock gel heavy tube for 20 minutes at 12,000 g. Supernatant containing the RNA was transferred to a new tube and RNA was precipitated with 500 μL of isopropyl alcohol. Pelleted RNA was washed in 75% ethanol and re-suspended in RNAase free water. Extracted total RNA was checked for quality and quantified using either an Agilent 2100 Bioanalyzer or Agilent 2200 TapeStation and stored at -80 C.

We prepared cDNA libraries from 1 μL high quality total RNA using the NEBNext Ultra RNA library Prep Kit for Illumina following the manufacturer’s instructions. Specifically, we isolated polyadenalated RNA with the NEB Poly(A) Magnetic Isolation Module (New England Biolabs) and fragmented resulting mRNA by heat fragmentation at 70 ∘ C for 14.5 minutes to attain an average size of approximately 370 bp. mRNA fragments were reverse transcribed to cDNA and each library was ligated with a unique combination of index primers and Illumina adapters. The cDNA libraries were amplified through PCR using the NEBNext High-Fidelity 2X Hot Start PCR Master Mix and 14 cycles of PCR. Amplified cDNA was purified with Agencourt AMPure XP PCR Purification beads. The resulting libraries were checked for quality, fragment size distribution, and concentration on either an Agilent 2100 Bioanalyzer or Agilent 2200 TapeStation. KAPA qPCR was additionally performed on each sample library to determine amplifiable concentrations. Libraries were then pooled in groups of twelve with equal representation of amplifiable cDNA for sequencing.

Sequencing took place on an Illumina HiSeq 2000 at the Florida State University College of Medicine’s Translational Science Laboratory. Combined libraries were multiplexed and sequenced with a 150 bp paired-end rapid run lane. Raw reads were demultiplexed and quality checked in FastQC [54]. To account for reads which may have been mis-assigned during demultiplexing, we used jellyfish v.2.2.6 [55] and KAT v.2.3.4 [56] to identify and filter reads with kmers that exhibited more than a 500 fold difference in occurrence between samples sequenced on the same lane. Adapter sequences and low quality bases were then trimmed using trim-galore v.0.4.4 [57]. Finally, to increase both quality and total length of read sequences, we used PEAR v 0.9.6 [58] to merge paired reads with a 3’ overlap of greater than 10 bp.

Transcriptome assembly and analyses

Previous transcriptome studies have shown the challenges associated with venom gland transcriptome assembly, due to the contrast in a proportionately low number of highly expressed toxin transcripts compared to the much broader, low expression of house keeping genes [59]. To overcome this, we performed three independent assemblies using Extender [53], the DNAstar NGen assembler v.15.0, and Trinity v.2.4.0 [60] per the strategy suggested in Holding et al. [59]. Sequence identities of toxins from each assembly were identified via local blastx search of SWISS-prot’s curated toxins database. Contigs with a blast match of greater than 90% identity were then clustered against a database of identified snake toxins to annotate coding regions of 90% similarity or greater. Coding regions of remaining toxin contigs were annotated manually in Geneious v.10.2.3 [61]. Contigs which were not identified as toxins were annotated by clustering against a database of previously identified snake nontoxins to annotate coding regions of 90% similarity or greater representing nontoxin transcripts used in later analyses. Annotated transcripts from independent assemblies were combined and duplicate sequences as well as coding regions with ambiguous sites were removed. The remaining transcripts were screened for chimeric or mis-assembled coding sequences by mapping merged reads against these sequences with bwa v.0.7.16 [62] and checking for uneven read distribution across sites. Specifically, sequences with sites where the mean number of bases per read on either side of a site differed by more than 50% of the mean read length were considered likely chimeras, checked manually, and removed accordingly. We clustered the remaining transcripts at a threshold of 98% similarity to account for toxin alleles or recent paralogs that may be present. This represented the final transcriptome for each individual. To account for variation among individuals in a species and for stochastic variation in the assembly process that may have resulted in failure to assemble specific toxins in a given individual, we combined final contig sets for individuals of the same species, removed duplicates, and clustered coding regions of 98% similarity to create a master transcriptome for each species. These species-specific master transcriptomes were then used for subsequent read mapping and expression analyses.

Expression analyses and ortholog identification

To determine relative expression of transcripts, we mapped reads from individuals to their species master transcriptome with Bowtie2 v2.3.2 and calculated relative expression with RSEM v.1.3.0 [63]. Intraspecific differences in expression were assessed using species-specific datasets for B. nigroviridis and B. nubestris. Because our limited intraspecific sampling precluded formal tests for differential expression within species, we generated pairwise null distributions of expression divergence for each species based on nontoxin expression to identify outlier toxins similar to [64]. Data were first centered log-ratio (clr) transformed to normalize the expression distributions while accounting for the compositional nature of relative expression values (e.g., TPM) using the cmultRepl function in the R package zCompositions [25, 65, 66]. Toxins whose pairwise divergence in expression fell outside the 99th percentile of the centered log-ratio transformed nontoxin distributions were considered outliers that are likely differential expression. RSEM can assign non-zero values to transcripts that may not be present in the transcriptome through mis-mapping of reads from other transcripts with regions of high similarity. To verify the extent to which toxins varied in presence or absence within species we aligned merged reads to the species-specific transcript sets to screen for poor read mapping. Toxins that had regions greater than 10% of the total sequence length with less than 5x coverage or highly anomalous read distributions (determined by manual review) were considered absent in the transcriptome of a given individual.

Toxin families in snakes are notorious for undergoing rapid expansions and losses, which is problematic for interspecific comparisons which assume orthology among matched transcripts. To overcome this we identified orthologous groups of transcripts using OrthoFinder v.2.3.1 [42] specifying multisequence alignments with mafft. OrthoFinder identifies groups of sequences derived from a single gene in the common ancestor of compared species (i.e., orthogroups), as well as identifies conserved orthologs within orthogroups. We classified transcripts as orthologs or paralogs by parsing the OrthoFinder “orthologs" output to identify single copy orthologs and one-to-one orthologs within orthogroups using a custom python script (orthocombiner.py). For interspecific comparisons, expression data for orthologous and paralogous transcripts were combined into a single dataset where paralogous transcripts were given an expression value of zero where absent for a given species. We used estimates of read counts from RSEM to test for differences in transcript expression with DESeq2 in R v.3.5.3 [67].

Network analyses

We performed weighted gene coexpression network analysis using the R package CEMitool [68] in R. A variance stabilizing transformation (vst) was used and transcripts were filtered to reduce correlation between variance and gene expression. We used pearson’s coefficient as the correlation method and a beta value of 10 was automatically selected. The minimum module size was set to 1 to allow the greatest flexibility in identifying modules of correlated expression. Because of the high variability in venom composition observed among B. nigroviridis (see above), we annotated samples as one of three venom types which correspond to venom phenotypes observed in rattlesnakes: B. nigroviridis Type A (CLP1864), B. nigroviridis Type A+B (CLP1856), and B. nubestris type B (CLP1859 and CLP1865).

Gene family analyses

To more closely examine how toxin family expansion, duplications, and loss have shaped venom composition, we constructed phylogenies for the four most highly expressed toxin families: C-type lectins (CTLs), PLA2s, snake venom serine proteases (SVSPs), and SVMPs. Alignments for each family were generated with mafft v.7.407 [69] and checked manually in Geneious. Partitioning schemes for each gene family were determined using PartitionFinder v.2 [70]. Phylogenies were then recovered with MrBayes v.3.2.6 [71]. MrBayes was run using one cold and three heated chains for 10 million generations with a variable rate prior. We then identified and mapped species-specific deletion and duplication events onto the trees based on the output of OrthoFinder. We considered toxins that were unassigned an ortholog to be indicative of gene loss in one species while one to many ortholog assignments indicated duplications within a species. We tested for differences in expression of one-to-one orthologs versus conserved and duplicated toxins with a two-way factorial with toxin type and species as factors in R. TPM values were used as the metric for expression and were centered log-ratio transformed to linearize the data while preserving their compositional nature [25, 65].

Sequence analyses

We compared divergence of orthologous toxin and nontoxin transcripts by calculating dN/dS ratios (ω). Orthologous transcripts were first aligned by codon using PRANK v.170427 [72]. PRANK alignments were then used as input to estimate ω, dS, and dN with codeml in paml v. 4.9 [73].

We compared ω, dS, and dN of toxin genes against a background of nontoxins as in [20] to discern if toxin genes exhibited higher synonymous and/or nonsynonymous substitution rates and if toxins displayed high rates of positive selection (i.e., higher values of ω). We excluded sequences with dS <0.001 due to the possibility of estimating excessively inflated values of ω, and sequences with dS >0.10 to reduce the risk of including misidentified orthologs. Statistical differences in ω, dS, and dN values between toxins and nontoxins were tested with a wilcoxon sign rank test in R.

Availability of data and materials

Raw sequence data and transcript sequences generated during the current study are available on the National Center for Biotechnology Information (NCBI) under accession numbers given in Table 1. Consensus transcriptomes have been to submitted to the NCBI Transcriptome Shotgun Assembly (TSA) database under GIBL00000000 (Bothriechis nigroviridis) and GIBM00000000 (B. nubestris). Scripts used in data analyses are available on GitHub at: https://github.com/masonaj157/Trait_differentiation_and_modular_expression_in_palm-pitvipers

Abbreviations

BPP:

Bradykinin potentiating peptide

CTL:

C-type lectin

PLA2 :

Phospholipase A2

SVMP:

Snake venom metalloproteinases

SVSP:

Snake venom serine proteinase

TIF:

Translation initiation factor

VEGF:

Snake venom vascular endothelial growth factor

WGCNA:

Weighted gene coexpression network analysis

References

  1. Wagner GP, Pavlicev M, Cheverud JM. The road to modularity. Nat Rev Genet. 2007; 8(12):921–31.

    Article  CAS  PubMed  Google Scholar 

  2. von Dassow G, Munro E. Modularity in animal development and evolution: elements of a conceptual framework for EvoDevo. J Exper Zool. 1999; 285(4):307–25.

    Article  CAS  Google Scholar 

  3. Yang AS. Modularity, evolvability, and adaptive radiations: a comparison of the hemi- and holometabolous insects. Evol Dev. 2001; 3(2):59–72.

    Article  CAS  PubMed  Google Scholar 

  4. Levine M, Davidson EH. Gene regulatory networks for development. Proc Nat Acad Sci. 2005; 102(14):4936–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Ferguson LC, Maroja L, Jiggins CD. Convergent, modular expression of ebony and tan in the mimetic wing patterns of Heliconius butterflies. Develop Genes Evol. 2011; 221(5-6):297–308.

    Article  Google Scholar 

  6. Joron M, Jiggins CD, Papanicolaou A, McMillan WO. Heliconius wing patterns: an evo-devo model for understanding phenotypic diversity. Heredity. 2006; 97(3):157–67.

    Article  CAS  PubMed  Google Scholar 

  7. Van Belleghem SM, Rastas P, Papanicolaou A, Martin SH, Arias CF, Supple MA, Hanly JJ, Mallet J, Lewis JJ, Hines HM, Ruiz M, Salazar C, Linares M, Moreira GRP, Jiggins CD, Counterman BA, McMillan WO, Papa R. Complex modular architecture around a simple toolkit of wing pattern genes. Nat Ecol Evol. 2017; 1(3):0052.

    Article  PubMed Central  Google Scholar 

  8. Mackessy SP. Handbook of Venoms and Toxins of Reptiles. Boca Raton, FL: CRC press; 2016.

    Book  Google Scholar 

  9. Casewell NR, Wüster W, Vonk FJ, Harrison Ra, Fry BG. Complex cocktails: the evolutionary novelty of venoms. Trends Ecol Evol. 2013; 28(4):219–29.

    Article  PubMed  Google Scholar 

  10. Margres MJ, Wray KP, Seavy M, McGivern JJ, Herrera ND, Rokyta DR. Expression differentiation is constrained to low-expression proteins over ecological timescales. Genetics. 2016; 202(1):273–83.

    Article  CAS  PubMed  Google Scholar 

  11. Strickland JL, Smith CF, Mason AJ, Schield DR, Borja M, Castañeda-Gaytán G, Spencer CL, Smith LL, Trápaga A, Bouzid NM, Campillo-García G, Flores-Villela OA, Antonio-Rangel D, Mackessy SP, Castoe TA, Rokyta DR, Parkinson CL. Evidence for divergent patterns of local selection driving venom variation in Mojave Rattlesnakes (Crotalus scutulatus). Sci Rep. 2018; 8(1):17622.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  12. Oron U, Bdolah A. Regulation of protein synthesis in the venom gland of viperid snakes. J Cell Biol. 1973; 56:177–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Fry BG, Casewell NR, Wüster W, Vidal N, Young B, Jackson TNW. The structural and functional diversification of the Toxicofera reptile venom system. Toxicon. 2012; 60(4):434–48.

    Article  CAS  PubMed  Google Scholar 

  14. Vonk FJ, Casewell NR, Henkel CV, Heimberg AM, Jansen HJ, McCleary RJR, Kerkkamp HME, Vos R. a., Guerreiro I, Calvete JJ, Wüster W, Woods AE, Logan JM, Harrison Ra, Castoe Ta, de Koning aPJ, Pollock DD, Yandell M, Calderon D, Renjifo C, Currier RB, Salgado D, Pla D, Sanz L, Hyder AS, Ribeiro JMC, Arntzen JW, van den Thillart GEEJM, Boetzer M, Pirovano W, Dirks RP, Spaink HP, Duboule D, McGlinn E, Kini RM, Richardson MK. The king cobra genome reveals dynamic gene evolution and adaptation in the snake venom system,. Proc Nat Acad Sci USA. 2013; 110(51):20651–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Schield DR, Card DC, Hales NR, Perry BW, Pasquesi GM, Blackmon H, Adams RH, Corbin AB, Smith CF, Ramesh B, Demuth JP, Betrán E, Tollis M, Meik JM, Mackessy SP, Castoe TA. The origins and evolution of chromosomes, dosage compensation, and mechanisms underlying venom regulation in snakes. Genome Res. 2019; 29(4):590–601.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Hargreaves AD, Swain MT, Hegarty MJ, Logan DW, Mulley JF. Restriction and recruitment—gene duplication and the origin and evolution of snake venom toxins. Genome Biol Evol. 2014; 6(8):2088–95.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Reyes-Velasco J, Card DC, Andrew AL, Shaney KJ, Adams RH, Schield DR, Casewell NR, Mackessy SP, Castoe TA. Expression of venom gene homologs in diverse python tissues suggests a new model for the evolution of snake venom. Mol Biol Evol. 2015; 32(1):173–83.

    Article  CAS  PubMed  Google Scholar 

  18. Glenn JL, Straight R. Mojave rattlesnake Crotalus scutulatus scutulatus venom: variation in toxicity with geographical origin. Toxicon. 1978; 16(1):81–4.

    Article  CAS  PubMed  Google Scholar 

  19. Glenn JL, Straight RC, Wolt TB. Regional variation in the presence of canebrake toxin in Crotalus horridus venom. Comparative Biochem Physiol Part C: Pharmacol, Toxicol Endocrinol. 1994; 107(3):337–46.

    CAS  Google Scholar 

  20. Rokyta DR, Wray KP, Margres MJ. The genesis of an exceptionally lethal venom in the timber rattlesnake (Crotalus horridus) revealed through comparative venom-gland transcriptomics. BMC Genomics. 2013; 14:394.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Soto JG, Perez JC, Minton SA. Proteolytic, hemorrhagic and hemolytic activities of snake venoms. Toxicon. 1988; 26(9):875–82.

    Article  CAS  PubMed  Google Scholar 

  22. Mackessy SP. Venom composition in rattlesnakes: trends and biological significance In: Hayes WK, Beaman KR, Cardwell MD, Bush SP, editors. The Biology of Rattlesnakes. Loma Linda, CA: Loma Linda University Press: 2008. p. 495–510.

    Google Scholar 

  23. Doley R, Zhou X, Kini RM. Snake venom phospholipase A2 enzymes. Handbook of venoms and toxins of reptiles. 2010; 1:173–205.

    Google Scholar 

  24. Gutiérrez JM, Lomonte B. Phospholipases A2: unveiling the secrets of a functionally versatile group of snake venom toxins. Toxicon. 2013; 62:27–39.

    Article  PubMed  CAS  Google Scholar 

  25. Rokyta DR, Margres MJ, Calvin K. Post-transcriptional mechanisms contribute little to phenotypic variation in snake venoms,. G3 (Bethesda, Md.) 2015; 5(11):2375–82.

    Article  CAS  Google Scholar 

  26. Fernández R, Edgecombe GD, Giribet G. Exploring phylogenetic relationships within Myriapoda and the effects of matrix composition and occupancy on phylogenomic reconstruction. Syst Biol. 2016; 65(5):871–89.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Mackessy SP, Williams K, Ashton KG. Ontogenetic variation in venom composition and diet of Crotalus oreganus concolor: a case of venom paedomorphosis?Copeia. 2003; 2003(4):769–82.

    Article  Google Scholar 

  28. Calvete JJ, Sanz L, Cid P, de la Torre P, Flores-Díaz M, Dos Santos MC, Borges A, Bremo A, Angulo Y, Lomonte B, Alape-Girón A, Gutiérrez JM. Snake venomics of the Central American Rattlesnake Crotalus simus and the South American Crotalus durissus complex points to neurotoxicity as an adaptive paedomorphic trend along Crotalus Dispersal in South America. J Proteome Res. 2010; 9(1):528–44.

    Article  CAS  PubMed  Google Scholar 

  29. Glenn JL, Straight RC. Intergradation of two different venom populations of the Mojave Rattlesnake (Crotalus scutulatus scutulatus) in Arizona. Toxicon. 1989; 27(4):411–8.

    Article  CAS  PubMed  Google Scholar 

  30. Fernández J, Lomonte B, Sanz L, Angulo Y, Calvete JJ. Snake venomics of Bothriechis nigroviridis reveals extreme variability among palm pitviper venoms : different evolutionary solutions for the same trophic purpose. J Proteome Res. 2010; 9:4234–41.

    Article  PubMed  CAS  Google Scholar 

  31. Lomonte B, Mora-Obando D, Fernández J, Sanz L, Pla D, María Gutiérrez J, Calvete JJ. First crotoxin-like phospholipase A2 complex from a New World non-rattlesnake species: Nigroviriditoxin, from the arboreal Neotropical snake Bothriechis nigroviridis. Toxicon. 2015; 93:144–54.

    Article  CAS  PubMed  Google Scholar 

  32. Doan TM, Mason AJ, Castoe TA, Sasa M, Parkinson CL. A cryptic palm-pitviper species (Squamata: Viperidae: Bothriechis) from the Costa Rican highlands, with notes on the variation within B. nigroviridis. Zootaxa. 2016; 4138(2):271–90.

    Article  PubMed  Google Scholar 

  33. Mason AJ, Grazziotin FG, Zaher H, Lemmon AR, Moriarty Lemmon E, Parkinson CL. Reticulate evolution in nuclear Middle America causes discordance in the phylogeny of palm-pitvipers (Viperidae: Bothriechis). J Biogeography. 2019; 46(5):833–44.

    Article  Google Scholar 

  34. Rokyta DR, Lemmon AR, Margres MJ, Aronow K. The venom-gland transcriptome of the Eastern Diamondback Rattlesnake (Crotalus adamanteus). BMC Genomics. 2012; 13(1):312.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Aird SD, Aggarwal S, Villar-Briones A, Tin MM-Y, Terada K, Mikheyev AS. Snake venoms are integrated systems, but abundant venom proteins evolve more rapidly. BMC Genomics. 2015; 16(647):1–20.

    Google Scholar 

  36. Strickland J, Mason A, Rokyta D, Parkinson C, Strickland JL, Mason AJ, Rokyta DR, Parkinson CL. Phenotypic variation in Mojave Rattlesnake (Crotalus scutulatus) venom is driven by four toxin families. Toxins. 2018; 10(4):135.

    Article  PubMed Central  CAS  Google Scholar 

  37. Hofmann EP, Rautsaw RM, Strickland JL, Holding ML, Hogan MP, Mason AJ, Rokyta DR, Parkinson CL. Comparative venom-gland transcriptomics and venom proteomics of four Sidewinder Rattlesnake (Crotalus cerastes) lineages reveal little differential expression despite individual variation. Sci Rep. 2018; 8(1):15534.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  38. Durban J, Juárez P, Angulo Y, Lomonte B, Flores-Díaz M, Alape-Girón A, Sasa M, Sanz L, Gutiérrez JM, Dopazo J, Conesa A, Calvete JJ. Profiling the venom gland transcriptomes of Costa Rican snakes by 454 pyrosequencing. BMC Genomics. 2011; 12:259.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Glenn JL, Straight RC, Wolfe MC, Hardy DL. Geographical variation in Crotalus scutulatus scutulatus (Mojave rattlesnake) venom properties. Toxicon. 1983; 21(1):119–30.

    Article  CAS  PubMed  Google Scholar 

  40. Angulo Y, Chaves E, Alape A, Rucavado A, Gutiérrez JM, Lomonte B. Isolation and characterization of a myotoxic phospholipase A2 from the venom of the arboreal snake Bothriechis (Bothrops) schlegelii from Costa Rica. Archives Biochem Biophys. 1997; 339(2):260–6.

    Article  CAS  Google Scholar 

  41. Tsai I-H, Chen Y-H, Wang Y-M, Tu M-C, Tu AT. Purification, sequencing, and phylogenetic analyses of novel Lys-49 phospholipases A2 from the venoms of rattlesnakes and other pit vipers. Archives Biochem Biophys. 2001; 394(2):236–44.

    Article  CAS  Google Scholar 

  42. Emms DM, Kelly S. OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy. Genome Biol. 2015; 16(1):157.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  43. Durban J, Sanz L, Trevisan-Silva D, Neri-Castro E, Alagón A, Calvete JJ. Integrated venomics and venom gland transcriptome analysis of juvenile and adult Mexican Rattlesnakes Crotalus simus, C. tzabcan, and C. culminatus revealed miRNA-modulated ontogenetic shifts. J Proteome Res. 2017; 16(9):3370–90.

    Article  CAS  PubMed  Google Scholar 

  44. Pla D, Sanz L, Whiteley G, Wagstaff SC, Harrison RA, Casewell NR, Calvete JJ. What killed Karl Patterson Schmidt? Combined venom gland transcriptomic, venomic and antivenomic analysis of the South African green tree snake (the boomslang), Dispholidus typus. Biochimica et Biophysica Acta. 2017; 1861(4):814–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Dowell NL, Giorgianni MW, Kassner VA, Selegue JE, Sanchez EE, Carroll SB. The deep origin and recent loss of venom toxin genes in rattlesnakes. Curr Biol. 2016; 26(18):2434–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Dowell NL, Giorgianni MW, Griffin S, Kassner VA, Selegue JE, Sanchez EE, Carroll SB. Extremely divergent haplotypes in two toxin gene complexes encode alternative venom types within rattlesnake species. Curr Biol. 2018; 28(7):1016–26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Sonenberg N, Hinnebusch AG. Regulation of translation initiation in Eukaryotes: mechanisms and biological targets. Cell. 2009; 136(4):731–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Wong ESW, Belov K. Venom evolution through gene duplications. Gene. 2012; 496(1):1–7.

    Article  CAS  PubMed  Google Scholar 

  49. Aird S, da Silva N, Qiu L, Villar-Briones A, Saddi V, Pires de Campos Telles M, Grau M, Mikheyev A. Coralsnake venomics: analyses of venom gland transcriptomes and proteomes of six Brazilian taxa. Toxins. 2017; 9(12):187.

    Article  PubMed Central  CAS  Google Scholar 

  50. Margres MJ, Bigelow AT, Lemmon EM, Lemmon AR, Rokyta DR. Selection to increase expression, not sequence diversity, precedes gene family origin and expansion in rattlesnake venom. Genetics. 2017; 206(3):1569–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Barua A, Mikheyev AS. Many options, few solutions: over 60 million years snakes converged on a few optimal venom formulations. Mole Biol Evol. 2019; 36(9):1964–74.

    Article  Google Scholar 

  52. Rotenberg D, Bamberger E, Kochva E. Studies on ribonucleic acid synthesis in the venom glands of Vipera palaestinae (ophidia, reptilia). Biochem J. 1971; 121(4):609–12.

    CAS  PubMed  PubMed Central  Google Scholar 

  53. Rokyta DR, Wray KP, Lemmon AR, Lemmon EM, Caudle SB. A high-throughput venom-gland transcriptome for the Eastern Diamondback Rattlesnake (Crotalus adamanteus) and evidence for pervasive positive selection across toxin classes. Toxicon. 2011; 57(5):657–71.

    Article  CAS  PubMed  Google Scholar 

  54. Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc.

  55. Marçais G, Kingsford C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics. 2011; 27(6):764–70.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  56. Mapleson D, Garcia Accinelli G, Kettleborough G, Wright J, Clavijo BJ. KAT: a K-mer analysis toolkit to quality control NGS datasets and genome assemblies. Bioinformatics. 2016; 33(4):663.

    Article  CAS  Google Scholar 

  57. Krueger F. Trim Galore!: A wrapper tool around Cutadapt and FastQC to consistently apply quality and adapter trimming to FastQ files. 2015. http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/.

  58. Zhang J, Kobert K, Flouri T, Stamatakis A. PEAR: a fast and accurate Illumina Paired-End reAd mergeR. Bioinformatics. 2014; 30(5):614–20.

    Article  CAS  PubMed  Google Scholar 

  59. Holding M, Margres M, Mason A, Parkinson C, Rokyta D, Holding ML, Margres MJ, Mason AJ, Parkinson CL, Rokyta DR. Evaluating the performance of de novo assembly methods for venom-gland transcriptomics. Toxins. 2018; 10(6):249.

    Article  PubMed Central  CAS  Google Scholar 

  60. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011; 29(7):644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C, Thierer T, Ashton B, Meintjes P, Drummond A. Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012; 28(12):1647–9.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009; 25(14):1754–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011; 12(1):323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Rokyta DR, Ward MJ. Venom-gland transcriptomics and venom proteomics of the black-back scorpion (Hadrurus spadix) reveal detectability challenges and an unexplored realm of animal toxin diversity. Toxicon. 2017; 128:23–37.

    Article  CAS  PubMed  Google Scholar 

  65. Aitchison J. The statistical analysis of compositional data. J Royal Stat Soc Ser B. 1982; 44(2):139–77.

    Google Scholar 

  66. Palarea-Albaladejo J, Martin-Fernandez JA. zCompositions – R package for multivariate imputation of left-censored data under a compositional approach. Chemometrics Intell Lab Syst. 2015; 143:85–96.

    Article  CAS  Google Scholar 

  67. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014; 15(12):550.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  68. Russo PST, Ferreira GR, Cardozo LE, Bürger MC, Arias-Carrasco R, Maruyama SR, Hirata TDC, Lima DS, Passos FM, Fukutani KF, Lever M, Silva JS, Maracaja-Coutinho V, Nakaya HI. CEMiTool: a Bioconductor package for performing comprehensive modular co-expression analyses. BMC Bioinformatics. 2018; 19(1):56.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  69. Katoh K, Standley DM. MAFFT Multiple Sequence Alignment Software Version 7: improvements in performance and usability. Mol Biol Evol. 2013; 30(4):772–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B. PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol Biol Evol. 2016; 34(3):772–3.

    Google Scholar 

  71. Ronquist F, Teslenko M, Van Der Mark P, Ayres DL, Darling A, Höhna S, Larget B, Liu L, Suchard Ma, Huelsenbeck JP. Mrbayes 3.2: Efficient bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012; 61(3):539–42.

    Article  PubMed  PubMed Central  Google Scholar 

  72. Löytynoja A, Goldman N. Phylogeny-aware gap placement prevents errors in sequence alignment and evolutionary analysis. Science. 2008; 320(5883):1632–5.

    Article  PubMed  CAS  Google Scholar 

  73. Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007; 24(8):1586–91.

    Article  CAS  PubMed  Google Scholar 

  74. Campbell JA, Lamar WW. Venomous Reptiles of the Western Hemisphere, 2nd ed. Ithaca, NY: Cornell University Press; 2004.

    Google Scholar 

Download references

Acknowledgements

We thank Fabian Bonilla and the Chacón and Sandi Harmon families for assistance collecting specimens used in this work.

Funding

This work was funded by National Science Foundation grants DUE 1161228, DEB 1638879, and DEB 1822417 to C.L.P. and DEB 1638902 to D.R.R. as well as the Southwestern Association of Naturalists McCarley Research Grant, the Theodore Roosevelt Memorial Fund through the American Museum of Natural History, and The Explorers Club Exploration Fund and Mamont Scholars program to A.J.M. The funding bodies played no role in study design; data collection, analysis, or interpretation; or writing the manuscript.

Author information

Authors and Affiliations

  1. Department of Biological Sciences, Clemson University, 190 Collings St., Clemson, SC, 29634, USA

    Andrew J. Mason, Mark J. Margres, Jason L. Strickland & Christopher L. Parkinson

  2. Department of Biological Sciences, Florida State University, Tallahassee, FL, 24105, USA

    Darin R. Rokyta

  3. Instituto Clodomiro Picado, Facultad de Microbiologia, Universidad de Costa Rica, San Jose, Costa Rica

    Mahmood Sasa

  4. Department of Forestry, and Environmental Conservation, Clemson University, Clemson, SC, USA

    Christopher L. Parkinson

Authors
  1. Andrew J. Mason
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Mark J. Margres
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Jason L. Strickland
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Darin R. Rokyta
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Mahmood Sasa
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Christopher L. Parkinson
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

A.J.M, M.S., D.R.R., and C.L.P. conceived the study. A.J.M and M.S. collected samples used here. A.J.M. processed data, performed analyses, and drafted the manuscript. M.J.M., J.L.S. provided analytical and conceptual input and contributed to drafting the manuscript. A.J.M., M.J.M, J.L.S, D.R.R., M.S., and C.L.P. reviewed and edited the manuscript. All authors reviewed and approved the final manuscript.

Corresponding author

Correspondence to Christopher L. Parkinson.

Ethics declarations

Ethics approval and consent to participate

The project and its protocols were submitted to and approved by University of Central Florida Institutional Animal Care and Use Committee (IACUC) protocol 16-17W and Clemson University IACUC protocol number 2017-067. Project protocols were also submitted to and approved locally in Costa Rica by the Universidad de Costa Rica’s Animal Care Office (Comité Institucional para el Cuidado y Uso de los Animales) permit CICUA-082-17. Samples were collected and exported from Costa Rica under Sistema Nacional de áreas conservación (SINAC) Investigation Permit 040-2016-INV-ACAT. Permission to collect samples used in this study was received from SINAC and property owners of sites from which animals were collected.

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1

This file contains supplemental table S1 and supplemental figure S1. Supplemental table S1 denotes module assignment of all coding sequence used in WGCNA analyses. Supplemental figure S1 shows phylogenetic trees of orthogroups containing toxins derived from amino acid alignments of coding sequences by OrthoFinder.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Mason, A.J., Margres, M.J., Strickland, J.L. et al. Trait differentiation and modular toxin expression in palm-pitvipers. BMC Genomics 21, 147 (2020). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-020-6545-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-020-6545-9

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us