The estimation of divergence dates between lineages has long been an area of interest in biogeographic investigations, as the timing of splits can provide significant insight into the evolutionary patterns and processes that have influenced taxonomic groups (Rambaut and Bromham 1998; Pyron and Burbrink 2013). For this research, the taxonomic group of interest is the genus Antrodiaetus (Araneae, Mygalomorphae, Antrodiaetidae), a group of folding-door (i.e., a type of trapdoor) spiders found throughout the Nearctic whose ecology and evolutionary history are deeply tied to their habitat—as is the case for almost all mygalomorph spiders. Mygalomorphs (i.e., tarantulas and trapdoor spiders) are an ancient group with fossil evidence dating to the Triassic Period (~250 million years before present [ybp]) (Dunlop et al. 2015), and their near-global distribution, linked with their restricted dispersal abilities, provides interesting insight into how biological diversity spread across the planet (Hendrixson and Bond 2007; Opatova et al. 2020; Biswas et al. 2023; Briggs et al. 2025).

Traditionally, the estimation of divergence dates from molecular data has used fossil date ranges as calibration points, with more fossils at different times across a phylogeny allowing for better estimation of evolutionary rates and divergence times (Donoghue and Benton 2007; Parham et al. 2012). Due to the burrowing lifestyle of mygalomorph spiders (i.e., living in a burrow reduces the likelihood of being trapped in layers of sediment), fossils are quite rare (Selden and Penney 2010). Being an old and diverse group, mygalomorph spiders have accumulated significant variation across their genomes, although this variation is not well understood (Král et al. 2013, 2019). Because of this, rates of molecular evolution likely vary across different taxonomic systems (Bromham et al. 2000; Bidegaray-Batista and Arnedo 2011; Azevedo et al. 2024; Briggs et al. 2026). Studies that have estimated rates of molecular evolution in mygalomorphs have used either a small number of mitochondrial or nuclear markers (Ayoub et al. 2007; Bond et al. 2012; Harrison et al. 2017; Rix et al. 2021) or, more recently, large genome-wide datasets composed of hundreds or thousands of loci (Garrison et al. 2016; Hedin et al. 2019; Opatova et al. 2020; Montes de Oca et al. 2022). Although genome-scale data have begun to be used in spider research, there is a dearth of knowledge of the rates of molecular evolution in mygalomorphs. This problem provides opportunities to investigate other means of calibrating divergence-dating methodologies.

One option is to use biogeographic events to calibrate the phylogeny, translating evolutionary rates from relative time into absolute time (Derkarabetian et al. 2010; Hipsley and Müller 2014; Landis 2017). Similar to using fossil calibrations, this approach requires certain assumptions about the taxonomic group in question. For example, one assumption is that abiotic forces strongly influence the evolution of the group, while another is that the biogeographic event can be reliably and independently estimated (Ho et al. 2015; Landis 2020). Therefore, knowledge of the taxonomic system and the history of a geographic region of interest is critically important when using such calibration points.

Antrodiaetus (Mygalomorphae, Antrodiaetidae) are fossorial spiders with low dispersal ability (i.e., no ballooning), long lives (i.e., estimated at 15–30 years), and long times to reach maturity; they exhibit strong habitat preferences at the species level and, such as most mygalomorph spiders (Mason et al. 2018), suffer from a depauperate fossil record. The only fossil that represents the family is Cretacattyma raveni (Eskov and Zohnstein 1990), a lineage related to the genus Atypoides (Dunlop et al. 2015)—the sister lineage to Antrodiaetus. This fossil places the group at least as far back as the Cretaceous Period (Eskov and Zohnstein 1990). Antrodiaetus possesses a disjunct distribution across the Holarctic, a pattern similar to that seen in other relatively old lineages from different taxonomic systems, and abiotic forces appear to have heavily influenced its evolution. The majority of species diversity occurs in the Pacific Northwest (PNW) of North America, which has a complex geological and ecological history (Graham 1993).

