Speciation with gene flow poses a central paradox: how do genome-wide barriers to gene exchange accumulate as recombination continually breaks down associations among selected loci? Although theory predicts that together recombination, selection, and genome structure shape reproductive isolation, empirical studies often report conflicting patterns, suggesting that these determinants change across the speciation continuum. Here we compare genomic landscapes of introgression across rattlesnake lineages spanning a range of divergence. We generated a chromosome-level reference genome for the Southwestern Speckled Rattlesnake (Crotalus pyrrhus) and analyzed whole genome data from 181 individuals across two species complexes with a history of gene flow upon secondary contact. We show that reproductive isolation is highly polygenic and dynamically structured. At early divergence, introgression is most reduced in high recombination regions, consistent with increased efficacy of selection against gene flow at few large-effect loci. As divergence progresses, linked selection against gene flow dominates, generating a positive relationship between recombination and introgression expected to occur through the genome-wide coupling of polygenic barrier effects. Introgression landscapes also become increasingly correlated across species pairs as divergence increases due to repeated evolution of barriers in the same genomic regions. Here, we infer that the Z chromosome plays a prominent role in reproductive isolation, harboring a disproportionate number of barrier loci and showing reduced introgression even at early divergence. Together, these results reveal how recombination, selection, and genome organization interact to shape speciation with gene flow upon secondary contact, reconciling empirical patterns with predictions of speciation theory.

Island environments impose distinct ecological pressures that shape the evolution of complex traits. In snakes, body and head size are often related to diet, as gape limitation constrains the size of prey an individual can consume. Many venomous snakes exhibit ontogenetic venom expression variation that parallels shifts in prey selection over development, reflecting the integration of venom and morphology into a complex feeding phenotype. On Isla Cerralvo, the Baja pared to its mainland counterpart (Crotalus enyo enyo), potentially constraining prey size. We hypothesized that reduced head size would correspond with a lack of ontogenetic shift in venom expression in C. e. cerralvensis, as smaller heads may limit ontogenetic prey shifts. We sampled venom from island and mainland populations of C. enyo and assessed venom expression patterns across different size classes in each population. Although island and mainland rattlesnakes significantly differed a single venom phenotype independent of size. The absence of ontogenetic venom variation in island snakes, coupled with their smaller head sizes, suggests the integration of venom expression and cranial morphology as part of an integrated feeding phenotype and highlights the parallel evolution of functionally linked traits in response to unique ecological conditions.

Snake venoms are complex biochemical phenotypes whose composition can diverge dramatically among taxa, yet the functional and biomedical consequences of such variation remain understudied for many lineages. Here, we characterized the functional activities of venom from multiple rattlesnake species distributed across insular and mainland localities in the Baja California Peninsula and Gulf of California. We quantified phospholipase A2 (PLA2), proteolytic, and fibrinogenolytic activities as well as median lethal doses (LD50) for 55 individuals across Crotalus enyo, C. mitchellii, C. pyrrhus sensu lato, and C. ruber. We also evaluated the neutralization efficacy of the commercial polyvalent antivenom Antivipmyn®. Venom functional activities differed significantly among species. Crotalus mitchellii expressed a neurotoxic phenotype characterized by the highest PLA2 activity, the lowest LD50 values, and minimal proteolysis. Crotalus pyrrhus and C. ruber expressed proteolytic phenotypes with high enzymatic activity and ubiquitous dual-chain fibrinogenolysis. Crotalus enyo exhibited intermediate functional activities with moderate enzymatic output and the most efficient antivenom neutralization among species. Although Antivipmyn® neutralized the lethal activity of all venoms, efficacy differed significantly among species in both mass-normalized (mgAV/mgV) and potency-scaled (LD50-equivalents per vial) metrics. Notably, insular and mainland populations within species showed qualitative variation in both functional activities and antivenom responses. These results establish that the focal species possess functionally distinct venom phenotypes with species-specific consequences for antivenom neutralization, underscoring the relevance of venom functional characterization for clinical management of envenomation across the Baja California region.

