Next Article in Journal
Dramatically Enhancing the Sensitivity of Immunoassay for Ochratoxin A Detection by Cascade-Amplifying Enzyme Loading
Previous Article in Journal
Role of Sesamia nonagrioides and Ostrinia nubilalis as Vectors of Fusarium spp. and Contribution of Corn Borer-Resistant Bt Maize to Mycotoxin Reduction
Previous Article in Special Issue
Widespread Evolution of Molecular Resistance to Snake Venom α-Neurotoxins in Vertebrates
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Toxins
Volume 13
Issue 11
10.3390/toxins13110782
toxins-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Varying Intensities of Introgression Obscure Incipient Venom-Associated Speciation in the Timber Rattlesnake (Crotalus horridus)

by
Mark J. Margres
1,*,
Kenneth P. Wray
2,
Dragana Sanader
2,
Preston J. McDonald
1,
Lauren M. Trumbull
1,
Austin H. Patton
3 and
Darin R. Rokyta
2
1
Department of Integrative Biology, University of South Florida, Tampa, FL 33620, USA
2
Department of Biological Science, Florida State University, Tallahassee, FL 32306, USA
3
Department of Integrative Biology, University of California, Berkeley, CA 94720, USA
*
Author to whom correspondence should be addressed.
Toxins 2021, 13(11), 782; https://doi.org/10.3390/toxins13110782
Submission received: 23 September 2021 / Revised: 14 October 2021 / Accepted: 25 October 2021 / Published: 5 November 2021
(This article belongs to the Special Issue Using Genomics to Understand Venom Evolution)
Browse Figures
Review Reports Versions Notes

Abstract

:
Ecologically divergent selection can lead to the evolution of reproductive isolation through the process of ecological speciation, but the balance of responsible evolutionary forces is often obscured by an inadequate assessment of demographic history and the genetics of traits under selection. Snake venoms have emerged as a system for studying the genetic basis of adaptation because of their genetic tractability and contributions to fitness, and speciation in venomous snakes can be associated with ecological diversification such as dietary shifts and corresponding venom changes. Here, we explored the neurotoxic (type A)–hemotoxic (type B) venom dichotomy and the potential for ecological speciation among Timber Rattlesnake (Crotalus horridus) populations. Previous work identified the genetic basis of this phenotypic difference, enabling us to characterize the roles geography, history, ecology, selection, and chance play in determining when and why new species emerge or are absorbed. We identified significant genetic, proteomic, morphological, and ecological/environmental differences at smaller spatial scales, suggestive of incipient ecological speciation between type A and type B C. horridus. Range-wide analyses, however, rejected the reciprocal monophyly of venom type, indicative of varying intensities of introgression and a lack of reproductive isolation across the range. Given that we have now established the phenotypic distributions and ecological niche models of type A and B populations, genome-wide data are needed and capable of determining whether type A and type B C. horridus represent distinct, reproductively isolated lineages due to incipient ecological speciation or differentiated populations within a single species.
Keywords:
adaptation; venom; genome; ecological speciation
Key Contribution: Genetic, proteomic, morphological, and environmental data show that incipient ecological speciation and putative reproductive isolation among Timber Rattlesnakes (Crotalus horridus) with different venom phenotypes vary substantially across the range.

1. Introduction

Natural selection can promote local adaptation and determine the geographical distributions of phenotypes within species [1,2]. Ecological speciation is the process by which populations become reproductively isolated due to divergent, ecologically based natural selection [3], and these divergent selective pressures can be due to differences in the abiotic environment [4], assemblages of interacting species [5], or other ecological factors. For example, a reduced immigrant or hybrid fitness mediated through ecological factors could reduce gene flow between populations [6,7], and reproductive isolation would evolve as a result. Evidence in support of ecological speciation is widespread from diverse taxa [3,8,9], but within species, divergence can occur without leading to speciation [10]. Selection can also reduce the extent of reproductive isolation by favoring introgression of mutually beneficial alleles across population and species boundaries (i.e., increased hybrid fitness via heterozygote advantage) [11,12]. Geography, history, ecology, genetics, selection, and chance determine, through unknown relative contributions, whether new species emerge or are absorbed, but the balance of evolutionary forces responsible is generally obscured by the inadequate assessment of demographic history and the genetics of traits under selection [13].
Snake venoms have emerged as a system for studying adaptive diversification because of their genetic tractability [14,15], contributions to fitness [16,17,18], and high evolutionary rates [19,20,21]. Speciation in venomous snakes is often associated with ecological diversification such as dietary shifts and corresponding changes in venom expression [22,23], and gene expression changes have been implicated in ecological speciation in other systems [24]. Rattlesnake venoms are broadly classified into two categories [25]: neurotoxic, which we refer to as type A, and hemorrhagic, which we refer to as type B [26,27]. This dichotomy reflects an inverse relationship between toxicity and tissue-damaging effects [28]. Type A venoms have high toxicity resulting from the presence of a set of homologous potent presynaptic heterodimeric phospholipase A 2 (PLA 2 ) neurotoxins (e.g., crotoxin, Mojave toxin, and canebrake toxin) and low levels of tissue-damaging snake-venom metalloproteinases (SVMPs). Type B venoms completely lack heterodimeric PLA 2 s and express high levels of SVMPs. Dowell et al. [29,30] showed that different sets of PLA 2 genes, rather than the differential regulation of shared loci, explained much of the differentiation between A and B venoms (but not all; see [15]).
Most of the >30 rattlesnake species have type B venoms, but type A venoms are known from at least ten species [25]. Nine of these ten species, however, are polymorphic for both venom types. The consistency of the A–B venom dichotomy suggests that intermediate venoms are disfavored by selection (although introgression, albeit rare, has been documented in some cases [31,32,33]), providing a mechanism for intrinsic or extrinsic postzygotic reproductive isolation. Although the selective pressures favoring one venom type over the other are unknown, venoms are costly to produce [34]; selection often favors reduced venom expenditures [18,35], suggesting that type A venoms might be generally advantageous due to their heightened toxicity. The predigestive benefits of type B venoms and SVMPs, however, may allow for the consumption of larger prey or digestion at lower temperatures [25,28], potentially expanding the range, altitude, and diet of type B species and populations. The distributions of A-B venoms in some species have been shown to be correlated with the annual environmental temperature [33], suggesting abiotic as well as biotic (e.g., prey) sources of selection may determine phenotype distributions in this system.
The Timber Rattlesnake (Crotalus horridus) occurs in the eastern half of North America from southern Ontario, southward to northern Florida, and westward to Texas and Minnesota [36], exhibiting both A and B venoms [26,27,30,37]. We previously examined venom-gland transcriptomic and venom proteomic divergence between one type A C. horridus from Florida (Osceola National Forest, herein referred to as ONF population) and one type B C. horridus from Georgia (Jones Center at Ichuaway, herein referred to as JC population) [37]. We found that complete changes in venom composition, rather than point mutations affecting venom gene coding sequences, explained much of the differences between venom types, consistent with the work of Dowell et al. [30]. Dowell et al. [30] used BAC cloning to show that type A C. horridus possessed six PLA 2 venom genes (including the acidic and basic subunits of canebrake toxin, the defining neurotoxin of type A venoms) whereas type B individuals possessed only one functional PLA 2 venom gene, highlighting the role gene loss can play in venom evolution. For SVMPs, however, type B C. horridus possessed 13 SVMP paralogs; type A C. horridus possessed nearly all of the same loci, indicating that, in this region, the differential gene regulation accounts for the divergent phenotypes, consistent with other type A species [15].
Given that the genetic bases of both venom phenotypes have been characterized in C. horridus [30,37], we can now begin to determine the roles geography, history, ecology, selection, and chance play in determining when and why new species emerge or are absorbed. Here, we aimed to (1) establish the phenotypic distributions of A and B venoms, (2) determine whether morphological characters were also associated with feeding such as head shape and fang length [33,38,39] vary by venom type, (3) identify the environmental differences between regions harboring different phenotypes, and (4) determine to what extent the venom dichotomy impacts gene flow. To accomplish these goals, we first sampled individuals from adjacent populations in northern Florida (ONF; type A transcriptome locality) and southern Georgia (JC; type B transcriptome locality) [37] and used proteomic, genetic, and morphological data to extensively characterize microgeographic A–B divergence (i.e., populations were ∼225 km apart). We then sampled C. horridus across their entire range to determine the spatial assortment of these phenotypes, explore abiotic characters that may explain these distributions, and test for reciprocal monophyly by phenotype. Overall, we used a combination of data to determine whether type A and type B C. horridus may represent distinct, reproductively isolated lineages due to incipient ecological speciation or differentiated populations within a single species. We recognize that mitochondrial-based inferences about species divergence and boundaries can be problematic (see [40]) and consider this study a first-step in exploring the evolutionary relationships between type A and type B C. horridus.

