As Old as the hills: montane scorpions in Southwestern North America reveal ancient associations between biotic diversification and landscape history.

BACKGROUND: The age of lineages has become a fundamental datum in studies exploring the interaction between geological transformation and biotic diversification. However, phylogeographical studies are often biased towards lineages that are younger than the geological features of the landscapes they inhabit. A temporally deeper historical biogeography framework may be required to address episodes of biotic diversification associated with geologically older landscape changes. Signatures of such associations may be retained in the genomes of ecologically specialized (stenotopic) taxa with limited vagility. In the study presented here, genetic data from montane scorpions in the Vaejovis vorhiesi group, restricted to humid rocky habitats in mountains across southwestern North America, were used to explore the relationship between scorpion diversification and regional geological history.
RESULTS: Strong phylogeographical signal was evident within the vorhiesi group, with 27 geographically cohesive lineages inferred from a mitochondrial phylogeny. A time-calibrated multilocus species tree revealed a pattern of Miocene and Pliocene (the Neogene period) lineage diversification. An estimated 21 out of 26 cladogenetic events probably occurred prior to the onset of the Pleistocene, 2.6 million years ago. The best-fit density-dependent model suggested diversification rate in the vorhiesi group gradually decreased through time.
CONCLUSIONS: Scorpions of the vorhiesi group have had a long history in the highlands of southwestern North America. Diversification among these stenotopic scorpions appears to have occurred almost entirely within the Neogene period, and is temporally consistent with the dynamic geological history of the Basin and Range, and Colorado Plateau physiographical provinces. The persistence of separate lineages at small spatial scales suggests that a combination of ecological stenotopy and limited vagility may make these scorpions particularly valuable indicators of geomorphological evolution.

Introduction

The strength of the interplay between geological transformation and biotic diversification—geobiotic history [1]—is an old and contentious issue in biology [2]. With the advent of molecular-based estimates of lineage relationships and ages, traditional dichotomies between views that posit either strong adherence to ‘Earth and life evolving together’ or ‘life diversifying largely without dependence on Earth history events’ have been tempered by studies showing that species assemblages often exhibit a temporally wide range of idiosyncratic responses to Earth history events [3], [4], [5]. Although the age of lineages has become a fundamental datum for investigating relationships between Earth and biotic history, studies of diversification within the popular phylogeographical [6] paradigm have been biased towards relatively young lineages [7], [8]. Expansion of the temporal context of a study beyond the relatively recent timeframes encompassing phylogeography could capture a record of biotic diversification associated with geologically older landscape changes. Such changes may be retained in the genomes of ecologically specialized (stenotopic) taxa with limited vagility.

Owing to their antiquity, morphological, and ecological conservatism, scorpions are excellent candidates for investigating such deep geobiotic histories [9]. Derived from amphibious ancestors that lived more than 425 million years ago (Ma), scorpions represent an ancient radiation of terrestrial arthropods, earning the title ‘living fossils’ [10]. Their conserved morphology implies relative stasis in ecological requirements over time [9]. Many scorpions demonstrate limited vagility and a high degree of stenotopy such that the distributions of phylogenetically related taxa are often predictably restricted to a narrow range of stable habitats that may have persisted over millions of years [9], [11]. As evidence, several evolutionary studies on scorpions have revealed patterns of deep species-level divergence and diversification during the Miocene and Pliocene [12], [13], [14], [15], [16]. The diversity and ecological specialization of North American scorpions [17], [18] makes them appropriate for exploring biogeographical patterns on the continent.

In the present study, the historical diversification of montane scorpions in the Vaejovis vorhiesi group (Vaejovidae) was examined as part of a comparative approach aimed at reconstructing the paleohistory of highlands in southwestern North America [19], [20]. The vorhiesi group, comprising 12 described and several undescribed species, appears to be restricted to humid, rocky habitats in mixed pine-oak-juniper woodlands [21], [22] across their known distribution from central Utah southwards through Arizona and New Mexico to Sonora in northern Mexico (Fig. 1). Diversification in these scorpions appears to be constrained by habitat specialization, and divergence may have been promoted by allopatric speciation in isolated rocky highland habitats, as proposed for scorpion taxa in other parts of the world [9], [11]. If ephemeral woodland corridors that developed during Pleistocene glacial periods and connected highland biota [20], [23] lacked suitable humid, rocky habitats, dispersal between adjacent mountain ranges inhabited by stenotopic scorpions in the vorhiesi group, would be limited.

