Snake venoms are complex cocktails of toxins acting synergistically to subdue and kill prey or predators. Venom composition is extremely variable not only between different species, but also within the same taxon causing grave problems to snakebite treatment. The underlying drivers and mechanisms of this variation remain poorly understood, particularly the relative importance of gene flow between populations. The Mohave rattlesnake, Crotalus scutulatus, displays extreme venom variation across a continuous distribution in Southwestern USA. Hence it represents an ideal model system to investigate both the genetic mechanisms underlying geographic venom variation and the extrinsic forces shaping it. Here we test whether venom variation in C. scutulatus (i) is the result of genetic differences at loci coding for toxin proteins, or (ii) is caused by differences in gene expression, and (iii) whether it reflects patterns of neutral genetic variation. We used proteomic analysis to characterize the venom phenotypes, a PCR assay to test for presence/absence of the most highly expressed toxin-encoding genes, and 16 microsatellite markers to infer population structure of C. scutulatus. Presence/absence of major venom proteins was strictly linked with presence/absence of the corresponding coding genes. Surprisingly, we found weak population structure and low genetic differentiation (Fst = 0.062), with no significant correlation between venom phenotypes and neutral genetic variation. These results suggest that forces other than neutral genetic drift are able to maintain marked differences of adaptive genes in the presence of gene flow. Furthermore, we revealed a genetic basis for snake venom variation across multiple toxin families.