Mountain organisms often exhibit complex distributions, patterns of variation, and evolutionary histories due to the topographical, climatic, and orogenic complexity of the areas in which they occur. Many montane taxa, however, lack the detailed sampling across space and the genome required to characterize diversity or estimate key evolutionary parameters. Here, we leverage comprehensive sampling, genome-wide data, and a draft genome assembly to investigate patterns of diversity and evolutionary history within a widespread but poorly studied snake, the Ridge-nosed Rattlesnake (Crotalus willardi), that inhabits the montane pine-oak woodlands in the Sierra Madre Occidental and Madrean Archipelago. We test for population genetic structure, measure gene flow and admixture between populations, and estimate phylogenetic history and divergence times among groups. We find substantial genetic structure, with up to seven differentiated populations, only three of which correspond with existing subspecific taxonomy. We identify one putative recent admixture event as well as evidence of historical gene flow among populations. We find that southern populations are generally older than those currently occurring in the Madrean Archipelago, which diverged within the last 0.6 million years. Together, our results reveal that current taxonomy fails to capture the substantial diversity and complex evolutionary processes at work in this species of montane snake and suggest that high-resolution sampling is crucial to understanding the biotic history of complex montane regions.

Understanding how human-mediated environmental change affects biodiversity is key for conserving evolvability. Because the most severe impacts are ongoing, such an understanding is proving exceptionally difficult to attain. Islands are natural, replicated experiments that serve as proxies for habitat fragmentation and, therefore, allow us to use historical changes in biodiversity under Island Biogeography Theory (IBT) to predict the consequences of immediate anthropogenic impacts on functional trait evolution. Rattlesnake venoms are molecular phenotypes that mediate interactions with prey, and diet and venom complexity are positively correlated. Consequently, rattlesnake venoms allow us to investigate how functional traits co-vary with changes in biodiversity according to IBT. We collected venom from 83 rattlesnakes across multiple species and 11 islands in the Gulf of California and estimated venom complexity using the Shannon Diversity Index. Using a mixed effects modeling approach, we found that the number of congenerics, island isolation, and island area best predicted venom complexity variability. All variables exhibited a negative relationship with venom complexity, contrary to predictions for island area under IBT. Larger islands with more congenerics exhibited reduced trait complexity, perhaps reflecting niche partitioning and venom specialization. Ultimately, we used a synthetic eco-evolutionary framework to predict functional trait evolution across fragmented landscapes.

The Crotalus molossus complex consists of five to seven phylogenetically related lineages of black-tailed rattlesnakes widely distributed in Mexico. While previous studies have noted venom variation within specific lineages of the Crotalus molossus complex, a comprehensive characterization of interspecific and ontogenetic venom variations, their functional implications, and the neutralizing ability of the Mexican antivenom against these variants remains largely unexamined. Herein, using two proteomic approaches for five lineages (C. basiliscus, C. m. molossus, C. m. nigrescens, C. m. oaxacus, and C. ornatus) of the C. molossus complex we characterized the number of toxins and their relative abundance in the venom of individuals of varying sizes. All five lineages undergo ontogenetic venom composition shifts associated with snake length. However, the pattern of ontogenetic shifts varied among lineages. In some lineages, these shifts led to significant differences in proteolytic, phospholipase A2, and fibrinogenolytic activities. Venom in smaller C. basiliscus, C. m. nigrescens, and C. m. oaxacus individuals had lower LD50 values (more lethal) in mice. Whereas the venom lethality of C. m. nigrescens (both juvenile and adult) and C. m. oaxacus (adult) was several times higher in a mammalian (mouse) model than in a reptilian (iguana) model. Antivipmyn® showed different neutralizing potencies toward venom pools. Overall, our results indicated that even among closely related rattlesnake lineages, venom phenotypes may vary greatly, impacting their function and the efficacy of antivenom neutralization.

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