Figure 1

Collection localities for genetic samples of scorpions in the Vaejovis vorhiesi group included in this study.

White squares indicate type localities of described species (see Table 1), and grey line indicates approximate known distribution of the group.

Table 1

Collection data for genetic samples of scorpions in the Vaejovis vorhiesi group.

Sample ID	Taxon	Locality
Ajos SON	Vaejovis bandido	Mexico: Sonora: Sierra de los Ajos
Chiricahua AZ	Vaejovis cashi	USA: Arizona: Chiricahua Mountains*
Chiricahua 2 AZ	Vaejovis cashi	USA: Arizona: Chiricahua Mountains*
Chiricahua 3 AZ	Vaejovis cashi	USA: Arizona: Chiricahua Mountains*
Prescott AZ	Vaejovis crumpi	USA: Arizona: Prescott*
Prescott 2 AZ	Vaejovis crumpi	USA: Arizona: Prescott*
Santa Catalina AZ	Vaejovis deboerae	USA: Arizona: Santa Catalina Mountains*
Pinaleño AZ	Vaejovis electrum	USA: Arizona: Pinaleño Mountains*
Pinaleño 2 AZ	Vaejovis electrum	USA: Arizona: Pinaleño Mountains*
Pinaleño 3 AZ	Vaejovis electrum	USA: Arizona: Pinaleño Mountains*
Pinos Altos NM	Vaejovis feti	USA: New Mexico: Pinos Altos*
Mt Ord AZ	Vaejovis halli	USA: Arizona: Mazatzal Mountains, Mount Ord*
Mt Ord 2 AZ	Vaejovis halli	USA: Arizona: Mazaztal Mountains, Mount Ord*
Wupatki AZ IT	Vaejovis jonesi	USA: Arizona: Wupatki National Monument*
Flagstaff AZ IT	Vaejovis lapidicola	USA: Arizona: Flagstaff*
Payson AZ IT	Vaejovis paysonensis	USA: Arizona: Payson*
Hualapai AZ	Vaejovis tenuipalpus	USA: Arizona: Hualapai Mountains*
Hualapai 2 AZ	Vaejovis tenuipalpus	USA: Arizona: Hualapai Mountains*
Huachuca AZ	Vaejovis vorhiesi	USA: Arizona: Huachuca Mountains*
Huachuca 2 AZ	Vaejovis vorhiesi	USA: Arizona: Huachuca Mountains*
Animas NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Animas Mountains
Bonito AZ	Vaejovis sp., vorhiesi group	USA: Arizona: NE Flagstaff, near Bonito
Castle HS AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Castle Hot Springs
Chuska AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Chuska Mountains
Chuska 2 AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Chuska Mountains
Control Rd AZ	Vaejovis sp., vorhiesi group	USA: Arizona: S side of Control Road
Cookes NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Cooke's Range
Cookes 2 NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Cooke's Range
Daze Canyon AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Daze Canyon
Doña Ana NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Aguirre Springs National Recreation Area
Dragoon AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Dragoon Mountains
El Oso AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Mazatzal Mountains, El Oso Mine
Elenita SON	Vaejovis sp., vorhiesi group	Mexico: Sonora: Sierra Elenita
Florida Mtns NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Little Florida Mountains
Four Peaks AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Four Peaks
Gila Cliff NM	Vaejovis sp., vorhiesi group	USA: New Mexico: nr. Gila Cliff Dwellings
Hatchet NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Big Hatchet Mountains
Hillsboro NM	Vaejovis sp., vorhiesi group	USA: New Mexico: W Hillsboro
La Ventana NM	Vaejovis sp., vorhiesi group	USA: New Mexico: La Ventana Arch
Moab UT	Vaejovis sp., vorhiesi group	USA: Utah: Moab
Mule Creek NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Mule Creek
Mule Creek 2 NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Mule Creek
Mule Mtns AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Mule Mountains
Navajo Mtn UT	Vaejovis sp., vorhiesi group	USA: Utah: Navajo Mountain
Oak Creek AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Oak Creek, SW Flagstaff
Oak Creek 2 AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Oak Creek, SW Flagstaff
Peloncillo 2 NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Peloncillo Mountains
Peloncillo NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Peloncillo Mountains
San Francisco Mtns NM	Vaejovis sp., vorhiesi group	USA: New Mexico: San Francisco Mountains, Blue Range Wilderness
San Mateo NM	Vaejovis sp., vorhiesi group	USA: New Mexico: San Mateo Mountains
Santa Rita AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Santa Rita Mountains
Santa Rita 2 AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Santa Rita Mountains
Sawmill AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Sawmill
Signal Peak AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Signal Peak
Star Valley AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Star Valley nr. Payson
Strawberry AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Strawberry
Sunflower AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Sunflower
Tres Hermanas NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Tres Hermanas Mountains
Whetstone AZ CO, IT	Vaejovis sp., vorhiesi group	USA: Arizona: Whetstone Mountains
White Mtns AZ	Vaejovis sp., vorhiesi group	USA: Arizona: White Mountains near Pinetop
Williams AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Williams
Yarnell AZ	Vaejovis sp., vorhiesi group	USA: Arizona: Yarnell
Zuni NM	Vaejovis sp., vorhiesi group	USA: New Mexico: Zuni Mountains
V. franckei	Vaejovis franckei	Mexico: Oaxaca: Cerro Corral del Piedra
V. granulatus	Vaejovis granulatus	Mexico: Morelos: W of Huitzilac
V. montanus CHIH	Vaejovis montanus	Mexico: Chihuahua: Zorillo*
V. montanus CHIH 2	Vaejovis montanus	Mexico: Chihuahua: El Vergel
V. montanus SON	Vaejovis montanus	Mexico: Sonora: Yécora
V. nayarit	Vaejovis nayarit	Mexico: Nayarit: Jesus María Corte*
V. pattersoni	Vaejovis pattersoni	Mexico: Baja California Sur: Sierra de la Laguna*
Vaejovis sp. “Mezquital” DGO	Vaejovis sp., mexicanus group	Mexico: Durango: Mezquital
Vaejovis sp. “Mimbres” DGO	Vaejovis sp., mexicanus group	Mexico: Durango: Hwy 40, Rio Mimbres
Vaejovis sp. “del Nido”	Vaejovis sp., mexicanus group	Mexico: Chihuahua: Sierra del Nido
Paruroctonus boreus	Paruroctonus boreus	USA: Nevada: North Fork, Humboldt River