To gain an understanding of the biogeographic history of the PNW, the focus is on one of these “species,” Antrodiaetus pacificus, currently a cryptic species complex that exhibits an interesting disjunct distribution. The Pacificus species complex is a mesic, montane “species” with one lineage living in the Coastal and Cascade mountain ranges from northern California to Alaska and the others living throughout the interior forests of the Rocky Mountains in Idaho (Fig. 1). These lineages, and the coastal and interior forests they inhabit, are separated by the inhospitable Columbia Basin and Oregon Highlands. Long suspected of possessing cryptic species (Coyle 1971; Hendrixson and Bond 2007), the Pacificus species complex spans a well-studied biogeographic break where several other dispersal-limited species of flora and fauna can be found in forests on both sides, separated by inhospitable xeric habitat in between (Carstens et al. 2005; Brunsfeld et al. 2007). Competing hypotheses have been proposed to explain the disjunct distribution of species across this biogeographic break, broadly classified as either a pattern of older ancient vicariance (AV) or more recent inland dispersal (ID) (Brunsfeld et al. 2000). The ID hypothesis suggests dispersal from ancient coastal populations into the Rocky Mountains via the brief appearance of marginal mesic forest habitat at the edge of glacial retreat during the Pleistocene—with dispersal occurring via either a northern route through the modern-day Okanogan Highlands or a southern route through the Oregon Highlands (Soltis et al. 1997; Smith et al. 2018a; Smith et al. 2022). The AV hypothesis proposes that the orogeny of the Cascades and Coastal Ranges was followed by a subsequent xerification of the Columbia Basin of eastern Washington and Oregon, fragmenting a previously contiguous forest system across the PNW (Nielson et al. 2001; Carstens et al. 2004; Richart and Hedin 2013). This latter event is markedly older, would have taken place during the late Miocene or early Pliocene, and implies the existence of Pleistocene refugia for Rocky Mountain populations. It is important to note that both hypotheses make predictions about the topology of species relationships across the geographic break, with AV suggesting that coastal populations are more closely related, while the ID hypothesis suggests that inland populations are more closely related to either southern or northern coastal populations, depending on the dispersal route.

Previous work (Ciaccio et al. 2024) demonstrated that the break between A. pacificus lineages found east and west of the Columbia Basin was a much deeper split in the Antrodiaetus phylogeny than suspected, with the two lineages (referred to as “pacificus WEST” and “pacificus EAST”) most likely representing different species. However, this previous research identified a small population of “pacificus WEST” residing at the southernmost extent of the “pacificus EAST” distribution. The fragmentation of this lineage, referred to as “pacificus SE,” from “pacificus WEST” also matches the biogeographic break discussed above and appeared to correspond with the southern dispersal route of the ID hypothesis. Interestingly, the distribution of “pacificus SE” is not typical of other taxa that exhibit this geographic break, with “pacificus SE” possessing a very small distribution that is isolated farther south than is typically seen. This unique distribution might point to the possibility that neither of the events or timeframes outlined in the AV/ID hypotheses explains the divergence between “pacificus WEST” and “pacificus SE”—providing unique biogeographic insight for the PNW.

Because Antrodiaetus evolution is hypothesized to have been strongly impacted by abiotic forces, the objective was to use biogeographic information to overcome the dearth of fossil data for this group and estimate molecular evolutionary rates and associated divergence dates for species in the genus. To do this, ultraconserved elements (UCEs) were employed to provide genome-wide phylogenetic insight into Antrodiaetus history. Specifically, the study asked whether some UCE loci were better suited (i.e., more “clock-like”) than others for rate estimation; whether the AV or ID model of divergence was a better fit for the data; and what the molecular rate of evolution was in the UCE loci under the best model of divergence. To answer these questions, the study focused on the Pacificus species complex because of its unique spread across the Antrodiaetus phylogeny and occupation of both sides of the PNW biogeographic break.

Samples used in DNA extraction, library preparation, and high-throughput sequencing came from recently collected samples curated in the Arthropod Molecular Systematics Lab at the University of Idaho and added to the recent dataset created by Ciaccio et al. (2024). DNA extractions were performed with Omniprep for Tissue DNA extraction kits (G-Biosciences), followed by qualitative and quantitative DNA assessments using gel electrophoresis and a Qubit 4 fluorometer (Life Technologies). Extractions were prepared for sequencing using an SRSLY NanoPlus Uracil+ Base Kit, ForShear Enzymatic Shearing Module, 96 UDI Primer Plate, and Clarefy DNA Purification Beads from Claret Bioscience. UCE loci were captured using the Spider 2K v.1 probeset kit from Arbor Biosciences. Illumina sequencing was carried out by Rapid Genomics on an SP flow cell (2 × 250 bp) and an S4 flow cell (2 × 150 bp), as well as by Novogene on the NovaSeq X Plus series (2 × 150 bp). Raw reads were cleaned and processed using the phyluce pipeline (Faircloth 2016); contigs were then assembled using SPAdes (Bankevich et al. 2012), sequence alignment was carried out in phyluce using MAFFT (Katoh and Standley 2013), and internal trimming was performed by Gblocks (Castresana 2000) with settings of --b1 0.5, --b2 0.5, --b3 5, and --b4 10. Poorly aligned sequence fragments were cleaned with Spruceup (Borowiec 2019)—the final dataset used a 0.9 global cutoff.

