Ethics statement

Scientific collecting permits in the United States were issued by the State of Arizona Game and Fish Department (SP628489, SP673390, SP673626, SP715023), the California Department of Fish and Wildlife (SC-12985), the New Mexico Department of Game and Fish (3563, 3576) and Texas Parks and Wildlife (SPR-0390-029). In Mexico, collection permits were issued by the Secretaria de Medio Ambiente y Recursos Naturales of the Estados Unidos Mexicanos (SEMARNAT: SGPA/DGVS/03562/15, SGPA/DGVS/01090/17, and FAUT-0015). Interactions with animals were approved by the University of Central Florida’s (UCF) Institutional Animal Care and Use Committee under protocol 13–17 W and followed the American Society of Ichthyologists and Herpetologists ethical guidelines.

Sample collection and DNA extraction

We collected representatives of all previously identified lineages38 of Mojave Rattlesnakes from throughout their distribution. For most samples collected in the field, we obtained venom and tissue. When possible, voucher specimens were created and deposited (Supplemental Table 1 with abbreviations following41). Tissues were stored in 95% ethanol or RNAlater and venom was collected and vacuum dried, frozen in liquid nitrogen, and/or stored at −80 °C. We collected a total of 216 individuals: 114 of these had tissue and venom, 34 had only venom, and 68 had only tissue (Supplemental Table 1). Whole genomic DNA was extracted from samples using the Serapure bead extraction protocol of Rohland and Reich42 following modifications in Faircloth43.

Reverse-phased High Performance Liquid Chromatography (RP-HPLC) to determine venom type

All venom was vacuum dried or lyophilized prior to use, resuspended in Millipore-filtered water, and centrifuged to remove insoluble debris. We determined protein concentration on a Qubit 3.0 Fluorometer (ThermoFisher Scientific) using the Qubit Protein Assay (ThermoFisher Scientific) following the manufacturer’s protocol. To determine the venom type (A, B, or A + B) of each individual, we used Reverse-phased High Performance Liquid Chromatography (RP-HPLC) based on the protocol in Margres et al.44. For RP-HPLC, we injected 100 μg of venom onto a Jupiter C18 column (250 × 2 mm; Phenomenex, Torrence, California, USA) using two solvents: 1 = 0.1% trifluoroacetic acid (TFA) in water and 2 = 0.075% TFA in acetonitrile. We used a Beckman System Gold HPLC (Beckman Coulter, Fullerton, California, USA) located in the Florida State University (FSU) Department of Biological Science Analytical Lab. The gradient started with 95% A and 5% B for 5 minutes followed by a 1% per minute linear gradient to 25% B, followed by a 0.25% per minute linear gradient to 55% B, a 2% per minute linear gradient to 75% B, a 14% per minute linear gradient to 5% B and then 5 minutes at the initial conditions all at a 0.2 mL/min flow rate. Run time was 180 minutes for each sample and the effluent was monitored at 220 and 280 nm45. Venoms were assayed using RP-HPLC for 148 individuals based on the presence of both subunits of MTX (Type A) and the presence of SVMPs (Type B) based on previous RP-HPLC profiles in C. scutulatus24,46 under these conditions5. When venom was available, RP-HPLC was used as the primary means of determining venom type.

Mojave toxin assay

To confirm venom type from RP-HPLC and to determine venom type when venom was unavailable, we used PCR assays for both subunits of Mojave Toxin (MTXA and MTXB). We amplified two fragments for each subunit using the primers designed by Zancolli et al.47. These primers were developed to determine MTX presence in C. scutulatus individuals from Arizona and New Mexico, USA and have also been shown to successfully amplify these fragments from individual C. scutulatus from Mexico22. PCR amplification was conducted on each DNA sample under the following conditions per 10 μL reactions: 3.5 μL PCR water, 1 μL 10X Sigma Buffer (Sigma-Aldrich, St. Louis, MO, USA), 1 μL of 25 mM MgCL 2 (Sigma-Aldrich), 1.3 μL of 2.5 mM each dNTPs, 0.5 μL of each primer at 10 μM, 0.2 μL of Taq Polymerase (Sigma-Aldrich), and 2 μL DNA. PCR was conducted on PTC200 Thermal Cycler (Bio-Rad, Hercules, CA, USA): 3.5 minutes at 94 °C, 35 cycles of 30 seconds at 94 °C, 1 minute at 57 °C, and 1 minute at 72 °C, and a final extension at 72 °C for 5 minutes. PCR product amplification was evaluated on a 2% agarorose gel using GelRed dye (Biotium, Fremont, CA, USA) to determine if the subunits were present. If one assay was positive, the individual was considered to have MTX and be Type A because all available data suggest that if these loci are present, they are also expressed in the venom22,47. We were unable to distinguish Type A and Type A + B individuals using the PCR assay because we could not determine the presence or absence of SVMPs without a venom sample.