All samples are deposited in the Ambrose Monell Cryocollection (AMCC) at the American Museum of Natural History, New York. Asterisks denote samples collected in the vicinity of the type localities. Samples in bold were used in multilocus phylogenetic estimates (letter abbreviations indicate samples with missing data: CO, COI; IT, ITS2). Sequences were deposited in GenBank (accession numbers JX909361–JX909616).

Based on their relatively specialized ecological requirements for humid, rocky habitats and therefore close association with geological features, montane scorpions of the vorhiesi group may show signatures of pre-Pleistocene divergence, contrary to the phylogeographical patterns observed in many co-distributed taxa [23], [24], [25], [26]. These scorpions may thus offer novel insights into the deeper geobiotic history of the highland landscapes of southwestern North America. A mitochondrial DNA (mtDNA) dataset was generated from 63 samples of vorhiesi group scorpions and their relatives to explore phylogeographical structure, and a species tree reconstructed using multilocus data (mtDNA and two nuclear genes). Divergence times were estimated across the mtDNA dataset and the multilocus species tree, and the temporal distribution of divergence events was modeled across southwestern North America. These reconstructions and models were used to explore the relationship between diversification and regional geological history. Diversification within the vorhiesi group appears to have occurred almost entirely within the Neogene period, and may be related to the dynamic geological history of southwestern North America.