Understanding the proximate and ultimate causes of phenotypic variation is fundamental in evolutionary research, as such variation provides the substrate for selection to act upon. Although trait variation can arise due to selection, the importance of neutral processes is sometimes understudied. We presented the first reference-quality genome of the Red Diamond Rattlesnake (Crotalus ruber) and used range-wide 'omic data to estimate the degree to which neutral and adaptive evolutionary processes shaped venom evolution. We characterized population structure and found substantial genetic differentiation across two populations, each with distinct demographic histories. We identified significant differentiation in venom expression across age classes with substantially reduced but discernible differentiation across populations. We then used conditional redundancy analysis to test whether venom expression variation was best predicted by neutral divergence patterns or geographically variable (a)biotic factors. Snake size was the most significant predictor of venom variation, with environment, prey availability, and neutral sequence variation also identified as significant factors, though to a lesser degree. By directly including neutrality in the model, our results confidently highlight the predominant, yet not singular, role of life history in shaping venom evolution.

Developmental phenotypic changes can evolve under selection imposed by age- and size-related ecological differences. Many of these changes occur through programmed alterations to gene expression patterns, but the molecular mechanisms and gene-regulatory networks underlying these adaptive changes remain poorly understood. Many venomous snakes, including the eastern diamondback rattlesnake (Crotalus adamanteus), undergo correlated changes in diet and venom expression as snakes grow larger with age, providing models for identifying mechanisms of timed expression changes that underlie adaptive life history traits. By combining a highly contiguous, chromosome-level genome assembly with measures of expression, chromatin accessibility, and histone modifications, we identified cis-regulatory elements and trans-regulatory factors controlling venom ontogeny in the venom glands of C. adamanteus. Ontogenetic expression changes were significantly correlated with epigenomic changes within genes, immediately adjacent to genes (e.g., promoters), and more distant from genes (e.g., enhancers). We identified 37 candidate transcription factors (TFs), with the vast majority being up-regulated in adults. The ontogenetic change is largely driven by an increase in the expression of TFs associated with growth signaling, transcriptional activation, and circadian rhythm/biological timing systems in adults with corresponding epigenomic changes near the differentially expressed venom genes. However, both expression activation and repression contributed to the composition of both adult and juvenile venoms, demonstrating the complexity and potential evolvability of gene regulation for this trait. Overall, given that age-based trait variation is common across the tree of life, we provide a framework for understanding gene-regulatory-network-driven life-history evolution more broadly.

Intraspecific phenotypic variation can be used as a window into the ecological differences among individuals of a species and lead to a better understanding of adaptive evolution. Adaptive traits, such as venom, that play an important ecological role for a species are useful models for understanding the sources of intraspecific variation. Intraspecific studies on front-fanged venomous snakes have offered deeper insights into the diverse mechanisms and adaptations that support the effectiveness of venom across species. Despite the extensive research on front-fanged venomous snakes, rear-fanged snakes, representing two-thirds of all snake species, have been largely overlooked. To test for sex and age-based intraspecific venom variation, we sequenced the messenger RNA from the Duvernoy’s gland of 9 male and 10 female Common Garter Snakes, Thamnophis sirtalis, of different sizes from a single location. Our data represent the most venom gland transcriptomes of any venomous snake species from a single location and represent the first Duvernoy’s venom gland transcriptomes for Thamnophis sirtalis. We found four toxin families dominate the Thamnophis sirtalis transcriptome: Snake Venom Metalloproteinases (SVMPs), Three-finger toxins (3FTxs), Cysteine-Rich Secretory Proteins (CRISPs), and C-type lectins (CTLs). Thamnophis sirtalis exhibits a unique balance in toxin expression, with approximately 30% each of neurotoxic (3FTx-dominated) and enzymatic (SVMP-dominated) components. No other published RFS Duvernoy’s gland transcriptome displays this ratio, rather they are dominated by one or the other. Additionally, venom expression varies with sex and size, with differences in toxin gene expression between males and females as they grow. Our study provides new insights on venom composition in a RFS species and highlights the amount of intraspecific variation possible among individuals from a single population.