Venom phylogeography of C. scutulatus

To determine if venom type correlated with phylogenetic lineage, we PCR amplified and sequenced the mitochondrial NADH4 (ND4) gene for any individual in our dataset not already sequenced in Schield et al.38. As outgroups, we included one sample each from C. viridis, C. cerberus, and C. oreganus, the sister species complex to C. scutulatus48 (Supplemental Table 1). We used the primers ND4 and Leu to sequence the partial ND4 gene as well as the adjacent His, Ser, and Leu tRNA genes49. PCR was conducted under the following conditions per 10 μL reaction: 3.8 μL PCR water, 1 μL 10X Sigma Buffer (Sigma-Aldrich), 1 μL of 25 mM MgCL 2 (Sigma-Aldrich), 0.8 μL of 2.5 mM each dNTPs (Invitrogen, Waltham, Massachusetts, USA), 0.5 μL of each primer at 10 μM, 0.4 μL of Taq Polymerase (Sigma-Aldrich) and 2 μL DNA. PCR was conducted on a PTC200 Thermal Cycler (Bio-Rad): 3.5 minutes at 94 °C, 35 cycles of 30 seconds at 94 °C, 1 minute at 53 °C, and 1 minute at 72 °C, and a final extension at 72 °C for 5 minutes.

Amplicons were purified by adding 0.5 μL FastAP (ThermoFisher Scientific #EF0651), 0.05 μL Exonuclease 1 (ThermoFisher Scientific #EN0581), and 7.45 μL PCR water to each 30 μL reaction and then placed on the thermal cycler for 30 minutes at 37 °C followed by 15 minutes at 85 °C. Sequencing was done in both directions at Eurofins Scientific (St. Charles, Missouri, USA) on an ABI 3730 genetic analyzer (Applied Biosystems, Waltham, MA, USA). Sequences were assembled and edited in Geneious v 10.1.2 (Biomatters Ltd., Auckland, New Zealand). Alignments were created with the MAFFT v 7.22 alignment algorithm50 implemented with default parameters in Geneious. We verified alignments by eye and trimmed low quality nucleotides and also checked to ensure there were no frameshift mutations. Our final alignment was 884 nucleotides for 190 ingroup and three outgroup taxa.

We used PartitionFinder 2.1.151 to determine the best-fit model of evolution for the ND4 alignment, partitioned by codon position and between protein-coding and tRNA regions. We used the “greedy” search algorithm, Bayesian Information Criterion (BIC), linked branch lengths, and only tested models available in MrBayes. The starting tree was generated using PhyML v 3.052. These analyses were run in the UCF Advanced Research Computing Center (ARCC) on the Stokes High Performance Computer (SHPC). Best fit models from PartitionFinder2 were HKY + I for the first codon position and the tRNA together, HKY for the second codon position, and HKY + Γ for the third codon position.

To determine the mitochondrial lineage of new samples in reference to Schield et al.38, we used Bayesian Inference (BI) in MrBayes v 3.2.653. Four independent MCMC runs were conducted for 107 generations and samples taken every 500 generations. The first 25% of each run was discarded as burnin and all runs were checked in Tracer v 1.654 to ensure stationarity was reached and that all ESS values for parameters from the individual and combined runs were ≥200. We combined the runs and generated a 50% majority rule tree.

Azocasein metalloproteinase assay

To determine if there are differences within and among venom types among populations, we performed an azocasein metalloproteinase assay on 146 samples in triplicate46. These assays were performed by incubating 20 μg of venom with 1 mg of azocasein substrate in buffer composed of 50 mM HEPES and 100 mM NaCl at a pH of 8.0 for 30 minutes at 37 °C. We stopped the reaction with 250 μL of 0.5 M trichloroacetic acid, vortexed, and brought it to room temperature. We then centrifuged it at 2000 rpm for 10 minutes. Sample absorbance was read at 342 nm and reported in Δ 342nm /min/mg of venom protein46. To determine the statistical significance of differences among samples at an α of 0.05, we used a Kruskal-Wallis test with venom type and lineage as factors implemented in R v. 3.4.355. If there was significance globally, we used a Nemenyi post hoc test implemented in the R package PMCMR to determine pairwise significance56.

Kallikrein-like serine protease assay

To test if there are differences within venom types for other toxin classes, we performed a kallikrein-like serine protease assay on 60 samples following the protocol of Mackessy57. This assay was conducted by adding 0.8 μg of whole venom to 373 μL of the same buffer as that used for azocasein metalloporoteinase assay described above. Samples were incubated for 3 minutes at 37 °C and then 50 μL of substrate (Bz-ProPheArg-pNA; Bachem, Torrance, CA, USA), the sample was vortexed and placed back at 37 °C for three minutes. The reaction was stopped with 50% acetic acid. Sample absorbance was read at 405 nm and the specific activity was calculated based on a standard curve of p-nitroaniline and reported as nanomoles of product produced per minute per mg of venom46. We used a Kruskal-Wallis test with venom type and lineage as factors in R to test for significant differences at an α of 0.05. If there was significance globally, we used a Nemenyi post hoc test implemented in the R package PMCMR to determine pairwise significance56.

SDS-PAGE protein gel electrophoresis

To characterize overall venom peptide diversity, we performed SDS-PAGE protein gel electrophoresis on 110 samples following the protocol of Smith and Mackessy46. We loaded 20 μg of whole venom into wells of a NuPAGE Novex bis-tris 12% acrylamide mini gel (Life Technologies, Grand Island, NY, USA) and elecrophoresed in MES buffer at 175 volts for 45 minutes. To estimate the molecular weight, we used 7 μL of Mark 12 standard. We stained gels overnight on a gentle shake with 0.1% Coomassie brilliant blue R-250 in 50% and 20% acetic acid (v/v). Gels were destained for approximately two hours in 30% methanol and 7% glacial acetic acid (v/v) in water until bands were clearly visible. Gels were gently shaken overnight at room temperature in 7% acetic acid (v/v) storage solution and imaged the following day using an HP Scanjet 4570c scanner.

Head morphology analysis

To determine if there are head morphological differences, we measured morphological characters for 57 individuals from Arizona, USA. We followed Margres et al.13 and measured SVL (snout to vent length), HL (head length), HW (head width), IF (interfang distance), and FL (fang length) for both fangs when not broken and then averaged them. We measured SVL with a tailor’s tape from the tip of the snout to the posterior end of the cloaca to the nearest 1 mm. We used IP54 digital calipers (iGaging, San Clemente, CA, USA) to measure HL, HW, IF, and FL to the nearest 0.01 mm. Head length was measured from the tip of the snout to the articular-quadrate joint, HW was the widest point behind the eye, IF was the distance between the two fang maxillae, and FL was from the top of the maxilla to the tip of the fang while folded. Average FL was determined if neither fang was broken. If one fang was broken, then the unbroken fang was used as the sole measurement; if both fangs were broken, that individual was not included in analyses of FL. All measurements were natural-log-transformed so they met the assumptions of normality and homoscedasticity. To size correct the data, we used the lnSVL value and subtracted each other value from it (ex. lnHL-lnSVL) to generate new columns for HL, HW, IF, and FL that were standardized based on the length of the snake. Using these data, we first compared all individuals based on presence or absence of Mojave Toxin which allowed us to include all samples. We then excluded animals in which we could not distinguish between Type A and Type A + B and compared individuals of the three venom types. All comparisons were done using Kruskal-Wallis tests at an α of 0.05. All transformations and analyses of morphological data were done in R v. 3.4.355. If there was significance globally, we used a Nemenyi post hoc test implemented in the R package PMCMR to determine pairwise significance56.

Niche modeling of venom type

To examine whether differences in venom type were correlated with differences in environmental variables, we constructed ecological niche models (ENMs) for the occurrence of Type A and Type B venom across the range of C. scutulatus. We used geographic localities for 81 type A and 68 type B C. scutulatus. Individuals that did not have venom but were positive for MTX were excluded from the analysis because we could not differentiate Type A vs Type A + B. ENMs were generated using MAXENT v 3.4.158 implemented in the R package dismo59. MAXENT uses occurrence records and a user-provided suite of environmental variables to predict the suitability of habitat and likelihood of occurrence across a landscape60. To limit the effect of sampling bias in the construction of EMS, we subsampled the total set of occurrence records to the same resolution as the our environmental data (30 arc seconds). This reduced the dataset to 72 and 54 representative Type A and Type B C. scutulatus occurrence records, respectively. For environmental data we used the 19 climatic variables collected in the WorldClim dataset v 1.461 as well as elevation, slope, and aspect with a 30 arc second resolution. To avoid biasing model fitting through inclusion of highly correlated inputs62, we removed 8 variables (BIO3, BIO5, BIO7, BIO10, BIO11, BIO13, BIO16, and BIO17) with a pair-wise Pearson’s correlation coefficient >0.90. BIO3 and BIO 7 were removed because they are derivatives of other bioclim variables and were correlated with them. BIO5 and BIO10 were removed because they were correlated with elevation. Elevation was kept because it was a variable specific to one of our hypotheses. BIO6 was removed in favor of BIO11 because temperature of the coldest month is less general than temperature of the coldest quarter. BIO13 and BIO16 were removed in favor of BIO12, the was more general. Lastly BIO14 was removed and BIO17 was kept, once again because BIO17 was more general. This left 14 remaining environmental variable characterizing southwestern North America. MAXENT models were run with 20 replicates and average model performance was evaluated by determining Area Under the Curve (AUC) for each model. ENMs were then generated for the Sonoran, Chihuahuan, and Central Mexican Plateau lineages independently following the same methods except only 10 replicates were used due to decreased sample size per lineage. This was done to further examine the correlations between environmental variables and venom types at finer geographic scales.

To test niche equivalency and niche similarity of Type A and Type B individual geographic distributions, we used the pseudoreplicate simulation method of Warren et al.63 as implemented in the R package phyloclim64. Both tests were run with 99 replicates to build a null distribution against which to test values of Schoener’s D and Warren’s I inferred from the full data models. To test niche equivalency, the occurrence points of two populations (e.g., individuals with Type A versus Type B venom) are combined and randomly partitioned into two datasets with sizes equal those of Type A and Type B. ENMs are generated for each dataset (e.g., a Type A model and a Type B model) and their similarity values (D and I) are calculated. This process is repeated to build a null distribution against which we test actual values of D and I inferred from each of the two full data models. The pseudoreplicate simulation method tests the null hypothesis that the two models (Type A and Type B) are not significantly different. In contrast, the test of niche similarity compares the ENMs of one venom type to an ENM of randomly selected subset of background cells (which include both presence and absence locations) of the other species. This is replicated 99 times and is then reversed such that the models for Type A venom are tested against background models for Type B venom and models for Type B venom are tested against background models for Type A venom. Comparison of the true D and I statistics to the pseudoreplicate distributions tests the hypothesis that one venom type’s niche model predicts the occurrences of the other venom type more than expected by chance. In addition, to more directly test for presence and absence differences as relating to specific environmental variables, we used logistic regression on the 14 variables used for ENMs. This approach was used to distinguish variables that may be important for C. scutulatus’ ecology, but may not be correlated with differences between venom types. Niche equivalency, niche similarity, and logistic regression of the model variables were then estimated for the Sonoran, Chihuahuan, and Central Mexican Plateau lineages independently to determine if the global pattern was the same at finer geographic scales.

Accession codes

MH883648-MH883754.