Methods

Ethics statements

Fieldwork in Mexico was conducted under permits granted by the Secretaría de Medio Ambiente y Recursos Naturales (SEMARNAT) to O.F. Francke, the late F. Mendoza-Quijano, and C. Solis-Rojas. The Navajo Nation and National Parks Service granted permits for collection in northern Arizona and New Mexico.

Genetic data

DNA sequence data were generated from 63 samples of vorhiesi group scorpions collected from throughout the known range (Fig. 1, Table 1). All described species in the group (V. bandido, V. cashi, V. crumpi, V. deboerae, V. electrum, V. feti, V. halli, V. jonesi, V. lapidicola, V. paysonensis, V. vorhiesi) were represented with the exception of V. bigelowi. The sister taxon to the group is uncertain, so additional samples from 10 other geographically proximate montane Vaejovis were included. Paruroctonus boreus was used as the outgroup [27].

Fragments of mitochondrial DNA were sequenced from a protein-coding gene (cytochrome c oxidase I, COI; 854 base pairs, or bp) and a ribosomal gene (16S rDNA, 16S; 404 bp). Two nuclear genes were also sequenced for a subset of samples (n = 39) representing geographically cohesive lineages within the vorhiesi group: 526 bp of 28S rDNA (28S) and 319 bp of the internal transcribed spacer region (ITS2) between the 5.8S rDNA and 28S rDNA. Taxon-specific mtDNA primers (Table S1) were used to avoid co-amplification of nuclear mitochondrial pseudogenes (numts) [28], and chromatograms examined for the presence of double peaks [29], indels, frameshifts, and premature stop codons. Primer sequences for nuclear genes were taken from Tully et al. [30] for 28S and Ji et al. [31] for ITS2. Genomic DNA was extracted from leg muscle tissue, and lab protocols provided in Bryson et al. [19] and Prendini et al. [32], [33] were followed to generate sequence data. Complete genetic data could not be obtained for several samples (Table 1). Sequence alignments for individual gene regions were performed with MAFFT v6 [34], [35] using default settings, the ‘1PAM/k = 2’ scoring matrix for nucleotide sequences, and the Q-INS-i algorithm for 16S and 28S data and the E-INS-i algorithm for ITS2 data. Heterozygous sites in nuclear segments were identified when two different nucleotides were present at the same position in electropherograms of both strands, with the weaker peak reaching at least 50% of the stronger signal. The gametic phases of the variants were computationally determined using PHASE v2.1.1 [36]. Five separate runs of 400 iterations were conducted for each nuclear dataset, and results with a probability threshold of 0.7 or higher accepted. All polymorphic sites with a probability less than 0.7 in both alleles were coded with the appropriate IUPAC ambiguity code. A partitioned homogeneity test with 1000 heuristic replicates was performed in PAUP* v4.0b10 [37] to test for conflicting phylogenetic signals between the mitochondrial and nuclear genes.

Phylogeographical estimation

The mtDNA dataset was analyzed to examine geographical structure and delineate geographically cohesive lineages within the vorhiesi group using Bayesian inference in MrBayes v3.1 [38]. Lineages were defined as genetically distinct geographical clusters with strong support values (≥0.95 Bayesian posterior probability), consistent with definitions of the term “phylogroup” [39], [40], [41]. Separate models were implemented for 16S and for each gene codon position of COI, the best partitioning scheme chosen based on Bayes Factors [42]. MrModeltest v2.1 [43] was used to select a best-fit model of nucleotide evolution, based on Akaike information criteria (AIC), for each partition. Bayesian settings included random starting trees, a variable rate prior, and heating temperature of 0.05. Analyses were run for 4 million generations, sampling every 100 generations, and output parameters visualized using the program TRACER v1.4 [44] to ascertain stationarity and whether the duplicated runs had converged on the same mean likelihood. All samples obtained during the first one million (25%) generations were discarded as burn-in, and a 50% majority-rule consensus phylogram with nodal posterior probability support estimated from the combination of the four runs post-burn-in.

