BackgroundSnake venoms have significant impacts on human populations through the morbidity and mortality associated with snakebites and as sources of drugs, drug leads, and physiological research tools. Genes expressed by venom-gland tissue, including those encoding toxic proteins, have therefore been sequenced but only with relatively sparse coverage resulting from the low-throughput sequencing approaches available. High-throughput approaches based on 454 pyrosequencing have recently been applied to the study of snake venoms to give the most complete characterizations to date of the genes expressed in active venom glands, but such approaches are costly and still provide a far-from-complete characterization of the genes expressed during venom production.ResultsWe describe the de novo assembly and analysis of the venom-gland transcriptome of an eastern diamondback rattlesnake (Crotalus adamanteus) based on 95,643,958 pairs of quality-filtered, 100-base-pair Illumina reads. We identified 123 unique, full-length toxin-coding sequences, which cluster into 78 groups with less than 1% nucleotide divergence, and 2,879 unique, full-length nontoxin coding sequences. The toxin sequences accounted for 35.4% of the total reads, and the nontoxin sequences for an additional 27.5%. The most highly expressed toxin was a small myotoxin related to crotamine, which accounted for 5.9% of the total reads. Snake-venom metalloproteinases accounted for the highest percentage of reads mapping to a toxin class (24.4%), followed by C-type lectins (22.2%) and serine proteinases (20.0%). The most diverse toxin classes were the C-type lectins (21 clusters), the snake-venom metalloproteinases (16 clusters), and the serine proteinases (14 clusters). The high-abundance nontoxin transcripts were predominantly those involved in protein folding and translation, consistent with the protein-secretory function of the tissue.ConclusionsWe have provided the most complete characterization of the genes expressed in an active snake venom gland to date, producing insights into snakebite pathology and guidance for snakebite treatment for the largest rattlesnake species and arguably the most dangerous snake native to the United States of America, C. adamanteus. We have more than doubled the number of sequenced toxins for this species and created extensive genomic resources for snakes based entirely on de novo assembly of Illumina sequence data.

Selection is predicted to drive diversification within species and lead to local adaptation, but understanding the mechanistic details underlying this process and thus the genetic basis of adaptive evolution requires the mapping of genotype to phenotype. Venom is complex and involves many genes, but the specialization of the venom gland toward toxin production allows specific transcripts to be correlated with specific toxic proteins, establishing a direct link from genotype to phenotype. To determine the extent of expression variation and identify the processes driving patterns of phenotypic diversity, we constructed genotype-phenotype maps and compared range-wide toxin-protein expression variation for two species of snake with nearly identical ranges: the eastern diamondback rattlesnake (Crotalus adamanteus) and the eastern coral snake (Micrurus fulvius). We detected significant expression variation in C. adamanteus, identified the specific loci associated with population differentiation, and found that loci expressed at all levels contributed to this divergence. Contrary to expectations, we found no expression variation in M. fulvius, suggesting that M. fulvius populations are not locally adapted. Our results not only linked expression variation at specific loci to divergence in a polygenic, complex trait but also have extensive conservation and biomedical implications. C. adamanteus is currently a candidate for federal listing under the Endangered Species Act, and the loss of any major population would result in the irrevocable loss of a unique venom phenotype. The lack of variation in M. fulvius has significant biomedical application because our data will assist in the development of effective antivenom for this species. N ATURAL selection can be a powerful force driving rapid diversification within species and is predicted to lead to local adaptation through the increase in frequency of mutations in gene-regulatory or protein-coding regions (Stern 2000;Hoekstra and Coyne 2007;Carroll 2008;Muller 2007). Expression variation at single loci has produced adaptive phenotypic divergence in the beaks of Darwin's finches (Abzhanov et al. 2004), coat color in mice (Manceau et al. 2011), and mimicry in butterflies (Reed et al. 2011). Most traits, however, are products of poorly characterized developmental pathways involving many loci. Many of the fundamental features of evolving systems, such as evolvability, epistasis, pleiotropy, and basic variational properties (Rokyta et al. 2008(Rokyta et al. , 2009(Rokyta et al. , 2011bWager 2008;Chou et al. 2011;Woods et al. 2011;Hill and Zhang 2012), result from the relationship between genotype and phenotype (Stadler et al. 2001;Hansen 2006), but the ability to study this relationship directly in polygenic traits is rare. Therefore, linking gene-regulatory changes to adaptive evolution in polygenic, complex phenotypes remains a challenge (Romero et al. 2012;Savolainen et al. 2013).Snake venoms are complex cocktails of 40-100 proteinaceous toxins (Boldrini-França et ...

Reconstructing species’ demographic histories is a central focus of molecular ecology and evolution. Recently, an expanding suite of methods leveraging either the sequentially Markovian coalescent (SMC) or the site-frequency spectrum has been developed to reconstruct population size histories from genomic sequence data. However, few studies have investigated the robustness of these methods to genome assemblies of varying quality. In this study, we first present an improved genome assembly for the Tasmanian devil using the Chicago library method. Compared with the original reference genome, our new assembly reduces the number of scaffolds (from 35,975 to 10,010) and increases the scaffold N90 (from 0.101 to 2.164 Mb). Second, we assess the performance of four contemporary genomic methods for inferring population size history (PSMC, MSMC, SMC++, Stairway Plot), using the two devil genome assemblies as well as simulated, artificially fragmented genomes that approximate the hypothesized demographic history of Tasmanian devils. We demonstrate that each method is robust to assembly quality, producing similar estimates of Ne when simulated genomes were fragmented into up to 5,000 scaffolds. Overall, methods reliant on the SMC are most reliable between ∼300 generations before present (gbp) and 100 kgbp, whereas methods exclusively reliant on the site-frequency spectrum are most reliable between the present and 30 gbp. Our results suggest that when used in concert, genomic methods for reconstructing species’ effective population size histories 1) can be applied to nonmodel organisms without highly contiguous reference genomes, and 2) are capable of detecting independently documented effects of historical geological events.