2. Materials and Methods

2.1. Sampling

We collected and/or were provided venom and blood samples from 22 C. horridus across Florida ( n = 7 ), Georgia ( n = 7 ), Texas ( n = 2 ), Oklahoma, ( n = 1 ), Alabama ( n = 1 ), and Ohio ( n = 4 ). Snout–vent length (SVL) and total length (TL) were recorded for each individual prior to release. Additional tissue was obtained from 236 preserved specimens that were collected by the authors and/or donated to us for the study (see Acknowledgments section). Samples collected by the authors were conducted under permits LSSC-13-00004, LSSC-09-0399, and LSSC-20-00037 from the Florida Fish and Wildlife Conservation Commission (FWC). All procedures were approved by the Florida State University Institutional Animal Care and Use Committee (IACUC) under protocols #0924 and #1333 and the University of South Florida IACUC under protocol #IS00008815.

2.2. Reversed-Phase High-Performance Liquid Chromatography

Reversed-phase high-performance liquid chromatography (RP-HPLC) was performed on a Beckman System Gold HPLC (Beckman Coulter, Fullerton, CA, USA) as previously described [41,42] for the 22 C. horridus venom samples described above. In total, 50 μ g of venom protein was injected onto a Jupiter C18 column (250 × 4.6 mm; Phenomenex, Torrence, CA, USA) using a solvent system of 0.1% trifluoroacetic acid (TFA) in water (solvent A) and 0.075% TFA in acetonitrile (solvent B). Following five minutes at 5% B, a 1% per minute linear gradient of A and B was run to 25% B; this was followed by a 0.25% per minute gradient from 25% to 65% B at a flow rate of 1 mL per minute. Column effluent was monitored at 220 and 280 nm.

2.3. Canebrake Toxin PCR Assay

Genomic DNA was extracted from whole blood samples drawn from the caudal vein or muscle/liver/scale tissue from 258 specimens using the Omega bio-tek E.Z.N.A Tissue DNA Kit according to the manufacturer’s protocol; scale extraction procedures were modified (e.g., increased tissue amounts and addition of 1 M diothiothreitol) as previously described [43]. All extractions were visualized on a 2% agarose gel and run for 30–60 min at 110 V to ensure the presence of high-quality DNA. The acidic and basic subunits of canebrake toxin were amplified from template DNA in 25–50 μL PCR reactions using the acidic subunit-sense (5 ′ GGT ATT TCG TAC TAC AGC TCT TAC GGA 3 ′ ), acidic subunit-antisense (5 ′ TGA TTC CCC CTG GCA ATT 3 ′ ), basic subunit-sense (5 ′ AAC GCT ATT CCC TTC TAT GCC TTT TAC 3 ′ ), and basic subunit-antisense (5 ′ CCT GTC GCA CTC ACA AAT CTG TTC C 3 ′ ) primers, respectively, under the following thermal cycling protocol: 95 ∘ C for five minutes, 35 cycles of 95 ∘ C for 30 s, 55 ∘ C for 30 s, and 72 ∘ C for five minutes, followed by 72 ∘ C for ten minutes [37,44]. Evidence for amplification was visualized by means of a 0.7% agarose gel and ethidium bromide imaged on a Bio Rad Gel Doc using Quantity One Version 4.1.1 or SYBR Safe DNA gel stain and photographed on Dual LED Blue/White Light Transilluminator.
Because these subunits have been shown to be present in type A individuals but absent in type B individuals [30], we were able to phenotype individuals from only genetic samples. If an assay for an individual amplified both subunits, that individual was considered type A. If an assay for an individual did not lead to amplification of either subunit, that individual was considered type B. No individuals amplified only a single subunit which, based on the genomic location of these genes [29,30], was expected. Because all DNA extractions were visualized on a 2% agarose gel to ensure the presence of high-quality DNA and reduce our false-negative rate, we were confident in all type B classifications; in some instances, amplification of cytochrome b (see below) was used to confirm the presence of amplifiable DNA for type B individuals. We do note, however, that this assay could not distinguish between type A and type A+B hybrid venoms, which we discuss below.

2.4. Morphological Analysis

We collected data for 5 morphological characters from 33 preserved C. horridus specimens (Table S1) as previously described [38]: (1) SVL, (distance from tip of snout to posterior edge of cloaca), (2) head length (HL, distance from tip of snout to articular-quadrate joint), (3) head width (HW, distance across widest point of head behind the eyes), (4) interfang distance (IF, distance between fangs at maxillae), and (5) fang length (FL, distance from anterior end of maxilla to the tip of fang as folded against roof of mouth). Snout–vent length was measured to the nearest 0.5 cm using a tailor’s tape, and HL, HW, IF, and FL were each measured to the nearest 0.01 mm using 150 mm digital calipers. Venom type was ascribed based on origin from one of the two focal populations (i.e., JC or ONF).
To remove the effect of size prior to statistical analysis, we subtracted each of the four natural-log-transformed head and fang measurements above (2–4) from the natural-log-transformed SVL [45]. Note that for these transformed values, larger transformed values represent smaller raw values and vice versa. We performed Kruskal–Wallis tests on the transformed data for each variable to determine whether type A and type B individuals significantly differed for any character. We then performed a linear discriminant function analysis using the lda function and leave-one-out crossvalidation from the MASS package [46] in R to assess group membership placement probabilities based solely on morphology across venom type. All statistical analyses were conducted using R v. 3.6.1 [47].

2.5. Niche Modeling

To determine if differences in venom type distributions were correlated with differences in abiotic environmental variables, we followed the approach of Strickland et al. [33] and constructed environmental niche models for each venom type independently using MaxEnt 3.4.1 [48] as implemented in dismo 1.3–3 [49]. We included locality information for 86 type A and 172 type B C. horridus and 19 bioclimatic variables as well as elevation at 30 arc second resolution [50]. All predictor variable rasters were cropped to the extent of the presence records and masked with the IUCN species range [51]. To address collinearity, we performed a principal component analysis and removed principal components with eigenvalues <1 [52]. The two highest loading predictors for each of the remaining principal components were included in the model; these variables comprised Minimum Temperature of Coldest Month, Maximum Temperature of Warmest Month, Precipitation of Warmest Quarter, Mean Temperature of Wettest Quarter, Mean Diurnal Range, Temperature Seasonality, Precipitation of Driest Month, Precipitation of Wettest Quarter, and Precipitation of Coldest Quarter. Models were evaluated by AUC with k = 4 (i.e., 25% of the data) crossvalidation. To reduce the effects of sampling bias on the model, presence records occurring in the same grid cell were removed. After this reduction, the Type A model used 51 training and 18 testing records, and the Type B model used 85 training and 29 testing records. Variable contribution was assessed by permutation and jackknife.
To test niche equivalency across venom types, we used the niche identity test [53] implemented in ENMTools 1.0.5 [54] to determine whether the distribution of type A and type B individuals was identical in geographic and environmental space. Here, the test computed three measures of similarity between our two environmental niche models, Schoener’s D, Warren’s I, and Spearman’s rank correlation, and compared these measures across both geographic and environmental space. For each statistic, null distributions were generated using 99 permutation replicates where the occurrence points of type A and type B individuals were randomly reassigned to one of the two models. The test statistic calculated using the original, unpermuted data was then compared with the null distribution to test for significance using a one-tailed test.

2.6. mtDNA Sequencing and Phylogenetic Reconstruction