Species tree and divergence date estimation

A species tree was reconstructed and divergence times within the vorhiesi group estimated from the multilocus dataset using *BEAST [45], a part of the BEAST v1.6.2 package. One to five exemplar samples were used from each geographical mtDNA lineage. Best-fit models of evolution were selected using MrModeltest, a Yule process speciation prior applied, and relaxed uncorrelated lognormal clocks specified for each gene tree. Trees were linked in analyses of the COI and 16S mtDNA data, which represent a single locus, and substitution and clock models were unlinked. The separation of the Cape region of Baja California from mainland Mexico between 7.5–14 Ma ([46], [47]; reviewed in [48]) was inferred as the causal mechanism for divergence between V. pattersoni and V. nayarit, putative sister species of montane Vaejovis isolated on opposite sides of the Sea of Cortez (Fig. 2). This split was used to calibrate the molecular clock. The xml file produced by BEAUTi was manually edited to set the divergence between V. pattersoni and V. nayarit to a normal distribution with a mean age of 10.75 Ma and standard deviation of 2, producing soft upper and lower bounds set to the maximum (14 Ma) and minimum (7.5 Ma) ages for the separation of the Cape from mainland Mexico. The ulcd.mean parameter for the four gene partitions was given wide uniform distributions (0.02 upper bound and 0.002 lower bound for the two mtDNA genes; 0.01 and 0.0001 for the two nuclear genes) to allow rates across the species tree to be estimated by constraints placed by the Baja calibration and not by rates placed on individual genes. Analyses were run for 200 million generations, with samples retained every 1000 generations. Results were displayed in TRACER to confirm acceptable mixing and likelihood stationarity, appropriate burn-in, and adequate effective sample sizes. The first 10% of generations were discarded as burn-in, and parameter estimates summarized on the maximum clade credibility tree using TreeAnnotator v1.6.2 [49]. This burn-in and visualization procedure was repeated for each of the three gene trees co-estimated by *BEAST.

Figure 2

Mitochondrial phylogeny of scorpions in the Vaejovis vorhiesi group.

Phylogeny inferred from Bayesian analyses of 1258 bp of concatenated COI and 16S mitochondrial sequence data. Posterior probability support values for nodes are indicated by coded dots explained in the figure legend. Inferred lineages are indicated by black bars. Samples used in species tree reconstruction are indicated in bold font. Localities of haplotypes are listed in Table 1. Inset depicts the geographical localities of putative sister species of montane Vaejovis isolated on opposite sides of the Sea of Cortez. The divergence of these two species caused by the split of the Cape region of Baja California from mainland Mexico was used to calibrate the molecular clock.

A relaxed Bayesian molecular clock framework implemented in BEAST was used in addition to estimate divergences from the complete mtDNA dataset. As with the mtDNA tree, the molecular clock was calibrated based on the split of the Cape region of Baja California from mainland Mexico. Vaejovis pattersoni and V. nayarit were grouped together as monophyletic, and the node representing the split of these two taxa was given a normal distribution with a mean age of 10.75 Ma and standard deviation of 2. Analyses were conducted for 40 million generations, with samples retained every 1000 generations, using a Yule tree prior. TRACER was used to confirm acceptable mixing and stationarity, appropriate burn-in, and adequate effective sample sizes. After discarding the first 4 million generations (10%) as burn-in, the parameter values of the samples from the posterior distribution were summarized on the maximum clade credibility tree using TreeAnnotator.

Diversification rate analysis