To corroborate the ontogenetic shift in the venom composition of the Mexican Black-tailed Rattlesnake (Crotalus molossus nigrescens) previously reported through the census approach, we evaluated the shift in the protein profile, lethality, and proteolytic and phospholipase activities of four venom samples obtained in 2015, 2018, 2019, and 2021 from one C. m. nigrescens individual (CMN06) collected in Durango, Mexico. We demonstrated that the venom of C. m. nigrescens changed from a myotoxin-rich venom to a phospholipase A2 and snake venom metalloproteinase-rich venom. Additionally, the proteolytic and phospholipase activities increased with age, but the lethality decreased approximately three times.

We describe a new species of Eleutherodactylus (subgenus Syrrhophus) from Guerrero, Mexico, based on morphological and molecular data, as well as advertisement call analysis. Eleutherodactylus franzi sp. nov. has unique features including widely expanded fingertips, indistinct, but visible lumbo-inguinal glands, an immaculate white venter and dark reticulations on a cream dorsal background colouration. The new species belongs to the Eleutherodactylus nitidus species group. Eleutherodactylus franzi sp. nov. is micro-endemic, restricted to a small range in karstic hillsides on the southern extreme of the Mexican Transverse Ranges. We discuss conservation needs of this species, which we provisionally classify as Critically Endangered.

Biological specialization reduces the size of niche space while increasing efficiency in the use of available resources. Specialization often leads to phenotypic changes via natural selection aligning with niche space constraints. Commonly observed changes are in size, shape, behavior, and traits associated with feeding. One often selected trait for dietary specialization is venom, which, in snakes, often shows variation dependent on diet across and within species. The Neotropical Blunt-headed Treesnake (Imantodes cenchoa) is a highly specialized, rear-fanged, arboreal, lizard hunter that displays a long thin body, enlarged eyes, and a large Duvernoy's gland. However, toxin characterization of I. cenchoa has never been completed. Here, we use RNA-seq and mass spectrometry to assemble, annotate, and analyze the venom gland transcriptomes of four I. cenchoa from across their range. We find a lack of significant venom variation at the sequence and expression levels, suggesting venom conservation across the species. We propose this conservation provides evidence of a specialized venom repertoire, adapted to maximize efficiency of capturing and processing lizards. Importantly, this study provides the most complete venom gland transcriptomes of I. cenchoa and evidence of venom specialization in a rear-fanged snake, giving insight into selective pressures of venom across all snake species.

Despite the medical significance to humans and important ecological roles filled by vipers, few high-quality genomic resources exist for these snakes outside of a few genera of pitvipers. Here we sequence, assemble, and annotate the genome of Fea's Viper (Azemiops feae). This taxon is distributed in East Asia and belongs to a monotypic subfamily, sister to the pitvipers. The newly sequenced genome resulted in a 1.56 Gb assembly, a contig N50 of 1.59 Mb, with 97.6% of the genome assembly in contigs >50 Kb, and a BUSCO completeness of 92.4%. We found that A. feae venom is primarily composed of phospholipase A(2) (PLA(2)) proteins expressed by genes that likely arose from lineage-specific PLA(2) gene duplications. Additionally, we show that renin, an enzyme associated with blood pressure regulation in mammals and known from the venoms of two viper species including A. feae, is expressed in the venom gland at comparative levels to known toxins and is present in the venom proteome. The cooption of this gene as a toxin may be more widespread in viperids than currently known. To investigate the historical population demographics of A. feae, we performed coalescent-based analyses and determined that the effective population size has remained stable over the last 100 kyr. This suggests Quaternary glacial cycles likely had minimal influence on the demographic history of A. feae. This newly assembled genome will be an important resource for studying the genomic basis of phenotypic evolution and understanding the diversification of venom toxin gene families.