Because divergence dating and estimation of molecular rates of evolution rely upon confidence in the inferred relationships, multiple approaches were used to build trees and assess confidence. Concatenated supermatrix phylogenies were inferred in IQTree2 (Minh et al. 2020) using 50%, 85%, and 95% occupancy matrices to compare with previous inferences of Antrodiaetus topology (Miller and Coyle 1996; Hendrixson and Bond 2007; Ciaccio et al. 2024) and to confirm the placement of new Antrodiaetus specimens within the tree. The data were partitioned by UCE loci using MODELFINDER (Kalyaanamoorthy et al. 2017) within IQTree2. To assess branch support, 1,000 ultrafast bootstrap and SH-aLRT replicates were performed. Individual gene trees were generated using IQTree2 for the same matrices. Species trees were inferred using ASTRAL-III (Zhang et al. 2018) and the 50%, 85%, and 95% datasets under the default heuristic search pattern. As Ciaccio et al. (2024) found a lack of concordance at certain nodes within the backbone of the Antrodiaetus phylogeny, a phylogenetic network was inferred using SplitsTreeApp (Huson and Bryant 2006) to visualize relationships in a non-bifurcating way. Networks were inferred using the 50% and 85% occupancy datasets, where a distance matrix with all tips was generated using P distance and a Neighbor Net approach, with the network then bootstrapped 100 times. The minimum percentage for registering a split used the default of 10%.

Selecting loci that are generally in agreement for parameters (e.g., topological arrangements) should increase the likelihood of a particular model, as variation among parameters should decrease (Doyle et al. 2015). Such a priori knowledge should improve the posterior probability in a Bayesian framework, as parameters will exhibit less variation in the data than without such data cleaning (Hedges et al. 1996; Kumar and Hedges 1998). Likewise, selecting loci that are more “clock-like” should decrease the uncertainty surrounding a posterior distribution of a molecular evolutionary rate (Jarvis et al. 2014). Recognizing that an assumption of “clock-likeness” is being introduced to the loci, the study chose to value a trade-off of capturing less of the overall UCE variation for greater specificity of signal.

The “gene shopping” approach described by Smith et al. (2018b) was employed, and the dataset was trimmed to include the most “clock-like” loci that reflected the preferred hypothesis of Antrodiaetus relationships. This data reduction and locus selection produced a dataset optimized for estimating a cohesive clock rate across taxa, as well as making analyses more computationally tractable. The analysis began with the 95% occupancy dataset, which reduced the overall number of loci while also decreasing the amount of missing data. Second, SortaDate (Smith et al. 2018b) was used to rank individual loci by user-determined parameters (i.e., root-to-tip variance, closeness of gene tree topology to a user-defined topology, and overall gene tree length). Loci were ranked by first matching the preferred topology, then minimizing root-to-tip variance, followed by selecting loci that possessed greater tree length (i.e., a discernible amount of evolution). The goal was to increase support for the divergence-dating inferences by removing loci with conflicting topologies and non-clock-like behavior. After the loci were ranked, the individual gene trees were visually examined, and the top loci were selected for the final dataset, as loci after the 20^th-best locus started to show higher root-to-tip variance.

Using the top 20 loci, a rate of molecular evolution was estimated, divergence dates were inferred, and disparate models were tested in the StarBeast3 package of BEAST 2 (Bouckaert et al. 2019; Douglas et al. 2022). Several runs were carried out initially to explore parameter settings, which were then organized into a set of models to be statistically compared through hypothesis testing. To evaluate the fit of these different models, the nested sampling (NS) package in BEAST 2 (Russel et al. 2019) was used to generate marginal likelihoods, as several of the prior parameters tested were mutually exclusive. Specifically, eight permutations of models were tested using the following parameters: two branching models (Yule and Fossilized Birth-Death [FBD]), a relaxed clock rate, a strict clock rate, an AV calibration for the time of divergence between “pacificus WEST” and “pacificus SE,” and an ID calibration for the same divergence (Fig. 4). Specifically, Model 1 used a strict clock, an ancient calibration, and an FBD model; Model 2 used a strict clock, an ancient calibration, and a Yule model; Model 3 used a relaxed clock, an ancient calibration, and an FBD model; Model 4 used a relaxed clock, an ancient calibration, and a Yule model; Model 5 used a strict clock, a recent calibration, and an FBD model; Model 6 used a strict clock, a recent calibration, and a Yule model; Model 7 used a relaxed clock, a recent calibration, and an FBD model; and Model 8 used a relaxed clock, a recent calibration, and a Yule model. For all models, the fossil calibration with Cretacattyma raveni (Eskov and Zohnstein 1990) was employed (parameters described below in the final run) using tip-date sampling. For all models, the biogeographic calibration (either ID or AV) created in BEAUTi was applied as a monophyletic node consisting of SE and the specimens immediately sister to it in the IQTree2 analysis (i.e., all WEST specimens excluding the Olympic Peninsula specimens and the lone southern Oregon specimen, as seen in Fig. 2). For the models using the addition of an ID hypothesis calibration to the fossil calibration, a lower bound of 0.001 (or 1,000 ybp) and an upper bound of 0.015 (or 15,000 ybp), with a mean of 0.005 (or 5,000 ybp), were used. For the models using the addition of an AV hypothesis calibration, a lower bound of 2 (or 2,000,000 ybp) and an upper bound of 20 (or 20,000,000 ybp), with a median of 5 (or 5,000,000 ybp), were used. The date ranges for both calibrations are based on previous investigation of the two hypotheses across several different taxa (Carstens et al. 2005) and are larger than those from Carstens et al. (2005) to allow for a thorough exploration of parameter space. Nested sampling was performed using the multi-threading option with four threads (one particle per thread), a subchain length optimally set for the data to 9,000 following (https://github.com/rbouckaert/NS-tutorial#nested-sampling-faq), and epsilon set to 1x10^–12.

Additional StarBeast3 runs were conducted to further test only the calibration parameter (i.e., AV and ID hypotheses) using Bayes factors (BF) following (https://www.beast2.org/2021/11/01/hypothesis-testing.html). This was performed by applying a biogeographic calibration prior that spanned the timeframe of both the AV and ID hypotheses, allowing direct comparison of the marginal likelihoods of a recent vs. ancient calibration point under either a relaxed or strict clock. For these runs, a uniform distribution bounded at 100 years (100 ybp) and 25,000,000 ybp was used to provide wide coverage of both hypotheses in parameter space. For one set of runs, a strict clock model was used, and for the other set of runs, a relaxed clock model was used—with both employing a Yule branching model. These analyses ran for 100,000,000 generations, sampling every 5,000 generations. When sampling from the prior-only and posterior runs for either the relaxed- or strict-clock model, samples in parameter space recorded for the most recent common ancestor (MRCA) between “pacificus WEST” and “pacificus SE” that were more recent than 1 million ybp were assigned as “hits” to the ID hypothesis. Any sampled point during the prior-only and posterior runs that was older than 5 million ybp was assigned to the AV hypothesis. Support for these hypotheses was assessed using Bayes factors.

Following model testing, a model for rate estimation and divergence dating was chosen after being parameterized accordingly (see Results and Discussion). To calibrate the phylogeny, the fossil Cretacattyma raveni was again used to bound the node that splits Antrodiaetus and Atypoides to the upper half of the Lower Cretaceous (~120–100 million ybp). The calibration used a log-normal distribution with a mean of 110 million ybp and a standard deviation of 0.5, with the mean in real space. The model was calibrated using the best-fitting biogeographic hypothesis (AV or ID) for “pacificus SE” and “pacificus WEST,” as determined from model testing and tree topology assessment. Two replicate runs of the final model were performed, with MCMC chains proceeding for 200,000,000 generations and sampling every 10,000 generations. All StarBeast3 runs were assessed in Tracer (Rambaut et al. 2018) to check for convergence and stationarity.

The topologies from the concatenated datasets (i.e., 50% and 85% occupancy) and the coalescent-based analysis (Fig. 2, Suppl. material 1) are generally in agreement and consistent with the findings of Ciaccio et al. (2024). For the concatenated dataset, support using ultrafast bootstraps closely follows that using SH-aLRT, with the exception of a few nodes within species near the tips of the tree. The additional samples added to the Ciaccio et al. (2024) dataset all fall within previously identified lineages and provide important increases in the understanding of Antrodiaetus species distributions across the PNW. The relationships of currently recognized species and novel lineages do not differ between the two methods of inference or dataset occupancy. Within species, there are differences between analyses. For example, in the species tree, the northern samples of “pacificus WEST” are more closely related to each other (with the exception of samples from the Olympic Peninsula) than to the southern “pacificus WEST” samples. The “pacificus SE” lineage is sister to both, whereas the Olympic Peninsula populations and the lone southern Oregon specimen are sister to all other “pacificus WEST” (Suppl. material 1). This contrasts with the concatenated topology, where northern samples of “pacificus WEST” do not form a clade (Fig. 2). Importantly, in both analyses, the “pacificus SE” lineage is always sister to the northern and southern “pacificus WEST” samples. The support for the relationships within these northern populations is relatively poor in both the concatenated and species-tree phylogenies (Fig. 2, Suppl. material 1). The position of the Olympic Peninsula specimens does not change, nor does that of the lone individual (ANTRO_0202) from southern Oregon between different phylogenetic inferences. The resulting phylogenetic network from the 85% matrix demonstrates the discordance (via reticulations) among UCE loci at certain nodes of the Antrodiaetus tree, inferring 413 splits. The network clearly identifies the main species lineages and species groups, although there are relatively short branches and numerous reticulations present. Specifically, many of these reticulations occur within what is now referred to as the “Inland PNW” group, which includes A. cerberus, A. pugnax, A. robustus, A. hageni, and “pacificus NE”—leading to lower certainty surrounding those nodes in both the concatenated and coalescent-based species trees (Fig. 2, Suppl. material 1).

Of 1,671 UCE loci recovered, a 95% occupancy dataset yielding 71 UCEs was produced. After visual examination of gene trees revealed a pattern of deviation from the preferred Antrodiaetus topology (i.e., the best estimation of relationships), the dataset was further reduced by ranking the loci with SortaDate based on increasing variation in tree length and root-to-tip variance. The goal was not to explicitly test the efficacy of “clock-like” UCEs vs. non-“clock-like” UCEs for divergence dating and rate estimation—only to find a subsample of loci with similar evolutionary properties that could feasibly be used for hypothesis testing of divergence dates and molecular evolutionary rate estimation due to computational restrictions. The top 20 loci were selected as the optimized dataset.

BF testing strongly supports the AV hypothesis (AV model with a strict clock, BF = 25.7; ID model with a strict clock, BF = 0.039)—see Kass and Raftery (1995) for an explanation of BF support. BF testing also strongly supports the AV hypothesis with a relaxed clock model (BF = 10.5), compared with weak support for the ID hypothesis with a relaxed clock (BF = 0.1). NS hypothesis testing of the eight different models ranked Model 3 as the best (with a marginal likelihood of –114,507), followed by Model 7 (–114,591) and Model 1 (–114,645) (Fig. 4). Model 3 (the NS best-ranked model) employs a relaxed molecular clock, an FBD branching model, and an AV hypothesis calibration. Although Models 3 and 7 were ranked higher, Model 1 was employed for the estimation of rates and dates for several reasons. First, stationarity could not be achieved with high enough confidence in the effective sample sizes (ESS) during StarBeast3 runs. Initial runs with those models showed no signs of stationarity at 100 million generations, and extending the runs to 200 million generations did not yield any marked improvement. Second, multiple runs of both models showed no convergence in the posterior distribution. Moreover, the runs using Models 3 and 7 output topologies that were markedly divergent from the known topological understanding. These topologies contained multiple polytomies, with species relationships that varied across the trees between runs, as well as trees containing negative branch lengths. These behaviors were seen only in those model permutations using a relaxed clock model and were particularly noticeable in Models 3 and 7. Therefore, Model 1, the next-best model, which did reach stationarity and achieve convergence between runs and did not exhibit the behaviors mentioned above, was used. Model 1 employs a strict clock, an FBD branching model, and an AV hypothesis calibration.

Under Model 1, a mean strict clock rate of 4.499 × 10^–4 substitutions per site per million years was estimated as the molecular rate of evolution for Antrodiaetus and Atypoides (with a 95% HPD of 2.747 × 10^–4–6.328 × 10^–4). Importantly, Model 3—the best-ranked model from NS hypothesis testing, although it did not achieve stationarity—yielded a similar rate (mean = 3.9693 × 10^–4, 95% HPD = 1.2639 × 10^–4–1.1811 × 10^–3; see Suppl. material 2). Under the AV hypothesis, the mean split between “pacificus WEST” and “pacificus SE” was estimated to be 8.5086 million ybp (95% HPD = 6.3372–10.6855).

The split between Antrodiaetus and Atypoides is deep, occurring during the Late Jurassic or Early Cretaceous, ~138 million ybp (95% HPD = 113.81–164.98). Divergence dates within Antrodiaetus are primarily Cenozoic (Fig. 3), occurring within the Eocene and later. The phylogeny placed A. roretzi as paraphyletic with respect to the rest of Antrodiaetus and sister to Atypoides. However, due to the unique morphology, age, and distribution of A. roretzi, no conclusions are drawn from this placement other than the need to more thoroughly sample this species. The Lincolnianus group diverged approximately 43.49 million ybp (95% HPD = 34.63–52.40), while A. montanus and “pacificus WEST” + “pacificus SE” diverged ~38.74 million ybp (95% HPD = 30.86–47.15). The Inland PNW group diverged from the A. unicolor + A. microunicolor + A. yesoensis group ~35.36 million ybp (95% HPD = 27.98–43.31). Ages of species splits between the main species lineages vary but, in general, are no more recent than the Miocene. Overall, the ranges of Antrodiaetus divergence date estimates (i.e., error bars) are not wide. For example, the 95% HPDs for divergence estimates range from ~4 million ybp (split between “pacificus WEST 1” and “pacificus SE,” 6.3–10.7 million ybp) to ~19 million ybp (“pacificus EAST” from the rest of Antrodiaetus, 37.3–56.3 million ybp).

Antrodiaetus evolution is hypothesized to have been strongly impacted by abiotic forces due to their long times to reach maturity, low dispersal abilities, strong habitat fidelity among species, and disjunct distributions. Due to a scarcity of fossils for Antrodiaetus or closely related lineages, the objective was to estimate a molecular evolutionary rate and divergence dates using biogeographic information in a hypothesis-testing framework. This research expanded upon Ciaccio et al. (2024) by adding twice as many samples, thus providing more coverage across species and their distributions. These new data confirm the previous work’s phylogenetic inference, as the relationships of the major species groups do not change (i.e., the Inland PNW, the Lincolnianus group, A. montanus + “pacificus WEST” + “pacificus SE,” A. unicolor + A. microunicolor + A. yesoensis, and “pacificus EAST”). The Japanese A. yesoensis is again confirmed as a lineage embedded within the North American Antrodiaetus phylogeny and sister to eastern Nearctic taxa. The “pacificus EAST” lineage is still an independent lineage sister to all other Antrodiaetus in North America, and A. robustus is still embedded within the Inland PNW group. These relationships are recovered in both the concatenated and coalescent-based analyses, with minor topological differences near the tips of the tree corresponding to populations within species (Fig. 2, Suppl. material 1).

The lack of concordance among UCE loci for some relationships, seen previously in Ciaccio et al. (2024), is also identifiable in the phylogenetic network (inset of Fig. 2). Non-bifurcating splits and reticulations at the base of the Inland PNW group, as well as with A. unicolor + A. microunicolor + A. yesoensis and the Lincolnianus group with “pacificus EAST,” highlight this—likely due to the short internal branches at these nodes. These short internal branches were identified previously in Ciaccio et al. (2024), and the expanded dataset has not meaningfully increased the resolution of those nodes across datasets or analyses. These types of short internal branches and lack of support (i.e., incomplete lineage sorting) are frequently inferred to result from rapid lineage splits that occurred so “quickly” (in an evolutionary sense) that even genomic-scale data do not harbor enough nucleotide differences for phylogenetics to effectively sort lineage splits (Degnan and Rosenberg 2006; Whitfield and Lockhart 2007; Philippe et al. 2011). Despite the discordance, the main species groups are cohesive (e.g., the Inland PNW group), as is the placement of “pacificus SE” with regard to the northern and southern populations of “pacificus WEST”—the pattern seen in the AV biogeographic hypothesis.

The inferred topologies support an ancient (AV) pattern of vicariance for “pacificus WEST” and “pacificus SE.” The placement of “pacificus SE” as sister to the “pacificus WEST” northern and southern populations does not change across phylogenetic analyses. Disagreement is primarily seen in the cohesiveness of the northern samples from the Coastal and Cascade Ranges in “pacificus WEST,” and the exact relationships of those populations remain uncertain. Furthermore, the notable Olympic Peninsula population, as well as a lone sample from southern Oregon, always appears as distinct groups that are supported in both the coalescent-based and concatenated analyses (Fig. 2, Suppl. material 1). Why these two populations are unique is difficult to determine at this time, and interpretation of the signal from them should be treated with caution until more sampling is added. An initial interpretation of the tree might suggest that the lone southern Oregon and Olympic Peninsula samples reflect contradictory support for the ID hypothesis; however, this is not considered to be the case. Previous research on this biogeographic break has shown that the topological patterns hypothesized by ID and AV can be less well-defined in empirical data than theory suggests, with this complexity possibly attributable to various causes (e.g., incomplete lineage sorting or relictual lineages) (Smith et al. 2018a). It is hypothesized that a similar situation has occurred in the Pacificus species complex. Given the age of the AV hypothesis, “pacificus WEST” is believed to have a complex phylogeographic history that generally follows the AV hypothesis but is punctuated by more nuanced breaks that could be attributable to smaller geographic events, a lack of time for reciprocal monophyly, or extirpation of local populations. As such, the data suggested that an older calibration point for this split was more appropriate for estimating an evolutionary rate from the UCE loci.

BF testing for the ID and AV calibration models strongly supports the fit of the AV model to the data. While the AV hypothesis under a relaxed clock model is not supported as strongly as under a strict clock, support for AV over ID is still clear. Testing values ranging from 100 ybp to 1 million ybp for a recent date and 5 million ybp to 25 million ybp for an ancient date would appear to skew the potential for the data to support the older prior; however, this was not the case. The posterior shows that much less sampling time is spent in parameter space at less than 1 million ybp. Most sampled points from the posterior (for both relaxed and strict clocks) are much higher, with the mean age of the posterior for the strict clock model being ~8.6 million ybp (95% HPD = 6.4–11.1 million ybp) and ~16.6 million ybp (95% HPD = 9.7–25 million ybp) for the relaxed clock model.

NS hypothesis testing and StarBeast3 runs also support an older age for calibration. NS ranking clearly supports a relaxed clock over a strict clock, an FBD branching process over a Yule branching process, and an ancient divergence over a recent split (Fig. 4). Unfortunately, StarBeast3 runs using a relaxed clock were not successful, as they could not achieve stationarity or convergence and generated trees with negative branch lengths, widely disparate topologies, and improbable topologies that made little sense given current knowledge of Antrodiaetus. Topological issues could be explained by discordance in the Antrodiaetus phylogeny at certain internal nodes, as the branches that often changed positions in these runs (relative to the preferred topology) were those that the phylogenetic network revealed to have many reticulations. The relaxed clock run with a recent calibration point (second-best model via NS; see Fig. 4) fared no better and inferred negative branch lengths along with a discordant topology. Relaxed clock trees, when unconstrained for topology, often estimated divergence dates similar to those of the strict clock model runs but with fluctuating, improbable topologies. Constraining the topology to the preferred topology caused the estimated dates to be much older, shifting them from Cenozoic events to Mesozoic events. Importantly, the top models of NS testing support an AV hypothesis over ID and an FBD branching model over a Yule model. Given the taxonomic diversity of the data, it is a fair assumption that the rate of evolution among these lineages is more similar than not. When the behaviors of Model 3 and Model 7 (described above) were evaluated in combination with the run times (even on a high-performance computing cluster), the authors determined that extending these runs (i.e., the number of generations) of these models would not improve the results and would diminish the tractability and reproducibility of the analysis. The authors suspect that the difficulty of using a relaxed clock model without strict assumptions about other parameters (e.g., more divergence-calibrated nodes) could be playing a role. An oversimplified interpretation would be that these models have “too many variables” to generate confidence. Therefore, Model 1 from NS testing, which applies a strict clock model under the AV calibration hypothesis and an FBD branching process, was chosen. The prior parameter for the “pacificus WEST”/“pacificus SE” split was broadly bounded to account for uncertainty.

The estimated rate is an order of magnitude slower than the nuclear rate that has previously been used for divergence dating in spiders (Ayoub et al. 2007; Opatova et al. 2013; Harvey et al. 2018), as well as the mitochondrial rate derived from Heliconius butterflies by Brower (1994) that has also been used for spiders (Hedin 2001; Hamilton et al. 2011; Rix and Harvey 2012). Other spider research that used different mitochondrial and/or nuclear markers has found faster rates than those inferred here (Ceccarelli et al. 2016). The estimated rate is also an order of magnitude slower than the average rate for UCEs found within Harmochirine jumping spiders (Araneae, Araneomorphae) (Azevedo et al. 2024), but similar to some of the individual UCE rates from those estimates. Azevedo et al. (2024) found higher variance in the individual rates of UCEs than that inferred here, whereas the evenness of the rates (see Suppl. material 2) is likely attributable to the “gene shopping” approach applied to the data and the data being focused on two sister genera (Antrodiaetus and Atypoides).

The slower rates inferred here are interesting, especially in comparison with the rates reported by Azevedo et al. (2024). Although UCEs have flanking regions of increasingly variable sites (i.e., more informative at the interspecific level), the “elements” are, by definition, quite conserved in nature. Furthermore, the data were restricted to the more conserved loci (i.e., the 95% matrix was the starting occupancy for gene shopping), and this might therefore be reflected in the data as a slower evolutionary rate. Research with UCEs in vertebrates and mammals has shown that UCE regions evolve an order of magnitude more slowly than coding regions (Cummins et al. 2024); however, as Hedin et al. (2019) pointed out, UCEs in arachnids (and arthropods in general) are primarily found within exonic regions (i.e., the core UCE region is more conserved and informative at deeper evolutionary relationships—unlike the flanking regions). A more likely explanation for the slower rates is that araneomorph spiders (such as the jumping spiders mentioned above) generally have very different life histories than mygalomorph spiders. For example, araneomorphs have much faster generation times than mygalomorphs, and they can also disperse via “ballooning,” which produces much more gene flow. These attributes certainly play a significant role in their faster molecular evolutionary rates.

The analyses suggest earlier divergences between many of the species of the Inland PNW group than previously hypothesized. Relatively rapid divergences between A. cerberus, A. pugnax, A. hageni, “pacificus NE,” and A. robustus occurred during the early Miocene or Oligocene/Miocene transition, a hypothesis more consistent with Hendrixson and Bond (2007) than with the later Pliocene hypotheses of Miller and Coyle (1996). The notable split of A. robustus from A. cerberus + A. pugnax across the Nearctic corresponds with the retreat of transcontinental forests that were once distributed across North America and were gradually replaced by grasslands throughout the Miocene. The splits between the other species in the Inland PNW occurred rapidly during the late Oligocene/early Miocene and are difficult to associate with particular geological events. During that period, a general global cooling transition from the hothouse Eocene to the cooler Pliocene was occurring, a time that saw the rise of grasses and savannas (Graham 1993). Because some Antrodiaetus species have evolved to live in lower-elevation, dry, xeric habitats (i.e., A. pugnax and A. hageni) compared with other more montane, mesic lineages (e.g., A. cerberus, “pacificus EAST), it is likely that the transitioning ecosystems both fragmented ancestral populations and presented new opportunities for niche adaptation.

The formation of the North American prairies and grasslands was a complex evolutionary event that included the xerification of central North America caused by the rain shadows due to the orogeny of the Rockies. This event likely contributed to the initial splitting of the large species groups within Antrodiaetus. These grasslands in the Nearctic are believed to have started forming during the Eocene, with modern prairies only becoming recognizable much later. This grassland expansion also probably explains the odd placement of A. robustus within the Inland PNW group. The common ancestor of that species with A. cerberus and A. pugnax likely fragmented into the extant species as the central plains expanded and created a large biogeographic barrier. Finally, this changing landscape also played a role in diversification within the Lincolnianus group itself. Consistent with Coyle (1971) and Hendrixson and Bond (2007), A. apachecus is not considered a recent split within that group but is considered to have diverged during the mid- to late Miocene as forest and savanna habitats were succeeded by grasslands and prairies (Retallack 1997; Herbert et al. 2016). It is likely, then, that A. apachecus was in the American Southwest before the orogeny of the Madrean Sky Islands—a unique montane ecosystem and biodiversity hotspot—subsequently becoming isolated in these mountains and the Colorado Plateau as the region gradually became drier, leaving it isolated on the “sky islands.”

The separation of “pacificus SE” from “pacificus WEST” occurs during the mid- to late Miocene, depending on the model of divergence, with the best model supporting a date of ~8 million ybp. This date is older than the Pliocene split suggested by the AV hypothesis but within the timeframe of the Cascade orogeny—the primary driver of the biogeographic break across the Columbia Basin. Throughout the Miocene, the PNW was geologically active, and the changing regional climate and habitats could account for an earlier (relative to AV) isolation of these lineages (Wolfe 1969). The Cascade orogeny and subsequent drying of the Columbia Basin and Oregon Highlands would then have reinforced any lineage split that was already occurring, as regional xerification (and global xerification more generally) would have cut off any habitat corridor between suitable forests in the Coastal/Cascade Ranges and the Rocky Mountains (Herbert et al. 2016).

The current distribution of “pacificus SE” is interesting because of its very southern distribution compared with other taxa that have similar biogeographic breaks. The results indicate that this divergence might not fit either hypothesis. If “pacificus SE” has existed in the Rocky Mountains for a significant period of time (as in the AV hypothesis), it is unclear why these spiders did not expand their distribution northward after glacial retreat. Perhaps competition from “pacificus EAST” for shared habitat precluded such expansion. Certainly, more sampling to determine the distribution extents of “pacificus SE” and “pacificus EAST” would help solve this conundrum.