Temporal shifts in diversification rates were analyzed using Maximum Likelihood-based diversification-rate analysis [50], using divergence dates estimated from the multilocus dataset with *BEAST and from the mtDNA dataset with BEAST. The fit of different birth–death models implementing two constant rates (pure-birth and birth–death) and three variable rates (exponential and logistic density-dependent, and two-rate pure-birth) was computed with LASER v2.3 [51]. Model fit was measured using AIC scores. The significance of the change in AIC scores (ΔAICrc) between the best rate-constant and best rate-variable model was determined by creating a null distribution for ΔAICrc. This was accomplished by simulating 1000 trees using yuleSim in LASER with the same number of nodes and same speciation rate estimated under the pure-birth model. Log likelihood and AIC values were calculated for three models (SPVAR, EXVAR, and BOTHVAR; described in [52]) that permit differential extinction and speciation rates. In addition, lineage-through-time plots were generated to visualize the pattern of accumulation of log-lineages over time. Another set of analyses were performed on the mtDNA-only divergence date estimates from BEAST to assess the influence of additional genetic structure present within lineages that may have been excluded from the diversification rate analyses above. Eight additional divergences that occurred within six of the inferred geographical lineages (see Fig. S1) were included for this dataset.

Results

Genetic Data

Sequences were deposited in GenBank (JX909361–JX909616). The mtDNA genes contained 418 parsimony-informative sites (283/854 bp, COI; 135/404 bp, 16S) in the full dataset and 413 sites in the reduced species tree dataset (279, COI; 134, 16S). The nuclear genes contained fewer parsimony-informative sites (30/526 bp, 28S; 128/319, ITS2). No significant conflict was found between the mitochondrial and nuclear genes (P = 0.8). GTR+I+G models of sequence evolution were selected for the CO1 and 16S partitions in the mtDNA dataset. In the species tree datasets, GTR+I+G was selected for CO1 and 16S, GTR+I for 28S, and GTR+G for ITS2.

Phylogeographical analyses

Monophyly of the vorhiesi group was strongly supported in the mtDNA tree (Fig. 2). Strong phylogeographical signal was evident, but basal relationships were poorly resolved. Twenty-seven geographically cohesive lineages (including singletons, henceforth referred to as ‘lineages’ for convenience) were inferred from the mtDNA tree. Eleven of these lineages represented the species V. bandido, V. cashi, V. deboerae, V. electrum, V. feti, V. halli, V. jonesi, V. lapidicola, V. paysonensis, V. tenuipalpus, and V. vorhiesi and are referred to hereafter by their species names. The remaining 16 lineages consisted of samples from single mountain ranges with two exceptions. The exceptions were one lineage, referred to as ‘Mogollon’, that included samples from several mountains, and one lineage with samples from the adjacent Florida and Tres Hermanas mountains (Fig. 1).

Species tree and divergence date estimates

The species tree reconstruction (Fig. 3) revealed relatively deep structure with long terminal branches similar to the mtDNA gene tree used for phylogeographical inference (Fig. 2). The three gene trees each revealed relatively heterogeneous tree topologies (Fig. S1). Six geographically cohesive groups of lineages were consistently inferred, suggesting that the species tree was not heavily biased by the mtDNA data. Differences between the species tree and mtDNA tree were confined to weakly supported nodes (see Fig. S2 for species tree posterior probability support values), resulting in topological rearrangements towards the base of the trees. The phylogenetic positions of the paysonensis and tenuipalpus lineages, in particular, differed among the topologies.

Figure 3

Tempo of diversification for scorpions in the Vaejovis vorhiesi group.

(a) Time-calibrated species tree reconstructed using multilocus data. Numbers denote major geographical clades. These clades are mapped across the landscape in Figure 4. Bars indicate 95% highest posterior densities. Divergence times and posterior probability support values are provided in Figure S2. (b) Lineage through time plots derived from mtDNA and multilocus estimates of divergence dates. Birth–death likelihood analyses suggest a rate-variable density-dependent rate of diversification though time for both datasets. Approximate timing of major geological and climatic events that changed the landscape of southwestern North American are delineated, with time shown in millions of years (Ma). Pleist. = Pleistocene, QT = Quaternary.

Figure 4

Geography of diversification among scorpions in the Vaejovis vorhiesi group.

Color-coded dots, representing each of the six major geographical clades, are shown on the multilocus species tree (see Figure 3). Three lineages not within inferred major clades are labeled on the map. Sample names for each locality are provided in Figure 1. Gray bar indicates the approximate location of the Mogollon Rim.