DNA was extracted and quality checked as described above. For 56 C. horridus, a 1079 bp fragment of cytochrome b was amplified in 25 μL PCR runs using the H16064 and L14910 primers and thermal cycling protocol previously described [55]. PCR products were purified using the QIagen QIAquick PCR Purification Kit. Sequencing was performed on the Applied Biosystems 3730 Genetic Analyzer.
Sequences were aligned and editing using Geneious Prime and adjusted manually. We then constructed two phylogenies using the same approach outlined below: (1) a ONF ( n = 17 ) and JC ( n = 9 ) phylogeny to match our initial venom and morphological comparisons, and (2) a range-wide phylogeny ( n = 56 ). For each, we fit 88 candidate models of nucleotide substitution using jModelTest 2.1 [56] using the CIPRES portal [57]; the best fit model was determined according to the Bayesian Information Criterion. Using BEAST2 [58], we inferred time-calibrated phylogenies. We used a birth–death tree prior as this prior leads to greater accuracy in the estimation of node times when datasets contain a mixture of intra- and inter-specific sampling [59]. For our clock model, we used an uncorrelated lognormal, relaxed molecular clock [60], and estimated, rather than specified, the rate. We specified two fossil calibrations [61] as previously described [38]. Here, we set the age as a mean of 0, SD of 1, and offset these prior distributions using the early boundary of the time period in which the fossils were confirmed. We constrained the outgroup, C. adamanteus (sequences from [38]), to be monophyletic, and specified their earliest time to the Most Recent Common Ancestor (MRCA) as 0.781 Ma, using a fossil sample from the Early Pleistocene (Irvingtonian I NALMA: [62]). We constrained the earliest time to MRCA for C. horridus to 0.126 Ma based on a fossil sample of this species from the Middle Pleistocene (Irvingtonian II NALMA; [62]). We conducted three independent runs for a length of 50 million generations each, sampling every 5000 generations following an initial burn-in of 10,000 generations. We visually assessed convergence and stationarity for all parameters in Tracer [63]. The three independent runs were subsequently combined, resampling from each every 50,000 generations for a final posterior distribution of 3000 samples. All effective sample sizes for these combined runs exceeded 300. We summarized the posterior distribution of trees using TreeAnnotator to obtain a maximum clade consensus tree using median clade heights for node heights. Trees were visualized in FigTree; clades with posterior probabilities >0.70 were annotated.

3. Results and Discussion

3.1. Fixed Venom Differences, Deep Mitochondrial Genetic Divergence, and Distinct Morphologies over Small Spatial Scales