Here we report, for the first time, a natural hybrid between Crotalus atrox and C. mictlantecuhtli based on intermediate characteristics of the external morphology and venom. Morphologically, the individual had characteristics of both parent species. The hybrid's venom exhibited an intermediate composition including the presence of crotoxin which has never been documented in C. atrox but is well documented in C. mictlantecuhtli. The hybrid's venom was highly toxic and showed an intermediate proteolytic activity between the parental species. The two Mexican antivenoms were able to neutralize the hybrid's venom's lethality.

Traits for prey acquisition form the phenotypic interface of predator-prey interactions. In venomous predators, morphological variation in venom delivery apparatus like fangs and stingers may be optimized for dispatching prey. Here, we determine how a single dimension of venom injection systems evolves in response to variation in the size, climatic conditions and dietary ecology of viperid snakes. We measured fang length in more than 1900 museum specimens representing 199 viper species (55% of recognized species). We find both phylogenetic signal and within-clade variation in relative fang length across vipers suggesting both general taxonomic trends and potential adaptive divergence in fang length. We recover positive evolutionary allometry and little static allometry in fang length. Proportionally longer fangs have evolved in larger species, which may facilitate venom injection in more voluminous prey. Finally, we leverage climatic and diet data to assess the global correlates of fang length. We find that models of fang length evolution are improved through the inclusion of both temperature and diet, particularly the extent to which diets are mammal-heavy diets. These findings demonstrate how adaptive variation can emerge among components of complex prey capture systems.

Document containing supplemental figures S1-S4

Intraspecific variation in snake venoms has been widely documented worldwide. However, there are few studies on this subject in Mexico. Venom characterization studies provide important data used to predict clinical syndromes, to evaluate the efficacy of antivenoms and, in some cases, to improve immunogenic mixtures in the production of antivenoms. In the present work, we evaluated the intraspecific venom variation of Crotalus basiliscus, a rattlesnake of medical importance and whose venom is used in the immunization of horses to produce one of the Mexican antivenoms. Our results demonstrate that there is variation in biological and biochemical activities among adult venoms and that there is an ontogenetic change from juvenile to adult venoms. Juvenile venoms were more lethal and had higher percentages of crotamine and crotoxin, while adult venoms had higher percentages of snake venom metalloproteases (SVMPs). Additionally, we documented crotoxin-like PLA2 variation in which specimens from Zacatecas, Sinaloa and Michoacán (except 1) lacked the neurotoxin, while the rest of the venoms had it. Finally, we evaluated the efficacy of three lots of Birmex antivenom and all three were able to neutralize the lethality of four representative venoms but were not able to neutralize crotamine. We also observed significant differences in the LD50 values neutralized per vial among the different lots. Based on these results, we recommend including venoms containing crotamine in the production of antivenom for a better immunogenic mixture and to improve the homogeneity of lots.

Variation in gene regulation is ubiquitous, yet identifying the mechanisms producing such variation, especially for complex traits, is challenging. Snake venoms provide a model system for studying the phenotypic impacts of regulatory variation in complex traits because of their genetic tractability. Here, we sequence the genome of the Tiger Rattlesnake, which possesses the simplest and most toxic venom of any rattlesnake species, to determine whether the simple venom phenotype is the result of a simple genotype through gene loss or a complex genotype mediated through regulatory mechanisms. We generate the most contiguous snake-genome assembly to date and use this genome to show that gene loss, chromatin accessibility, and methylation levels all contribute to the production of the simplest, most toxic rattlesnake venom. We provide the most complete characterization of the venom gene-regulatory network to date and identify key mechanisms mediating phenotypic variation across a polygenic regulatory network.