Six geographically cohesive groups of lineages (referred to herein as ‘clades’; Clades 1–6, Figs. 3 and 4) and three lineages (paysonensis, tenuipalpus, and Doña Ana) formed the framework of the species tree. However, only the relationship between two clades (Clades 5 and 6) was strongly supported (Fig. S2). The geographical distribution of five of the six clades (Clades 1 and 4–6) and one of the two individual lineages (tenuipalpus) were anchored in the transition zone between the Colorado Plateau and the Basin and Range geological provinces, along the Mogollon Rim (Fig. 4). Two clades of Madrean ‘sky island’ endemics (Clades 2 and 3) were confined to the Basin and Range province. The tenuipalpus and Doña Ana lineages were respectively situated at the eastern and western extremes of the distribution of the vorhiesi group.

Divergence date estimates inferred from the Cape-calibrated multilocus and mtDNA datasets using *BEAST and BEAST were similar (Fig. S2). Mean divergence dates were generally within 2 Ma, whereas posterior credibility intervals were slightly narrower in the multilocus analyses. Divergence dates estimated across the species tree (Fig. 3 and Fig. S2) suggested that early diversification in the vorhiesi group began around 20 Ma. Although 95% posterior credible intervals were wide, most divergences estimated within the vorhiesi group were within the Neogene period from 25 to 2.6 Ma. Only five of the 26 inferred divergences incorporated the Pleistocene within credible intervals, and all mean divergence date estimates pre-dated the Pleistocene.

Tempo of diversification

Birth–death likelihood analyses of lineage diversification rates rejected the null hypothesis of rate-constancy for both datasets (P = 0.003, multilocus; P<0.001, mtDNA-only). The rate-variable model that best fit the datasets was the logistic density-dependent (DDL) model. Under this model, diversification rate in the vorhiesi group gradually decreased through time, with diversification rate estimated at 0.276 (multilocus *BEAST) or 0.370 (mtDNA-only BEAST) divergences per million years (Fig. 3). Models considering variable rates of extinction and speciation did not provide a better fit to the datasets. Similar results (rejection of the null hypothesis, P = 0.01; DDL best rate-variable model) were obtained when more recent divergences within lineages were accounted for (Fig. S2).

Discussion

Lineage diversification across southwestern North America

The time-calibrated phylogenetic hypothesis presented here suggests that scorpions in the vorhiesi group have had a long history in the highlands of southwestern North America. Results reveal a striking pattern of lineage diversification almost entirely confined within the bounds of the Neogene period (Fig. 3). Based on mean estimates, all 26 inferred cladogenetic events probably occurred prior to the Pleistocene epoch. Only five divergences incorporated the Pleistocene within credible intervals (Fig. S2).

Geological processes that cumulatively acted to separate the Basin and Range province from the Colorado Plateau occurred during a period known as the mid-Miocene transition [53], [54]. This episode of geological activity, 24–12 Ma, probably created the landscape in which the vorhiesi group originated. The transformation of an originally more continuous landscape to one of isolated mountain ranges associated with Basin and Range extension may have provided opportunities for rapid early diversification. The largely unsupported basal divergences in the species tree that occurred around 18–13 Ma (Fig. S2) are consistent with a rapid series of divergences. It is difficult to resolve these specific relationships with the current dataset, and whether this represents a case of simultaneous divergence from a common ancestor (i.e., hard polytomy) or estimation uncertainty (i.e., soft polytomy) remains unknown. Regardless, it is clear that the early divergences within these scorpions occurred in rapid succession over a narrow window of time (Fig. S2).

Following early diversification in the vorhiesi group, potentially linked with the mid-Miocene transition, the development of regional clades between about 12–7 Ma (Fig. 3) may have initiated in response to continued landscape deformation coincident with paleoclimatical change [55]. Basin and Range extension in southern Arizona may have spread rapidly northeastward across southern Arizona into the Rio Grande rift of southern New Mexico between about 13–9 Ma [56], [57], [58]. This development was synchronous with climatic change in the area that was punctuated by a dramatic shift in atmospheric conditions at 7–5 Ma [55], [59]. The distributions of four of the six inferred regional clades (Clades 2, 3, 5, and 6) and one lineage (Doña Ana) extend across southern Arizona and New Mexico (Fig. 4). The two remaining clades to the north (Clades 1 and 4) are centered along the Mogollon Rim. Late Miocene formation of localized volcanic fields across the Mogollon Rim [60] and tilting and faulting, around 13–6 Ma during the extension of the Basin and Range [54], may have contributed to diversification at the base of these two clades.