To characterize A–B divergence at the microgeographic scale, we used proteomic, genetic, and morphological data to compare the ONF and JC populations, which were only ∼225 km apart. We used RP-HPLC to show that each population was monomorphic for their respective venom phenotype (Figure 1A), with all individuals from JC possessing type B venoms, and all individuals from ONF possessing type A venoms. These fixed differences in venom phenotype were consistent with previous results showing that the two canebrake toxin subunits were present in the genomes of animals from ONF but not in the genomes of animals from JC [37]. Overall, the fixed venom phenotypes between JC and ONF populations suggested little-to-no gene flow, reinforcing the hypothesis that intermediate venoms may be disfavored by selection in this system (but see [31,32,33]) and, therefore, may provide a mechanism for intrinsic or extrinsic postzygotic reproductive isolation.
To determine whether A–B venom differences reflected population structure based on neutral markers at this scale, we constructed a time-calibrated Bayesian phylogeny based on cytochrome b sequence for 17 ONF (type A) and 9 JC (type B) individuals (Figure 1B). The venom phenotypes were reciprocally monophyletic with a median divergence time estimate between type A and type B individuals of ∼1.23 Ma (95% CI = 0.249–2.37 Ma). Type A ONF individuals comprised a strongly supported monophyletic clade (posterior probability = 0.99 ) distinct from that of the type B individuals (for which relationships were less well supported). Based on mtDNA, these populations were genetically distinct and putatively reproductively isolated.
To determine whether venom expression differences were integrated with morphological characters also associated with feeding, we measured HL, HW, IF, FL and SVL for 33 preserved C. horridus (18 JC and 15 ONF individuals) and performed Kruskal–Wallis rank sum tests on each morphological character by venom type as previously described [38]. Type A individuals had significantly narrower heads ( p = 0.023 ) and significantly shorter fangs (p < 0.001) than type B individuals (Figure 2A). Although type A snakes also had shorter heads, the difference was not significant ( p = 0.096 ); we identified no differences in IF ( p = 0.745 ). A linear discriminant function analysis based on these characters alone accurately placed 80.0% of type A snakes and 94.4% of type B snakes (Figure 2B), indicating a strong association between head morphology and venom type, consistent with previous results in rattlesnakes [38]. Coefficients of linear discriminants indicated that fang length ( − 12.859 ) was the most important predictor when discriminating among groups (IF = 8.059 , HW = − 4.217 , HL = 6.278 ). Although venom plays a primary role in prey incapacitation and predation in rattlesnakes, morphological characters such as gape and fang length may also determine dietary breadth and feeding ecology [33,38,39]. Margres et al. [38] previously showed that fang length positively covaried with myotoxin (i.e., crotamine) expression in Crotalus adamanteus and posited that the increased depth of injection for myotoxin-rich venoms may increase predation success. Here, the significant reductions in head size and fang length were inversely correlated with overall toxicity, suggesting that the (1) depth of injection may not matter for type A venoms (or at least may be less important relative to type B venoms) and (2) narrower heads may be a consequence of prey differences and/or a reduction in the volume of venom needed given the increased toxicity of type A venoms. Interestingly, similar work in type A and type B Mojave Rattlesnakes (C. scutulatus) did not find significant differences in HW, HL, IF, or FL between type A and type B individuals, although trends did exist for type A individuals to have longer fangs and potentially wider heads [33], opposite of what was found here. The discordance between A–B C. horridus and C. scutulatus individuals suggests that these morphological differences may be prey specific rather than venom-type specific.
Overall, we identified fixed venom differences, a deep phylogenetic split dating back >1 Ma (although we recognize that mitochondrial gene trees can be misleading when delimiting species (see [40]) and distinct head morphologies across a rather small geographic area. Previous work by ourselves [37] and others [30] showed that many of the major venom differences between these populations, such as the lack of heterodimeric PLA 2 neurotoxins in type B venoms, reflected an absence in the genome of the corresponding genes rather than the differential regulation of shared genes. Together, these results suggest that type A ONF and type B JC C. horridus may represent independently evolving lineages, perhaps due to incipient ecological speciation. Alternatively, these differences could reflect incipient allopatric speciation as these populations occur on either side of the Suwannee River, a known biogeographic barrier for continental and peninsular Florida organisms [64], including rattlesnakes [13,38]. The Suwannee Straits ran through the Okefenokee Trough from the Gulf of Mexico to the Atlantic Ocean, repeatedly separating peninsular Florida from the continent [65]. Although the Straits likely closed ∼4 Ma, some evidence suggests that the Straits may have opened briefly as recently as ∼1.75 Ma [38,64]. With our current dataset, we could not distinguish between incipient allopatric and ecological speciation at this scale. Crotalus horridus, however, has a very large geographic range, and ecological speciation comes with explicit predictions. Isolation by environment (IBE) should be much greater than isolation by distance (IBD) if ecological differences are leading to reproductive isolation [66]. To determine whether phenotypes were spatially distinct or overlapping, we next sampled more broadly across the range to identify the distributions of each phenotype.

3.2. Type A Venom Dominates along the Southern Periphery of the Range

We assayed 258 individuals across the range of C. horridus and detected 86 type A individuals and 172 type B individuals (Figure 3). Type A venoms were largely present in the southern periphery of the range, and type B venoms dominated the northern and interior parts of the range. Type A venoms may be common along the western range edge as well, but sampling in this region was too sparse to make a determination. Both phenotypes co-occurred in the coastal plain regions of South Carolina, Georgia (see inset in Figure 3), and Orleans Parish, Louisiana, although sampling in the latter was sparse (three type Bs and one type A). We note that our assay could not distinguish between type A and type A+B “hybrid” venoms as both venom types would possess the canebrake toxin subunits. Therefore, some of the individuals classified as type A in Figure 3 could represent A+B venoms, and one type A individual included here from Jekyll Island, Georgia was a C. horridus × C. adamanteus hybrid [67].
Venom-based analyses showed that both type A and type B phenotypes were similar across the range (e.g., Florida and Texas type A venoms were very similar; Figure 1 and Figure 4), although there appeared to be more variation within type B venoms based on our limited venom-based sampling (Figure 4). Given (1) the much larger geographic distribution of type B venoms and, therefore, greater variation in abiotic and biotic selective pressures and (2) that type B venoms are much more complex than type A venoms (Figure 1 and Figure 4; [15,23]) and have many more axes on which they can vary, this result was expected.
Overall, venom types were spatially sorted across the range, suggesting an ecological divergence between type A and type B individuals. To test whether differences in venom-type distributions were correlated with differences in abiotic environmental variables, we next constructed environmental niche models for each venom type independently.

3.3. Type A and B Individuals Occupy Significantly Different Niches

The distributions of type A and type B individuals were significantly different in both geographic ( D G e o = 0.29 , I G e o = 0.58 , Spearman’s rank correlation = − 0.03 ; p < 0.01 for all tests) and environmental ( D E n v = 0.19 , I E n v = 0.44 , Spearman’s rank correlation = 0.33 ; p < 0.01 for all tests) space (Figure 5). For the type A model (AUC: training 0.948, test 0.960; Figure 5A), the mean temperature of the wettest quarter (72.2%) and minimum temperature of the coldest month (14.2%) were the most informative environmental variables. Permutation importance (mean temperature of the wettest quarter p e r m = 26.1 %; minimum temperature of the coldest month p e r m = 60.8 %) and jackknife testing confirmed these results, with mean temperature of the wettest quarter being the most informative variable during jackknife testing. For the type B model (AUC: training 0.826, test 0.810; Figure 5B), precipitation of warmest quarter (47.4%) and temperature seasonality (24.1%) were the most informative environmental variables. Results were again confirmed by permutation importance (precipitation of warmest quarter p e r m = 37.3 %; temperature seasonality p e r m = 26.4 %) and jackknife testing, with precipitation of warmest quarter being the most informative variable during jackknife testing.
Temperature appeared to be a key factor in determining the distributions of type A and B venoms in C. horridus. Type A habitat suitability was positively correlated with mean temperature of the wettest quarter and minimum temperature of the coldest month, the two most informative variables in the model, whereas these variables had relatively low contributions in the type B model (<10%) and were negatively associated with type B habitat suitability. Previous work showed that the distributions of A–B venoms in Mojave Rattlesnakes (Crotalus scutulatus) was correlated with minimum temperature of the coldest quarter and mean temperature of wettest quarter [33], consistent with our results here. Given that type A venoms within polymorphic species appear to occur in warmer parts of the range, perhaps the beneficial predigestive benefits of type B venoms and SVMPs are no longer required for effective predation/digestion along the southern periphery [25,28]. In this case, the more toxic type A venoms may be advantageous in the absence of digestive constraints, but when these constraints do apply, type B venoms may be favored, thus limiting the northward expansion of the type A phenotype.
Type A and B venoms, however, could also be locally adapted for different prey species or genetically distinct populations of the same species, and the relationship with temperature may simply correlate with local prey community composition. For example, previous work in Pacific Rattlesnakes (C. oreganus) found that biotic information about local prey species explained more variation in venom expression than abiotic environmental data [68]. These hypotheses, however, are not mutually exclusive; here, type B venoms may be adapted for a certain species of prey and pro-digestive effects whereas type A venoms may be adapted for a different species of prey and also be free of digestive constraint. Although diet information for most pit vipers is sparse, we examined 22 published articles, books, and conference proceedings [69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90] containing dietary information for C. horridus across their range to compare prey items in regions where C. horridus was monomorphic for type B to regions where C. horridus was monomorphic for type A (regions where C. horridus was polymorphic in venom type were excluded). We found 130 type B records but only 10 type A records (Table S2). Mammals (∼76% of type B records, 60% of type A records) and birds (∼17% of type B records, 30% of type A records) comprised the majority of the diet for each venom type. Although the small sample size in type A prevented statistical comparisons, the largest presence/absence difference was found in rodents belonging to the family Cricetidae. Cricetid rodents (e.g., Peromyscus, Microtus, and Sigmodon) were ∼33% of all prey items recorded for type B animals and were notably absent from the type A records. Type B venoms, therefore, may be locally adapted for Cricetid rodents, and perhaps type A venoms are locally adapted to a different prey species and/or abiotic conditions. Toxicity and digestive assays quantifying these differences in natural prey items are needed to test these hypotheses; although interesting and ultimately necessary to characterize the source of ecological divergence (if present), such work was beyond the scope of the present study.
The identified significant differences in niche space, accompanied with differences in multiple, putatively adaptive phenotypes in venom and head morphology, again suggested incipient ecological speciation. Assessing the degree of reproductive isolation is required when investigating putative species boundaries in wide-ranging taxa, especially along contact zones [40,91]. Therefore, we next sequenced cytochrome b across the range, with a focus on the contact zone in the southeastern coastal plain, to determine whether venom types retained reciprocal monophyly (as found in the JC and ONF populations) or if reproductive isolation between venom types broke down at larger geographic scales.

3.4. A Lack of Reciprocal Monophyly Suggests Gene Flow between Venom Types across the Range

To determine whether the A–B reciprocal monophyly detected for the ONF and JC populations (Figure 1B) was found range wide, we sequenced cytochrome b for 30 type A and 26 type B individuals across the range and constructed a Bayesian phylogeny (Figure 6). Although we identified two well-supported monophyletic clades of type A animals (one including individuals from Florida, Georgia, and South Carolina; one including all three Texas type As), neither venom type was monophyletic, and no well-supported monophyletic clades of only type B animals were identified (Figure 6). This pattern suggested that venom types were not reproductively isolated and is consistent with introgression occurring primarily from type A populations into type B populations. However, more extensive sampling of nuclear genetic data is necessary to test this hypothesis.
Given the lack of reciprocal monophyly based on venom type, we next looked at the spatial distribution of the population structure identified in the phylogeny (Figure 6). We identified six clades with varying levels of support: (1) a well-supported type A-only clade containing all 17 ONF individuals, as well as one individual from Georgia and one from South Carolina (pink clade in Figure 6), (2) a well-supported western clade of mostly type B animals (seven type Bs and two type As) along the Mississippi drainage (cyan clade in Figure 6), (3) a well-supported type A-only clade containing three Texas individuals (black clade in Figure 6), (4) a less supported Georgia–South Carolina clade containing all JC individuals and four other individuals (11 type Bs and 2 type As; orange clade in Figure 6), (5) a well-supported South Carolina clade (three type Bs and four type As; gray clade in Figure 6), and (6) a less supported type B only clade in the southeast (yellow clade in Figure 6). Surprisingly, the type A-only pink clade that contained all 17 ONF individuals and 1 individual each from Georgia and South Carolina was more closely related to the western lineages (cyan Mississippi Drainage and black Texas lineages in Figure 6) than the other three lineages identified in the southeast (yellow, orange, and gray lineages in Figure 6); this relationship, however, was not strongly supported (posterior probability = 0.74 ). While this could indicate that the type A phenotype arose in the west and has since spread east, denser sampling across the genome and the range is needed to determine the degree of reproductive isolation between venom types/populations as well as the origin(s) of the type A phenotype in C. horridus. Overall, these results revealed substantial geographic population structure within C. horridus that was only weakly associated with venom type.
How putatively distinct lineages interact along contact zones and potential species boundaries is essential for determining whether such lineages are evolving independently as distinct species or collectively as geographically differentiated populations [40]. Based on a mitochondrial marker, type A and type B C. horridus were not independently evolving lineages and appeared to readily interbreed in coastal Georgia and South Carolina (Figure 6) despite isolation elsewhere (Figure 1). These equivocal results support varying intensities of introgression between type A and type B individuals across the range, and this variance may be the result of variance in ecologically divergent selective pressures. Here, ecologically divergent selection would be strong between the ONF and JC populations, thus producing substantial genetic divergence and putative reproductive isolation, whereas selection would be weaker in the coastal plain of South Carolina and Georgia (i.e., orange and gray clades in Figure 6) where the evidence of introgression was strong. Although the ecologically divergent selection pressures may act on the venom phenotype as predicted in this study, we could not rule out the possibility of other loci in linkage disequilibrium with the PLA 2 venom genes being the actual targets of selection, perhaps for a trait not related to predation and feeding ecology. Additionally, neutral processes due to historical biogeographic events (e.g., the formation of the Suwannee Straits as discussed above [65]) could also produce similar patterns. Given that mitochondrial gene trees are often misleading when making inferences about species boundaries [40], nuclear genomic and ecological data across the range with targeted, dense sampling along contact zones is required to better understand the population structure and adequately assess reproductive isolation. Nevertheless, our current mitochondrial data do not support type A and type B C. horridus as distinct lineages range-wide.

3.5. Conclusions

We used a combination of genetic, proteomic, morphological, and environmental data to test whether type A and type B C. horridus may represent distinct, reproductively isolated lineages due to incipient ecological speciation or differentiated populations within a single species. Although we found significant genetic, proteomic, morphological, and environmental differences between type A and type B C. horridus at microgeographic scales (i.e., between ONF and JC populations), the lack of reciprocal monophyly among venom types range-wide suggested varying intensities of gene flow and an overall lack of reproductive isolation across the range. Speciation is a continuous process [92], and speciation with gene flow may be common [93], but resolving species relationships among lineages with complex, reticulate evolutionary histories remains challenging [91]. To address this issue, Marshall et al. [40] recently proposed the use of genome-wide nuclear data and intensive sampling along contact zones to determine whether the lineages of interest are evolving mostly independently (i.e., reproductively isolated with narrow zones of F 1 hybrids) or as geographically differentiated populations (i.e., intergradation across wider zones with F 2 + hybrids). Such an approach here with increased sampling in the coastal plain regions of South Carolina and Georgia as well as along the Gulf Coast is necessary to disentangle the complex evolutionary relationships between type A and type B C. horridus.

Supplementary Materials

The following is available online at https://www.mdpi.com/article/10.3390/toxins13110782/s1, Table S1: morphology data; Table S2: diet data. The data used to construct the niche models are available on request from the corresponding author. The data are not publicly available to protect specific locality information for a species that is protected and imperiled across much of its range.

Author Contributions

M.J.M., D.R.R. and K.P.W. conceived of and designed the study. D.S., P.J.M., L.M.T., K.P.W. and M.J.M. generated the data. M.J.M., K.P.W., P.J.M., A.H.P. and D.R.R. analyzed the data. All authors wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This work was funded by the National Science Foundation (DEB 1145987 to DRR) and the University of South Florida (to MJM).

Institutional Review Board Statement

This study was conducted under permits LSSC-13-00004, LSSC-09-0399, and LSSC-20-00037 from the Florida Fish and Wildlife Conservation Commission (FWC). All procedures were approved by the Florida State University Institutional Animal Care and Use Committee (IACUC) under protocols #0924 and #1333 and the University of South Florida IACUC under protocol #IS00008815.

Informed Consent Statement

Not applicable.

Data Availability Statement

All cytochrome b sequences used in our phylogenetic analyses were deposited on NCBI under accessions MZ394520–MZ394575.

Acknowledgments

We thank Jennifer and Daniel Abney, Richard Bartlett, Joseph Colbert and the Jekyll Island Authority, Brian Crother, Matthew Holding, Mark S. and Robin Margres, and Lora Smith and Jennifer Howze at the Joseph W. Jones Ecological Research Center at Ichuaway for help in acquiring venom and tissue samples. We thank Paul Moler for donating hundreds of tissue samples to the project as well as Max Nickerson and the Florida Museum of Natural History for access to preserved specimens for morphological analysis. We thank Rhett M. Rautsaw and Erich P. Hofmann for assistance in collating diet information.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hendry, A.P.; Nosil, P.; Rieseberg, L.H. The speed of ecological speciation. Funct. Ecol. 2007, 21, 455–464. [Google Scholar] [CrossRef] [Green Version]
  2. Nosil, P.; Funk, D.J.; Ortiz-Barrientos, D. Divergent selection and heterogeneous genomic divergence. Mol. Ecol. 2009, 18, 375–402. [Google Scholar] [CrossRef]
  3. Schluter, D. Evidence for ecological speciation and its alternative. Science 2009, 323, 737–741. [Google Scholar] [CrossRef] [Green Version]
  4. Keller, I.; Seehausen, O. Thermal adaptation and ecological speciation. Mol. Ecol. 2012, 21, 782–799. [Google Scholar] [CrossRef]
  5. Öhlund, G.; Bodin, M.; Nilsson, K.A.; Öhlund, S.O.; Mobley, K.B.; Hudson, A.G.; Peedu, M.; Brännström, Å.; Bartels, P.; Præbel, K.; et al. Ecological speciation in European whitefish is driven by a large-gaped predator. Evol. Lett. 2020, 4, 243–256. [Google Scholar] [CrossRef] [PubMed]
  6. Garant, D.; Forde, S.E.; Hendry, A.P. The multifarious effects of dispersal and gene flow on contemporary adaptation. Funct. Ecol. 2007, 21, 434–443. [Google Scholar] [CrossRef]
  7. Zhang, L.; Hood, G.R.; Roush, A.M.; Shzu, S.A.; Comerford, M.S.; Ott, J.R.; Egan, S.P. Asymmetric, but opposing reductions in immigrant viability and fecundity promote reproductive isolation among host-associated populations of an insect herbivore. Evolution 2021, 75, 476–489. [Google Scholar] [CrossRef]
  8. Ryan, P.G.; Bloomer, P.; Moloney, C.L.; Grant, T.J.; Delport, W. Ecological speciation in South Atlantic island finches. Science 2007, 315, 1420–1423. [Google Scholar] [CrossRef] [Green Version]
  9. Shafer, A.B.; Wolf, J.B. Widespread evidence for incipient ecological speciation: A meta-analysis of isolation-by-ecology. Ecol. Lett. 2013, 16, 940–950. [Google Scholar] [CrossRef]
  10. Nosil, P.; Harmon, L.J.; Seehausen, O. Ecological explanations for (incomplete) speciation. Trends Ecol. Evol. 2009, 24, 145–156. [Google Scholar] [CrossRef] [PubMed]
  11. Song, Y.; Endepols, S.; Klemann, N.; Richter, D.; Matuschka, F.R.; Shih, C.H.; Nachman, M.W.; Kohn, M.H. Adaptive introgression of anticoagulant rodent poison resistance by hybridization between old world mice. Curr. Biol. 2011, 21, 1296–1301. [Google Scholar] [CrossRef] [Green Version]
  12. Griffiths, J.S.; Kawji, Y.; Kelly, M.W. An experimental test of adaptive introgression in locally adapted populations of splash pool copepods. Mol. Biol. Evol. 2020, 38, 1306–1316. [Google Scholar] [CrossRef] [PubMed]
  13. Margres, M.; Patton, A.; Wray, K.; Hassinger, A.; Ward, M.; Lemmon, E.; Lemmon, A.; Rokyta, D. Tipping the scales: The migration-selection balance leans toward selection in snake venoms. Mol. Biol. Evol. 2019, 36, 271–282. [Google Scholar] [CrossRef]
  14. Rokyta, D.; Margres, M.; Calvin, K. Post-transcriptional mechanisms contribute little to phenotypic variation in snake venoms. G3 Genes Genomes Genet. 2015, g3, 115. [Google Scholar] [CrossRef] [Green Version]
  15. Margres, M.J.; Rautsaw, R.M.; Strickland, J.L.; Mason, A.J.; Schramer, T.D.; Hofmann, E.P.; Stiers, E.; Ellsworth, S.A.; Nystrom, G.S.; Hogan, M.P.; et al. The Tiger Rattlesnake genome reveals a complex genotype underlying a simple venom phenotype. Proc. Natl. Acad. Sci. USA 2021, 118, e2014634118. [Google Scholar] [CrossRef] [PubMed]
  16. Daltry, J.C.; Wüster, W.; Thorpe, R.S. Diet and snake venom evolution. Nature 1996, 379, 537–540. [Google Scholar] [CrossRef]
  17. Holding, M.; Biardi, J.; Gibbs, H. Coevolution of venom function and venom resistance in a rattlesnake predator and its squirrel prey. Proc. R. Soc. B 2016, 283, 20152841. [Google Scholar] [CrossRef]
  18. Margres, M.; Wray, K.; Hassinger, A.; Ward, M.; McGivern, J.; Lemmon, E.; Lemmon, A.; Rokyta, D. Quantity, not quality: Rapid adaptation in a polygenic trait proceeded exclusively through expression differentiation. Mol. Biol. Evol. 2017, 34, 3099–3110. [Google Scholar] [CrossRef] [Green Version]
  19. Lynch, V.J. Inventing an arsenal: adaptive evolution and neofunctionalization of snake venom phospholipase A2 genes. BMC Evol. Biol. 2007, 7, 2. [Google Scholar] [CrossRef] [Green Version]
  20. Margres, M.; Wray, K.; Seavy, M.; McGivern, J.; Herrera, N.; Rokyta, D. Expression differentiation is constrained to low-expression proteins over ecological timescales. Genetics 2016, 202, 273–283. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  21. Margres, M.; Bigelow, A.; Lemmon, E.; Lemmon, A.; Rokyta, D. Selection to increase expression, not sequence diversity, precedes gene family origin and expansion in rattlesnake venom. Genetics 2017, 206, 1569–1580. [Google Scholar] [CrossRef] [Green Version]
  22. Gibbs, H.L.; Rossiter, W. Rapid Evolution by Positive Selection and Gene Gain and Loss: PLA2 Venom Genes in Closely Related Sistrurus Rattlesnakes with Divergent Diets. J. Mol. Evol. 2008, 66, 151–166. [Google Scholar] [CrossRef] [PubMed]
  23. Holding, M.L.; Strickland, J.L.; Rautsaw, R.M.; Hofmann, E.P.; Mason, A.J.; Hogan, M.P.; Nystrom, G.S.; Ellsworth, S.A.; Colston, T.J.; Borja, M.; et al. Phylogenetically diverse diets favor more complex venoms in North American pitvipers. Proc. Natl. Acad. Sci. USA 2021, 118, e2015579118. [Google Scholar] [CrossRef]
  24. Pavey, S.; Collin, H.; Nosil, P.; Rogers, S. The role of gene expression in ecological speciation. Ann. N. Y. Acad. Sci. 2010, 1206, 110–129. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Mackessy, S.P. Venom composition in rattlesnakes: Trends and biological significance. In The Biology of Rattlesnakes; Hayes, W.K., Beaman, K.R., Cardwell, M.D., Bush, S.P., Eds.; Loma Linda University Press: Loma Linda, CA, USA, 2008; pp. 495–510. [Google Scholar]
  26. Straight, R.C.; Glenn, J.L.; Wolt, T.B.; Wolfe, M.C. Regional differences in content of small basic peptide toxins in the venoms of Crotalus adamanteus and Crotalus horridus. Comp. Biochem. Physiol. 1991, 100, 51–58. [Google Scholar] [CrossRef]
  27. Glenn, J.L.; Straight, R.C.; Wolf, T.B. Regional variation in the presence of canebrake toxin in Crotalus horridus venom. Comp. Biochem. Phys. C 1994, 107, 337–346. [Google Scholar] [CrossRef]
  28. Mackessy, S.P. Evolutionary trends in venom composition in the Western Rattlesnakes (Crotalus viridis sensu lato): Toxicity vs. tenderizers. Toxicon 2010, 55, 1463–1474. [Google Scholar] [CrossRef]
  29. Dowell, N.; Giorgianni, M.; Kassner, V.; Selegue, J.; Sanchez, E.; Carroll, S. The deep origin and recent loss of venom toxin genes in rattlesnakes. Curr. Biol. 2016, 26, 2434–2445. [Google Scholar] [CrossRef] [Green Version]
  30. Dowell, N.; Giorgianni, M.; Griffin, A.; Kassner, V.; Selegue, J.; Sanchez, E.; Carroll, S. Extremely divergent haplotypes in two toxin gene complexes encode alternative venom types within rattlesnake species. Curr. Biol. 2018, 28, 1016–1026. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  31. Wilkinson, J.A.; Glenn, J.L.; Straight, R.C.; Sites, J.W., Jr. Distribution and genetic variation in venom A and B populations of the Mojave rattlesnake Crotalus scutulatus scutulatus in Arizona. Herpetologica 1991, 47, 54–68. [Google Scholar]
  32. Zancolli, G.; Baker, T.G.; Barlow, A.; Bradley, R.K.; Calvete, J.J.; Carter, K.C.; De Jager, K.; Owens, J.B.; Price, J.F.; Sanz, L.; et al. Is hybridization a source of adaptive venom variation in rattlesnakes? A test, using a Crotalus scutulatus × viridis hybrid zone in southwestern New Mexico. Toxins 2016, 8, 188. [Google Scholar] [CrossRef] [Green Version]
  33. Strickland, J.; Smith, C.; Mason, A.; Schield, D.; Borja, M.; Cataneda-Gaytan, G.; Spencer, C.; Smith, L.; Trapaga, A.; Bouzid, N.; et al. Evidence for divergent patterns of local selection driving venom variation in Mojave Rattlesnakes (Crotalus scutulatus). Sci. Rep. 2018, 8, 17622. [Google Scholar] [CrossRef]
  34. McCue, M.D. Cost of Producing Venom in Three North American Pitviper Species. Copeia 2006, 2006, 818–825. [Google Scholar] [CrossRef]
  35. Barlow, A.; Pook, C.E.; Harrison, R.A.; Wüster, W. Coevolution of diet and prey-specific venom activity supports the role of selection in snake venom evolution. Proc. R. Soc. Lond. B 2009, 276, 2443–2449. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. McDiarmid, R.W.; Campbell, J.A.; Touré, T. Snake Species of the World: A Taxonomic and Geographic Reference, Vol. 1; Herpetologists’ League: Washington, DC, USA, 1999. [Google Scholar]
  37. Rokyta, D.; Wray, K.; McGivern, J.; Margres, M. The transcriptomic and proteomic basis for the evolution of a novel venom phenotype within the Timber Rattlesnake (Crotalus horridus). Toxicon 2015, 98, 34–48. [Google Scholar] [CrossRef] [PubMed]
  38. Margres, M.; Wray, K.; Seavy, M.; McGivern, J.; Sanader, D.; Rokyta, D. Phenotypic integration in the feeding system of the eastern diamondback rattlesnake (Crotalus adamanteus). Mol. Ecol. 2015, 24, 3405–3420. [Google Scholar] [CrossRef]
  39. Cleuren, S.G.; Hocking, D.P.; Evans, A.R. Fang evolution in venomous snakes: Adaptation of 3D tooth shape to the biomechanical properties of their prey. Evolution 2021, 75, 1377–1394. [Google Scholar] [CrossRef]
  40. Marshall, T.L.; Chambers, E.A.; Matz, M.V.; Hillis, D.M. How mitonuclear discordance and geographic variation have confounded species boundaries in a widely studied snake. Mol. Phylogenet. Evol. 2021, 162, 107194. [Google Scholar] [CrossRef]
  41. Margres, M.; McGivern, J.; Wray, K.; Seavy, M.; Calvin, K.; Rokyta, D. Linking the transcriptome and proteome to characterize the venom of the eastern diamondback rattlesnake (Crotalus adamanteus). J. Proteom. 2014, 96, 145–158. [Google Scholar] [CrossRef]
  42. Margres, M.; McGivern, J.; Seavy, M.; Wray, K.; Facente, J.; Rokyta, D. Contrasting modes and tempos of venom expression evolution in two snake species. Genetics 2015, 199, 165–176. [Google Scholar] [CrossRef] [Green Version]
  43. Campos, P.; Gilbert, T. DNA Extraction from Keratin and Chitin. Methods Mol. Biol. (Clifton N.J.) 2012, 840, 43–49. [Google Scholar] [CrossRef]
  44. Wooldridge, B.J.; Pineda, G.; Banuelas-Ornelas, J.J.; Dagda, R.K.; Gasanov, S.E.; Rael, E.D.; Lieb, C.S. Mojave rattlesnakes (Crotalus scutulatus scutulatus) lacking the acidic subunit DNA sequence lack Mojave toxin in their venom. Comp. Biochem. Phys. B 2001, 130, 169–179. [Google Scholar] [CrossRef]
  45. Lombarte, A.; Lleonart, J. Otolith size changes related with body growth, habitat depth and temperature. Environ. Biol. Fish. 1993, 37, 297–306. [Google Scholar] [CrossRef]
  46. Venables, W.N.; Ripley, B.D. Modern Applied Statistics with S, 4th ed.; Springer: New York, NY, USA, 2002; ISBN 0-387-95457-0. [Google Scholar]
  47. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2019. [Google Scholar]
  48. Phillips, S.J.; Anderson, R.P.; Schapire, R.E. Maximum entropy modeling of species geographic distributions. Ecol. Model. 2006, 190, 231–259. [Google Scholar] [CrossRef] [Green Version]
  49. Hijmans, R.J.; Phillips, S.; Leathwick, J.; Elith, J.; Hijmans, M.R.J. Package ‘dismo’. Circles 2017, 9, 1–68. [Google Scholar]
  50. Fick, S.E.; Hijmans, R.J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 2017, 37, 4302–4315. [Google Scholar] [CrossRef]
  51. Hammerson, G. Crotalus horridus. The IUCN Red List of Threatened Species 2007: e.T64318A12765920. 2007. Available online: http://dx.doi.org/10.2305/IUCN.UK.2007.RLTS.T64318A12765920.en (accessed on 22 September 2021).
  52. Dormann, C.F.; Elith, J.; Bacher, S.; Buchmann, C.; Carl, G.; Carré, G.; Marquéz, J.R.G.; Gruber, B.; Lafourcade, B.; Leitão, P.J.; et al. Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography 2013, 36, 27–46. [Google Scholar] [CrossRef]
  53. Warren, D.L.; Glor, R.E.; Turelli, M. Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution 2008, 62, 2868–2883. [Google Scholar] [CrossRef] [PubMed]
  54. Warren, D.L.; Matzke, N.J.; Cardillo, M.; Baumgartner, J.B.; Beaumont, L.J.; Turelli, M.; Glor, R.E.; Huron, N.A.; Simões, M.; Iglesias, T.L.; et al. ENMTools 1.0: An R package for comparative ecological biogeography. Ecography 2021, 44, 504–511. [Google Scholar] [CrossRef]
  55. Burbrink, F.T.; Lawson, R.; Slowinski, J. Mitochondrial DNA phylogeography of the polytypic North American rat snake (Elaphe obsoleta): A critique of the subspecies concept. Evolution 2000, 54, 2107–2118. [Google Scholar] [CrossRef]
  56. Darriba, D.; Taboada, G.; Doallo, R.; Posada, D. jModelTest2: More models, new heuristics and parallel computing. Nat. Methods 2012, 9, 772. [Google Scholar] [CrossRef] [Green Version]
  57. Miller, M.A.; Pfeiffer, W.; Schwartz, T. The CIPRES science gateway: Enabling high-impact science for phylogenetics researchers with limited resources. In Proceedings of the 1st Conference of the Extreme Science and Engineering Discovery Environment: Bridging from the Extreme to the Campus and Beyond, Chicago, IL, USA, 16–20 July 2012; pp. 1–8. [Google Scholar]
  58. Bouckaert, R.; Vaughan, T.G.; Barido-Sottani, J.; Duchêne, S.; Fourment, M.; Gavryushkina, A.; Heled, J.; Jones, G.; Kühnert, D.; De Maio, N.; et al. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 2019, 15, e1006650. [Google Scholar] [CrossRef] [Green Version]
  59. Ritchie, A.M.; Lo, N.; Ho, S.Y. The impact of the tree prior on molecular dating of data sets containing a mixture of inter-and intraspecies sampling. Syst. Biol. 2017, 66, 413–425. [Google Scholar] [CrossRef] [PubMed]
  60. Drummond, A.; Ho, S.; Phillips, M.; Rambaut, A. Relaxed phylogenetics and dating with confidence. PLoS Biol. 2006, 4, e88. [Google Scholar] [CrossRef]
  61. Ho, S. Calibrating molecular estimates of substitution rates and divergence times in birds. J. Avian Biol. 2007, 38, 409–414. [Google Scholar] [CrossRef]
  62. Holman, J. Fossil Snakes of North America: Origin, Evolution, Distribution, Paleoecology; Indiana University Press: Bloomington, IN, USA, 2000. [Google Scholar]
  63. Rambaut, A.; Drummond, A. Tracer v1.4. 2007. Available online: http://beast.bio.ed.ac.uk (accessed on 22 September 2021).
  64. Bert, T. Speciation in western Atlantic stone crabs (genus Menippe): The role of geological processes and climatic events in the formation and distribution of species. Mar. Biol. 1986, 93, 157–170. [Google Scholar] [CrossRef]
  65. Randazzo, A.; Jones, D. (Eds.) The Geology of Florida; University Press of Florida: Gainesville, FL, USA, 1997. [Google Scholar]
  66. Wang, I.J.; Bradburd, G.S. Isolation by environment. Mol. Ecol. 2014, 23, 5649–5662. [Google Scholar] [CrossRef]
  67. Harrison, C.M.; Colbert, J.; Richter, C.J.; McDonald, P.J.; Trumbull, L.M.; Ellsworth, S.A.; Hogan, M.P.; Rokyta, D.R.; Margres, M.J. Using morphological, genetic, and venomic analyses to present current and historic evidence of Crotalus horridus × Crotalus adamanteus hybridization on Jekyll Island GA. Southeast. Nat. 2021, in press. [Google Scholar]
  68. Holding, M.L.; Margres, M.J.; Rokyta, D.R.; Gibbs, H.L. Local prey community composition and genetic distance predict venom divergence among populations of the Northern Pacific rattlesnake (Crotalus oreganus). J. Evol. Biol. 2018, 31, 1513–1528. [Google Scholar] [CrossRef]
  69. Babcock, H.L. Food habits of the timber rattlesnake. Bull. Boston Soc. Nat. Hist. 1929, 51, 12–14. [Google Scholar]
  70. King, W. A Survey of the Herpetology of Great Smoky Mountains National Park. Am. Midl. Nat. 1939, 21, 531–582. [Google Scholar] [CrossRef]
  71. Uhler, F.M.; Cottam, C.; Clarke, T.E. Food of snakes of the George Washington National Forest, Virginia. In Proceedings of the Fourth National American Wildlife Conference, Detroit, MI, USA, 13–15 February 1939; pp. 605–662. [Google Scholar]
  72. Smyth, T. Notes on the Timber Rattlesnake at Mountain Lake, Virginia. Copeia 1949, 1949, 78. [Google Scholar] [CrossRef]
  73. Barbour, R.W. The Reptiles of Big Black Mountain, Harlan County, Kentucky. Copeia 1950, 1950, 100–107. [Google Scholar] [CrossRef]
  74. Martin, J.R.; Wood, J.T. Notes on the Poisonous Snakes of the Dismal Swamp Area. Herpetologica 1955, 11, 237–238. [Google Scholar]
  75. Hamilton, W.J.; Pollack, J.A. The Food of Some Colubrid Snakes from Fort-Benning, Georgia. Ecology 1956, 37, 519–526. [Google Scholar] [CrossRef]
  76. Bush, F.M. Foods of Some Kentucky Herptiles. Herpetologica 1959, 15, 73–77. [Google Scholar]
  77. Keenlyne, K.D. Sexual Differences in Feeding Habits of Crotalus horridus horridus. J. Herpetol. 1972, 6, 234–237. [Google Scholar] [CrossRef]
  78. Brown, E.E. Some Snake Food Records From The Carolinas, USA. Brimleyana 1979, 1, 113–124. [Google Scholar]
  79. Fitch, H.S. Resources of a Snake Community in Prairie-Woodland Habitat of Northeastern Kansas. In Herpetological Communities: A Symposium of the Society for the Study of Amphibians and Reptiles and the Herpetologists’ League, August 1977; Scott, N.J., Jr., Ed.; U.S. Fish and Wildlife Service: Washington, DC, USA, 1982; pp. 83–97. [Google Scholar]
  80. Reinert, H.K.; Cundall, D.; Bushar, L.M. Foraging Behavior of the Timber Rattlesnake, Crotalus horridus. Copeia 1984, 1984, 976–981. [Google Scholar] [CrossRef]
  81. Palmer, W.M.; Braswell, A.L. Reptiles of North Carolina; The University of North Carolina Press: Chapel Hill, NC, USA, 1995. [Google Scholar]
  82. Trauth, S.E.; McAllister, C.T.; Worth, F. Vertebrate Prey of Selected Arkansas Snakes. J. Ark. Acad. Sci. 1995, 49, 188–192. [Google Scholar]
  83. Parmley, D.; Parmley, A.M. Food habits of the Canebrake Rattlesnake (Crotalus horridus atricaudatus) in Central Georgia. Ga. J. Sci. 2001, 59, 172. [Google Scholar]
  84. Sajdak, R.A.; Bartz, A.W. Crotalus horridus (timber rattlesnake). Arboreality, diet. Herpetol. Rev. 2004, 35, 60–61. [Google Scholar]
  85. Wittenberg, R.D. Foraging Ecology of the Timber Rattlesnake (Crotalus horridus) in a Fragmented Agricultural Landscape. Herpetol. Conserv. Biol. 2012, 7, 449–461. [Google Scholar]
  86. Hammerstein, G.; Steen, D.A.; Stevenson, D.J. Crotalus horridus (timber rattlesnake). Diet. Herpetol. Rev. 2015, 46, 640–641. [Google Scholar]
  87. Goetz, S.M.; Petersen, C.E.; Rose, R.K.; Kleopfer, J.D.; Savitzky, A.H. Diet and Foraging Behaviors of Timber Rattlesnakes, Crotalus horridus, in Eastern Virginia. J. Herpetol. 2016, 50, 520–526. [Google Scholar] [CrossRef] [Green Version]
  88. Schalk, C.M.; Trees, T.; Pierce, J.B.; Rudolph, D.C. Food habits of sympatric pitvipers from the West Gulf Coastal Plain, USA. Herpetol. Rev. 2018, 49, 1–5. [Google Scholar]
  89. Hewlett, J. Crotalus horridus (timber rattlesnake). Diet. Herpetol. Rev. 2019, 50, 589–590. [Google Scholar]
  90. Pynne, J.T.; Castleberry, S.B.; Conner, L.M.; Parsons, E.I.; Gitzen, R.A.; Austin, J.D.; Duncan, S.I.; McCleery, R.A. Timber Rattlesnake (Crotalus horridus) Predation on a Southeastern Pocket Gopher (Geomys pinetis). Southeast. Nat. 2019, 18, N34–N36. [Google Scholar] [CrossRef]
  91. Rautsaw, R.M.; Schramer, T.D.; Acuña, R.; Arick, L.N.; DiMeo, M.; Mercier, K.P.; Schrum, M.; Mason, A.J.; Margres, M.J.; Strickland, J.L.; et al. Genomic Adaptations to Salinity Resist Gene Flow in the Evolution of Floridian Watersnakes. Mol. Biol. Evol. 2021, 38, 745–760. [Google Scholar] [CrossRef]
  92. De Queiroz, K. Species concepts and species delimitation. Syst. Biol. 2007, 56, 879–886. [Google Scholar] [CrossRef] [Green Version]
  93. Nosil, P. Speciation with gene flow could be common. Mol. Ecol. 2008, 17, 2103–2106. [Google Scholar] [CrossRef] [PubMed]
Toxins 13 00782 g001 550
Figure 1. The type B Jones Center (JC) population showed significant venom expression and genetic divergence from the type A Osceola National Forest (ONF) population, suggesting incipient venom-associated speciation. (A) Reversed-phase high-performance liquid chromatography (RP-HPLC) profiles for seven adult individuals from each population are shown. The highlighted region from ∼100–120 min includes the snake venom metalloproteinases (SVMPs), and the two phospholipase A 2 subunits of canebrake toxin were found prior to 60 min [37]. All individuals from JC were monomorphic for type B venom, and all individuals from ONF were monomorphic for type A venom. (B) A time-calibrated BEAST analysis of cytochrome b sequences from both populations showed reciprocal monophyly and deep divergence between the populations. Median divergence time estimates (with 95% credible intervals) and clade support values (posterior probabilities) are shown for nodes where the posterior probability >0.70. Crotalus adamanteus (Cadam) was used as the outgroup. Type A individuals are represented by red dots, and type B individuals are represented by blue dots.
Figure 1. The type B Jones Center (JC) population showed significant venom expression and genetic divergence from the type A Osceola National Forest (ONF) population, suggesting incipient venom-associated speciation. (A) Reversed-phase high-performance liquid chromatography (RP-HPLC) profiles for seven adult individuals from each population are shown. The highlighted region from ∼100–120 min includes the snake venom metalloproteinases (SVMPs), and the two phospholipase A 2 subunits of canebrake toxin were found prior to 60 min [37]. All individuals from JC were monomorphic for type B venom, and all individuals from ONF were monomorphic for type A venom. (B) A time-calibrated BEAST analysis of cytochrome b sequences from both populations showed reciprocal monophyly and deep divergence between the populations. Median divergence time estimates (with 95% credible intervals) and clade support values (posterior probabilities) are shown for nodes where the posterior probability >0.70. Crotalus adamanteus (Cadam) was used as the outgroup. Type A individuals are represented by red dots, and type B individuals are represented by blue dots.
Toxins 13 00782 g001
Toxins 13 00782 g002 550
Figure 2. The type B Jones Center (JC) and the type A Osceola National Forest (ONF) populations differed significantly in head morphology. We measured interfang distance, head width, head length, and fang length and adjusted the values on the basis of snout-to-vent length. (A) We found significant differences in head width and fang length. The type A population had shorter fangs and narrower heads. Note that due to our transformation of the data, larger values on the y-axis represent smaller raw values. (B) A linear discriminant function analysis clearly distinguished the two populations on the basis of head morphology. Representative head and fang morphology for each venom type are shown. * p < 0.05 .
Figure 2. The type B Jones Center (JC) and the type A Osceola National Forest (ONF) populations differed significantly in head morphology. We measured interfang distance, head width, head length, and fang length and adjusted the values on the basis of snout-to-vent length. (A) We found significant differences in head width and fang length. The type A population had shorter fangs and narrower heads. Note that due to our transformation of the data, larger values on the y-axis represent smaller raw values. (B) A linear discriminant function analysis clearly distinguished the two populations on the basis of head morphology. Representative head and fang morphology for each venom type are shown. * p < 0.05 .
Toxins 13 00782 g002
Toxins 13 00782 g003 550
Figure 3. Type A venoms were found largely along the southern periphery of the range. We used a PCR assay to detect the presence of the two subunits of canebrake toxin to establish the presence of type A venom. The grey shading indicates the historical range of C. horridus. Our two focal populations (JC and ONF) are indicated. Some points were jittered in an attempt to show all individuals from a single locality. The inset map in the bottom right corner shows a clear contact zone between type A and type B individuals in the coastal plain region of South Carolina and Georgia.
Figure 3. Type A venoms were found largely along the southern periphery of the range. We used a PCR assay to detect the presence of the two subunits of canebrake toxin to establish the presence of type A venom. The grey shading indicates the historical range of C. horridus. Our two focal populations (JC and ONF) are indicated. Some points were jittered in an attempt to show all individuals from a single locality. The inset map in the bottom right corner shows a clear contact zone between type A and type B individuals in the coastal plain region of South Carolina and Georgia.
Toxins 13 00782 g003
Toxins 13 00782 g004 550
Figure 4. Representative reversed-phase high-performance liquid chromatography (RP-HPLC) profiles of venoms from additional populations showed that the A/B dichotomy was consistent across the range. Venoms either had peaks after 100 min and were therefore type B or peaks before 60 min that corresponded to the two canebrake toxin subunits.
Figure 4. Representative reversed-phase high-performance liquid chromatography (RP-HPLC) profiles of venoms from additional populations showed that the A/B dichotomy was consistent across the range. Venoms either had peaks after 100 min and were therefore type B or peaks before 60 min that corresponded to the two canebrake toxin subunits.
Toxins 13 00782 g004
Toxins 13 00782 g005 550
Figure 5. Environmental niche models for Type A and Type B individuals revealed significant differences in geographic and environmental distributions. Habitat suitability predictions produced by MaxEnt environmental niche models using 69 Type A (A) and 114 Type B (B) presences, indicated by points, are shown. The heatmap corresponds to habitat suitability scores and can be interpreted as a percentage.
Figure 5. Environmental niche models for Type A and Type B individuals revealed significant differences in geographic and environmental distributions. Habitat suitability predictions produced by MaxEnt environmental niche models using 69 Type A (A) and 114 Type B (B) presences, indicated by points, are shown. The heatmap corresponds to habitat suitability scores and can be interpreted as a percentage.
Toxins 13 00782 g005
Toxins 13 00782 g006 550
Figure 6. Venom types were not reciprocally monophyletic, suggesting gene flow between venom types in some geographic regions. We sequenced cytochrome b from 30 type A and 26 type B individuals and constructed a time-calibrated BEAST phylogeny. We identified two well-supported monophyletic clades of only type A animals, but we found no well-supported pure type B clades. The dot to the right of each tip corresponds to venom type (red type A, blue type B), and branch color corresponds to clade (distribution colored on map above). Top map inset shows clade color; bottom inset shows venom type. Median divergence time estimates (with 95% credible intervals) are shown at key nodes. Posterior probabilities >0.70 are shown. Crotalus adamanteus (Cadam) was used as the outgroup.
Figure 6. Venom types were not reciprocally monophyletic, suggesting gene flow between venom types in some geographic regions. We sequenced cytochrome b from 30 type A and 26 type B individuals and constructed a time-calibrated BEAST phylogeny. We identified two well-supported monophyletic clades of only type A animals, but we found no well-supported pure type B clades. The dot to the right of each tip corresponds to venom type (red type A, blue type B), and branch color corresponds to clade (distribution colored on map above). Top map inset shows clade color; bottom inset shows venom type. Median divergence time estimates (with 95% credible intervals) are shown at key nodes. Posterior probabilities >0.70 are shown. Crotalus adamanteus (Cadam) was used as the outgroup.
Toxins 13 00782 g006
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Margres, M.J.; Wray, K.P.; Sanader, D.; McDonald, P.J.; Trumbull, L.M.; Patton, A.H.; Rokyta, D.R. Varying Intensities of Introgression Obscure Incipient Venom-Associated Speciation in the Timber Rattlesnake (Crotalus horridus). Toxins 2021, 13, 782. https://doi.org/10.3390/toxins13110782

AMA Style

Margres MJ, Wray KP, Sanader D, McDonald PJ, Trumbull LM, Patton AH, Rokyta DR. Varying Intensities of Introgression Obscure Incipient Venom-Associated Speciation in the Timber Rattlesnake (Crotalus horridus). Toxins. 2021; 13(11):782. https://doi.org/10.3390/toxins13110782

Chicago/Turabian Style

Margres, Mark J., Kenneth P. Wray, Dragana Sanader, Preston J. McDonald, Lauren M. Trumbull, Austin H. Patton, and Darin R. Rokyta. 2021. "Varying Intensities of Introgression Obscure Incipient Venom-Associated Speciation in the Timber Rattlesnake (Crotalus horridus)" Toxins 13, no. 11: 782. https://doi.org/10.3390/toxins13110782

APA Style

Margres, M. J., Wray, K. P., Sanader, D., McDonald, P. J., Trumbull, L. M., Patton, A. H., & Rokyta, D. R. (2021). Varying Intensities of Introgression Obscure Incipient Venom-Associated Speciation in the Timber Rattlesnake (Crotalus horridus). Toxins, 13(11), 782. https://doi.org/10.3390/toxins13110782

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop