Trait differentiation and modular toxin expression in palm-pitvipers

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 A₂ (PLA₂) 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 PLA₂ 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.

Keywords: Bothriechis; Gene family evolution; Modularity; Transcriptomics; Venom.

Fig. 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

Fig. 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

Fig. 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 PLA₂ and BPPs

Fig. 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

Fig. 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 PLA₂ 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

Fig. 6

Toxin family phylogenies and expression plots of a C-type lectins (CTLs), b phospholipase A₂s (PLA₂s), 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)

Fig. 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

Fig. 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