Fine-scale ancient endemism

Mountainous regions are known to generally harbor greater species endemism than adjacent lowland areas [61], [62], [63]. Much of this endemism resulted from Quaternary climatic oscillations driving fragmentation and isolation of highland habitats (e.g., [23]). Scorpions in the vorhiesi group display remarkable patterns of fine-scale endemism that appear to have developed much earlier than the Pleistocene. This is unsurprising given their stenotopic habitat requirements [21]. Diversification within the vorhiesi group was probably driven by allopatry in isolated rocky highland habitats: ephemeral woodland corridors may have connected highland biota during the Pleistocene but lacked suitable humid, rocky habitats required by these scorpions for dispersal between adjacent mountain ranges. The persistence of separate lineages at small spatial scales suggests that a combination of ecological stenotopy and limited vagility make the vorhiesi group particularly valuable as indicators of geomorphological evolution in southwestern North America, as observed in other scorpion taxa elsewhere in the world [9], [11].

Diversification rate in the vorhiesi group appears to have decreased slowly through time (Fig. 3). Given the restricted distributions and stenotopic requirements of vorhiesi group scorpions for insular montane habitats, this finding is consistent with studies suggesting that diversification rate slows as ecological opportunities diminish [52], [64], [65], [66], and that declining diversification rates may be caused by the saturation of geographical space [67]. The decreased diversification rate of these scorpions is probably related to the dynamic geological history of southwestern North America previously outlined, which appears sufficient to have formed new habitats and dispersal routes, and to have isolated older habitats throughout the Neogene. Divergence times in the vorhiesi group suggest that these highland habitats were occupied and persisted throughout the Neogene and that, by the onset of Quaternary, few unoccupied habitats would have been available. Quaternary climate change, while dramatically shifting vegetation communities across southwestern North America [68], may have had little impact on diversification in these stenotopic scorpions.

Montane scorpions as windows to pre-Quaternary history and orogeny

A range of genetic divergences across a Neogene temporal continuum, combined with fine-scale endemism, suggests that montane scorpions offer the potential to assess correlations between Earth history and biotic diversification. Millions of years of dynamic change across a heterogeneous landscape chiseled a genetic footprint in scorpions of the vorhiesi group that has not eroded through potential Pleistocene habitat connectivity. These scorpions can thus provide unique insight into the impact of pre-Quaternary events on the generation of biodiversity. Perhaps more importantly, these events can be inferred at a relatively fine geographical scale.

Agreement on a cohesive geohistorical framework for mountain formation across southwestern North America remains elusive despite decades of research [69]. Ecologically specialized ‘living fossils’, like the scorpions of the vorhiesi group, represent model organisms for tracking orogeny, as evidenced in the present study by the amount of genetic structure observed across the relatively small distances between adjacent mountain ranges. Other stenotopic scorpions may offer similar insights regarding mountain formation and the generation of diversity on montane islands.

Three gene trees for scorpions in the group embedded within the shared species tree. Strongly supported nodes (≥0.95 posterior probability support) are indicated with black dots.

(PDF)

Click here for additional data file.

Chronograms with estimated divergence times (in millions of years, Ma) for scorpions in the group. Estimates are shown for multilocus species tree (top) and mtDNA gene tree (bottom) analyses. Posterior probability support values for nodes are indicated by coded dots explained in the figure legend. Nodes that received <0.50 support are not indicated with dots. Means and 95% highest posterior densities (in brackets) are shown for each node. Asterisks indicate eight additional divergences included in diversification rate analyses (see Methods and Materials). QT = Quaternary.

(PDF)

Click here for additional data file.

Scorpion-specific mtDNA primers used in study on the group.

(DOC)

Click here for additional data file.