Specimens representing the genus Ectatosticta Simon, 1892 from China, the sister genus to Hypochilus (and the only other hypochilid genus described), were used to root all phylogenies (UCE data from Ramírez et al. 2021). The monophyly of the family, and two constituent genera, is strongly supported by both morphology (Forster et al. 1987; Alberti and Coyle 1991; Catley 1994; Li et al. 2021) and phylogenomics (Fernández et al. 2018; Ramírez et al. 2021). Phylogenomic data were gathered for forty-three Hypochilus specimens, representing all ten described species (see Suppl. material 1; Fig. 1). A priori species identifications were based on male morphology in combination with geography (Forster et al. 1987; Catley 1994), and we included many samples from prior genetic research (Keith and Hedin 2012). Multiple specimens per species were sampled, chosen to maximize the breadth of geographic coverage within species (Fig. 1). Our southern Rocky Mountain samples included collections from near the respective type localities (Gertsch 1964; Catley 1994) of both H. bonneti Gertsch, 1964 and H. jemez Catley, 1994, plus two geographically intermediate populations from southern Colorado of uncertain species identity. For species delimitation focal taxa (H. pococki and H. petrunkevitchi), larger sample sizes were used to obtain more fine-scale genetic data on potentially cryptic species (Fig. 1). For H. pococki, sampling included representatives of the five mitochondrial haplogroups identified in Keith and Hedin (2012). For H. petrunkevitchi, specimens were sampled from multiple drainage basins, including the Merced River (YOSE), San Joaquin River (SAN), Kings River (KING), Kaweah River (KAW), Tule River (TULE), and Cedar Creek (CEDAR) drainages.

For almost all specimens (with tissues stored at -80 °C), DNA extraction was performed using a Qiagen DNEasy kit from leg tissues. Sequence capture libraries were prepared using an ultraconserved elements capture protocol for arachnids (Starrett et al. 2017; Hedin et al. 2019), with arachnid probes designed by Faircloth (2017). Sequencing was done at the Brigham Young University DNA Sequencing Center on an Illumina HiSeq 2500 150 cycle paired-end sequencing platform. Published data for two H. pococki specimens (H595 and H232, from Starrett et al. 2017), one H. kastoni Platnick, 1987 (G2519, Hedin et al. 2019), and Ectatosticta (Ramírez et al. 2021) were used from prior studies. Raw reads were filtered using the illumiprocessor wrapper (Faircloth 2013) within PHYLUCE v1.6 (Faircloth 2016), after which cleaned reads were assembled using Trinity v2.0.6 (Grabherr et al. 2011) and Velvet v1.0.19 (Zerbino and Birney 2008) on the HPC Cluster at UC Riverside. These assemblies were combined and resulting contigs were matched to probes with minimum identity and minimum coverage values of 80 (--min-identity 80 --min-coverage 80). UCE loci were aligned and trimmed within PHYLUCE using MAFFT (Katoh and Standley 2013) and Gblocks (Castresana 2000; Talavera and Castresana 2007) with relatively liberal settings for GBLOCKS (--b1 0.5 --b2 0.5 --b3 10 --b4 8).

A 50 percent occupancy matrix (623 loci) was generated from the pipeline above (here called the “unfiltered” matrix). A second matrix was further filtered to remove duplicate and potentially non-homologous sequences. Previous work has shown that arthropod UCEs are mostly located within exonic regions (Bossert and Danforth 2018), and the arachnid probe set is no exception (Hedin et al. 2019). In fact, the arachnid probe set can target separate exons from the same protein as separate loci (Hedin et al. 2019), perhaps violating assumptions such as independence and linkage of loci. An annotated list of the arachnid UCEs from Hedin et al. (2019) was used to identify duplicate loci. A total of 73 duplicate loci was found, with the longest locus in each instance retained while the rest were discarded. Using the annotated list, those loci identified by Hedin et al. (2019) as including potential paralogs were also identified and discarded (n = 2), leaving a matrix containing a total of 550 loci. Finally, this filtered dataset was trimmed again with more stringent Gblocks settings (--b1 0.5 --b2 0.85 --b3 4 --b4 8), resulting in a “filtered and trimmed” matrix. This last step was conducted to isolate as much of the purely exonic region as possible for each locus.

Interspecific relationships were reconstructed using all three UCE data matrices (“unfiltered”, “filtered”, and “filtered and trimmed”), utilizing both concatenation and coalescent-model approaches. Data were partitioned by locus, with optimal models selected using PartitionFinder2 (Lanfear et al. 2012) on the CIPRES portal HPC Cluster. Concatenated maximum likelihood trees were reconstructed using RAxML v8.2.12 (Stamatakis 2014) and IQ-TREE v1.5.0 (Nguyen et al. 2015). Using IQ-TREE 2 we also calculated gene (gCF) and site concordance (sCF) factors. For every node of a reference tree, gCF can be defined as the percentage of “decisive” gene trees containing that node, while sCF can be defined as the percentage of decisive sites (in an alignment) supporting a node (Minh et al. 2018). The latter support metric is particularly useful when individual gene trees are uncertain, perhaps because individual alignments are short. Concordance factor calculations were performed using the topology from IQ-TREE which has an identical topology to the RAxML tree and mostly agrees with the SVDquartets reconstruction (see Results). The latter is a coalescent-model topology reconstructed using SVDquartets v1.0 (Chifman and Kubatko 2014) implemented in PAUP* (Swofford 2003), set for 1M quartets with 500 bootstrap replicates for all runs.

Nuclear single nucleotide polymorphism (SNP) data were extracted from UCE loci following the methods of Harvey et al. (2016) and a combination of tools and methods from vcf tools (Danecek et al. 2011) and a modified version of the best practices approach for variant isolation with GATK v4.0.0.0 (Van der Auwera et al. 2013). Separate datasets were created for Californian H. petrunkevitchi plus H. bernardino specimens, and for samples of H. pococki. These matrices started with cleaned read data containing only relevant samples and used a highest coverage reference specimen (H_petrunkevitchi_G2543, H_pococki_H551). Genetic structure and sample clustering were explored using k-means clustering in STRUCTURE v2.3.4 (Pritchard et al. 2000), with unlinked SNPs. Multiple K values (K = 1–10) were run for 100,000 generations with each K value replicated 10 times. Optimal K values were determined following Pritchard et al. (2000) and Evanno et al. (2005), using the online resource CLUMPAK (Kopelman et al. 2015).

Using multiple data sources (phylogenomic results, STRUCTURE results, geography for H. petrunkevitchi, and mitochondrial haplogroup membership for H. pococki), alternative species models were generated and compared using nuclear SNP datasets in the program SNAPP (Tables 1, 2; Bryant et al. 2012). These analyses also included outgroup data (other Hypochilus species), allowing us to test hypotheses of H. pococki as a single species (current taxonomy) and H. petrunkevitchi lumped with H. bernardino (see Tables 1, 2). The averaged marginal likelihoods of duplicate runs were compared for alternative species models using Bayes Factor Delimitation for genomic data (*BFD) (Leaché et al. 2014); we followed the recommendations of Kass and Raftery (1995) in interpreting Bayes factor values.

A rooted triplet species delimitation approach was also implemented using the Python2 compatible version of the program TR2 (Fujisawa et al. 2016). Here, nuclear gene trees are decomposed into partially rooted triplets and congruence is assessed among triplet topologies using a likelihood model testing framework. Input gene trees were constructed in RAxML using rapid bootstrap analysis (-f a) with 200 bootstrap replicates for each gene tree with 550 UCE loci from the “filtered and trimmed” dataset. With the intent to detect patterns of increasing support for increasing species number (i.e., over-splitting), models in which every tip was categorized as a species were included in both SNAPP and TR2 runs.

Mitochondrial Cytochrome c oxidase subunit I (CO1) data for the California taxa were captured from UCE “by-catch” for purposes of phylogenetic and distance analyses, particularly considering the extreme CO1 distances observed in Hedin (2001). A consensus reference sequence was created from specimens of H. petrunkevitchi and H. bernardino (from Hedin 2001), which was then used as a custom database for a BLASTN search for extracting CO1 sequences from UCE contigs. The BLASTN search and subsequent alignment was performed using Geneious Prime (2020.2). Sequence alignments were partitioned by codon position using MODELFINDER (Kalyaanamoorthy et al. 2017) and a phylogeny was constructed using IQ-TREE with 1000 ultrafast bootstrap replicates. A pairwise mitochondrial distance matrix was generated in PAUP* (Swofford 2003) using a Kimura two-parameter model of nucleotide substitution (Kimura 1980).

A Generalized Mixed Yule Coalescent (GMYC) model was used to algorithmically delimit Californian species using CO1 data, in a manner similar to the approach of Keith and Hedin (2012) for Appalachian H. pococki. A GMYC model assumes a difference in intra and inter-specific branching patterns under maximum likelihood and estimates a threshold for the transition from intraspecific (a coalescent process) to interspecific (speciation) branching on an ultrametric tree. Both a single threshold and multiple thresholds can be estimated in different approaches with the model. To this end, an ultrametric CO1 tree was generated using the chronopl function from the APE library in R version 4.0. 2 (Paradis and Schliep 2019), with both single and multiple threshold models performed on the GMYC web server (https://species.h-its.org/gmyc/; Fujisawa and Barraclough 2013).

Although Hypochilus is a strongly morphologically conserved genus, current species were described and are diagnosed using subtle morphological variation, mostly in male genitalia (Hoffman 1963; Forster et al. 1987; Catley 1994). Using SEM we examined male palpal morphology for California taxa (H. kastoni, H. bernardino, multiple lineages of H. petrunkevitchi), and for all primary lineages of H. pococki. We lacked adult males for the southern Rocky Mountain populations of uncertain species identity, so could not study them at this time. Male pedipalps preserved in 80% EtOH were transferred to pure EtOH (200 PF, ≥ 99.5%) for at least 10 min. Following EtOH dehydration, samples were placed in a multilayered sample holder and dried to the critical point using a Tousimis SAMDRI790. Transitional and intermediate fluids used were CO₂(l) and pure EtOH, respectively. Dried palps were mounted on aluminum stubs fitted with sticky carbon conductive spectro-grade tabs. Following mounting, samples were run through an EMS Quorum Q150T sputter coater and covered with 6 nm platinum nanoparticles. SEM micrographs of retrolateral and prolateral views were taken using a vertical stage on an FEI Quanta 450 FEG. Female spermathecal organs were imaged using a Visionary Digital Imaging System, comprising a Canon EOS 5D Mk II DSLR mounted to an Infinity Optics microscope tube. Spermathecal organs were dissected from specimens using fine forceps, immersed for 2–5 minutes in BioQuip specimen clearing fluid (http://www.bioquip.com), then imaged in this fluid on depression slides.

For taxonomic descriptions, morphological measurement details follow Catley (1994: figs 1–4):

PTW/PTL maximum width of male pedipalpal tibia in retrolateral view/length of tibia in retrolateral view;

CdL male palpal conductor length in retrolateral view;

AME diameter of anterior median eye pupil;

PTaL length of male palpal tarsus in retrolateral view;

CTpr number of promarginal cheliceral teeth;

CTre number of retromarginal cheliceral teeth.

Measurements were taken from alcohol-preserved specimens using an Olympus SZ40 dissecting microscope fitted with an ocular micrometer, and converted to millimeters; raw measurements are provided in Suppl. material 2, and summarized in Table 5.

Original UCE raw reads have been submitted to the SRA (PRJNA760946), with summary statistics presented in the Suppl. material 1. Data matrices and resulting .tre files are deposited at Dryad (https://doi.org/10.5061/dryad.g1jwstqsd). The “unfiltered” matrix included 623 UCE loci with an average length of 783 base pairs (bp) and a concatenated length of 734,881 bp (189,387 parsimony informative (PI) sites). The “filtered” matrix included 550 loci with an average length 921 bp and concatenated length of 506,689 bp (145,037 PI sites), while the “filtered and trimmed” matrix contained 550 loci with an average length of 591 bp and concatenated length of 325,452 bp (81,911 PI).

Concatenated ML and SVDquartets analyses of the above three matrices recover nearly identical Hypochilus relationships, except for some nodes in the Appalachian and Rocky Mountain clades (Figs 2, 3; https://doi.org/10.5061/dryad.g1jwstqsd). These analyses confirm regional faunas as monophyletic, rejecting CA paraphyly, and recover Appalachia and California clades as sister taxa. Within Appalachia, a sister species relationship between the geographically disjunct H. thorelli Marx, 1888 and H. coylei Platnick, 1987 is strongly supported across all analyses, consistent with mitochondrial evidence (Hedin 2001; Keith and Hedin 2012), but contrary to morphological evidence which groups H. coylei with the geographically adjacent H. sheari Platnick, 1987 (Huff and Coyle 1992; Catley 1994). The monophyly of all currently recognized Appalachian species is supported, contradicting mitochondrial results which had previously suggested H. pococki paraphyly. Certain parts of the Appalachian topology include low bootstrap values, low gene and site concordance values, and discordant topologies among analyses. Whereas ML analyses place H. pococki as sister to other Appalachian species, SVDquartets nests H. pococki well within the clade with H. gertschi Hoffman, 1963 as sister to other Appalachian species (Fig. 2). A California clade is consistently recovered, with H. kastoni sister to remaining lineages. Hypochilus bernardino is nested prominently within H. petrunkevitchi, and for this reason analyses examining species boundaries included samples of all three species (see below).

Figure 2. 

UCE phylogeny. Phylogeny reconstructed from partitioned RAxML analysis of “filtered and trimmed” UCE matrix. Bootstrap support values are 100 for all nodes unless otherwise indicated; second support values from IQ-TREE. Inset – SVDquartets UCE tree with bootstrap support values less than 95 shown.

Figure 3. 

Concordance Factor values. Phylogeny reconstructed from partitioned RAxML analysis of “filtered and trimmed” UCE matrix, with gCF / sCF values. The lowest gCF and sCF values, for two unresolved nodes which collapse to a four lineage polytomy within the Appalachian clade, are highlighted by red text.

Nuclear SNP datasets included 670 unlinked SNPs for the California sample (allowing 20% missing data) and 655 unlinked SNPs for the Appalachian sample (19% missing data). Overall, STRUCTURE analyses reveal strong population structure for both samples, with inferred genetic clusters congruent with phylogenomic clades (Fig. 4). Within California there is little evidence for mixed ancestry, with genetic populations of H. petrunkevitchi also appearing to be structured by river basin (Fig. 4). Best K as determined by the Pritchard et al. (2000) method is decisive for K = 5 while the Evanno method is less conclusive, supporting a scenario of K = 2 with almost equal but slightly lower support for K = 4 (https://doi.org/10.5061/dryad.g1jwstqsd). When H. kastoni samples were included in STRUCTURE runs, the Pritchard method recovered the same population structure scheme with the addition of a separate H. kastoni “population”, while the Evanno method recovered stronger support for K = 2 in which H. kastoni and H. bernardino comprised a single genetic population. Because we viewed this latter K = 2 result as spurious (see Janes et al. 2017), we generally preferred optimal K values as inferred by the Pritchard method.

Figure 4. 

UCE STRUCTURE results, with optimal and suboptimal K clusters shown in relation to UCE RAxML phylogeny and geography. For H. pococki, the distribution of mitochondrial “microclades” follows Keith and Hedin (2012).

There is more evidence for mixed ancestry in H. pococki STRUCTURE analyses, using a best K from both Pritchard (K = 4) and Evanno (K = 3) methods. The Elk and Northeast “microclades” (ELK and NE) are lumped as a single genetic population, while one member from the mitochondrial Central (CENT) group clusters with the WEST population (Fig. 4), discordant with previously delineated mitochondrial groups (Keith and Hedin 2012). The geographically isolated BONE population, not sampled in previous mitochondrial studies, clusters with disjunct WEST specimens (Fig. 4).

Alternative species hypothesis models were generated and compared using SNAPP and TR2 for both H. petrunkevitchi (five models) and H. pococki (four models). Pritchard-based best K schemes were used for the STRUCTURE derived species models. The most-favored SNAPP model for H. petrunkevitchi (Table 1) is one where every tip represents a species; the next best supported model is the “Basins”, 8-taxon model. SNAPP results for H. pococki were similar and favored every tip as a species as the best model, with a pattern of more speciose models being more favored (Table 2). The rooted triplet TR2 approach also found this pattern of favoring more species-rich models over current taxonomy with the most species-rich model, every tip as a distinct taxon, being the most favored (Tables 3, 4).

The CO1 by-catch phylogeny, using H. kastoni and H. bernardino as possible outgroups, shows strong support (BP = 100) for a clade including H. bernardino and H. petrunkevitchi together (Fig. 5). Within this clade, recovered mitochondrial lineages are consistent with UCE optimal K = 5 STRUCTURE lineages, but the interrelationships of these mitochondrial lineages are not resolved. Considering this phylogenetic uncertainty (i.e., collapsing poorly-resolved nodes), the mitochondrial results are not strictly inconsistent with nuclear results. Pairwise mitochondrial distance values are extremely high among primary lineages (Fig. 5 inset), ranging from 12%–15%. Divergence values within lineages are lower, except for the combined SAN + KING + KAW lineage (> 12% divergence); this obviously reflects significant mitochondrial divergence across drainage basins within this more broadly-distributed lineage. Similarly, there is evidence for structuring across drainage basins within the combined TULE + CEDAR lineage (Fig. 5 inset). As shown previously for Appalachian H. pococki (Keith and Hedin 2012), implementation of a GMYC model using CO1 data appears to severely over-split species. Specimens from the same geographic location are collapsed as the same species, but all unique geographic locations are delimited as distinct species (multi = 13, single = 14; Fig. 5 inset).

Figure 5. 

Maximum likelihood mitochondrial tree. Black circles designate clusters of sequences collapsed as the same species in multi-threshold GMYC analyses; all other branches supported as separate species (e.g., n = 13 for multi-threshold model). Inset – K2P distances within and among primary mitochondrial lineages.

Morphological data for H. pococki and California taxa are presented below in the Discussion and Taxonomy sections, respectively.

Our results confirm the Catley (1994) hypothesis of a sister relationship between highly disjunct California and Appalachian faunas, sister to a more geographically central Rocky Mountain clade. Support was unequivocal for this topology, recovered in all analyses from all UCE matrices (Figs 2, 3; https://doi.org/10.5061/dryad.g1jwstqsd). These regional relationships in Hypochilus are contrary to patterns seen in the salamander genus Aneides, a taxon which also includes montane-associated species from California, the southern Rocky Mountains, and Appalachia. Molecular phylogenetic research in Aneides has recovered a sister relationship between California and southern Rocky Mountain species, sister to Appalachian taxa (Vieites et al. 2007). Inter-regional divergences in Aneides are estimated to have occurred roughly in the timeframe spanning the Eocene to the Oligocene, perhaps coincident with periods of global warming allowing for intercontinental dispersal events (Vieites et al. 2007).

We hypothesize that Hypochilus has a more ancient history, and that this timing difference might also explain the unique phylogenetic topology seen in Hypochilus. Divergence time estimates for Hypochilus are hindered by a lack of direct fossil evidence, with current age estimates for the genus derived from broader examinations of diversification dates for spiders. For example, using published transcriptome data, Magalhaes et al. (2020) estimated divergence between H. gertschi and H. pococki (both Appalachian taxa) during the early Paleogene, with a very large confidence interval. Given approximately polytomous relationships within Appalachia (see below), this point estimate would correspond to a crown group age for the Appalachian radiation, and thus implies older divergences at the base of Hypochilus, perhaps during the Cretaceous. We note here that many other non-entelegyne araneomorph spider lineages are at least this ancient (both within extant families and sometimes within extant genera), as estimated from molecular clock analyses (e.g., sicariids – Magalhaes et al. 2019; leptonetids – Ledford et al. 2021), but also known directly from Upper Cretaceous Burmese amber fossils (e.g., psilodercids – Magalhaes et al. 2021). Also, the combination of Cretaceous fossil evidence in the context of living spider families suggests that non-entelygyne araneomorph lineages (akin to Hypochilus and Hypochilidae) dominated spider diversity at this time (Wunderlich 2008; Magalhaes et al. 2020).

Hypothesized Cretaceous-age divergences for Hypochilus are complicated by the presence of the Western Interior Seaway of North America, a major transcontinental marine barrier in place from ~ 105–65 mya (Blakey and Ranney 2017). Although an east / west vicariance hypothesis associated with the Western Interior Seaway seems attractive, the Rocky Mountain orogeny took place after the withdrawal of the Western Interior Seaway (Blakey and Ranney 2017). We speculate that Hypochilus is either much older than imagined (ages exceeding 105 mya), or that diversification (and dispersal) was spurred soon after the withdrawal of the Western Interior Seaway. Lacking the discovery of relevant fossils, future work could aim to more precisely estimate rates of nuclear gene molecular evolution in Hypochilus, in order to better understand the origin and timing of Hypochilus diversification events. Also, inclusion of a transcriptome representing the Rocky Mountain clade could be incorporated into the well-calibrated Magalhaes et al. (2020) dataset, allowing for a crown group age estimate for the entire genus.

Within Appalachia, nuclear UCE data strongly support a monophyletic H. pococki, contra mitochondrial paraphyly as in Keith and Hedin (2012). Although all currently described species in Appalachia are recovered as monophyletic, and the geographically disjunct H. thorelli and H. coylei are strongly supported as sister species (see also Hedin 2001; Keith and Hedin 2012), our nuclear datasets otherwise do not resolve species relationships, with an overall topology consistent with a four-lineage polytomy. Gene and site concordance factor values at two key unresolved interspecific nodes take lower values than seen anywhere else in Hypochilus, including all nodes within species (Fig. 3). Gene CF values are 9–13 for these nodes, meaning that only ~ 10% of the UCE alignments support these nodes. Site CF values that hover around minimum values (30%) illustrate that the data are essentially equivocal regarding three possible resolutions of a quartet for both of these unresolved nodes (Lanfear 2018). This incongruence and lack of resolution possibly points to a non-adaptive radiation where lineages became rapidly isolated from one another due to environmental factors, but the nature of incongruence requires more study.

Nuclear STRUCTURE results confirm distinct genetic groups within H. pococki (Fig. 4); however, these genetic groups do not correspond exactly to the previously described mitochondrial “microclades”. In particular, the Alarka Mountain specimen from the mitochondrial Central (CENT) clade of Keith and Hedin (2012) instead groups with the nuclear WEST genetic cluster (Fig. 4). This is important because the geographic boundaries of previously defined mitochondrial clades were hypothesized to coincide with riverine barriers (e.g., the Little Tennessee River separating the CENT versus WEST mitochondrial clades, etc.; see Keith and Hedin 2012: fig. 2). In fact, most phylogeographic studies in the southern Appalachians have primarily relied upon mitochondrial evidence to define geographic groupings (e.g., references in Keith and Hedin 2012). Future studies that include dense geographic sampling of nuclear lineages will be important here, with H. pococki representing a prime candidate. More generally, if gene flow across cryptic lineages is promoting mitonuclear discordance, this system might provide interesting insight into how cryptic lineages interact at areas of contact. Also, areas of parapatric contact can be used to understand the degree to which gene flow is restricted across cryptic lineage boundaries, providing strong and direct tests of species status (see Singhal et al. 2018).

Although SNAPP and TR2 show higher support for increasing species numbers within H. pococki (i.e., a many species hypothesis), nuclear STRUCTURE results and consideration of male pedipalp morphology suggest more conservative species numbers. In their diagnosis, Forster et al. (1987) stated that H. pococki males “can be recognized by the flaplike tip of the palpal conductor”. We thus focused particular attention on this structure when searching for morphological differences that might distinguish primary genomic lineages (e.g., VA, ELK + NE, CENT, WEST), and included representatives of all such lineages in our SEM surveys. We did not examine female variation (e.g., in spermathecal morphology), but this is another character system to search for morphological differentiation. We observed minimal differences in male palpal morphology across H. pococki populations and genomic lineages (Figs 6–8). One possible difference is the shorter secondary coil of the conductor tip observed in ELK specimens (Fig. 6), but we note that ELK itself is well nested within the primary K = 4 lineages (Fig. 4). The discord between nuclear genomic data (and analytical results) which suggest many species, versus morphology which suggests few to one species, is a conspicuous example of the cryptic species challenge, and also focuses attention on patterns of morphological stasis. Despite high genomic divergences and ample evolutionary time, morphological change in Hypochilus remains conservative. We might expect conserved Hypochilus somatic morphology because of selective constraints on both niche evolution and morphological differentiation, under a model of phylogenetic niche conservatism (Keith and Hedin 2012; Fišer et al. 2018). The fact that we also observe similar conservatism in genitalia, where at least genetic drift in isolated populations is expected, is compelling.

Figure 6. 

H. pococki male palp (conductor) comparison. For each specimen, left panel = prolateral view, right panel = retrolateral view. NE lineage A, B Green Mtn (MCH 01_162) C, D Boone Fork (MCH 01_159); ELK lineage E, F Elk River (MCH 01_155) G, H Linville Gorge (MCH 01_165). Shorter secondary coil for ELK specimens highlighted by arrows; VA lineage I, J Cliff Mtn (MCH 04_028) K, L Guest River (MCH 04_027); CENT lineage M, N Hickory (MCH 01_144) O, P Wagon Road Gap (MCH 01_181); WEST lineage Q, R Alarka (MCH 02_168) S, T Starr Mtn (MCH 02_156) U, V Backbone Rock (MCH 04_025) W, X Chunky Gal Mtn (MCH 02_142). Detailed specimen information provided in Suppl. material 2.

Figure 7. 

H. pococki male palp comparison, prolateral views. NE lineage A Green Mtn (MCH 01_162) B Boone Fork (MCH 01_159); ELK lineage C Elk River (MCH 01_155) D Linville Gorge (MCH 01_165); VA lineage E Cliff Mtn (MCH 04_028) F Guest River (MCH 04_027); CENT lineage G Hickory (MCH 01_144) H Wagon Road Gap (MCH 01_181); WEST lineage I Alarka (MCH 02_168) J Starr Mtn (MCH 02_156) K Backbone Rock (MCH 04_025) L Chunky Gal Mtn (MCH 02_142). Detailed specimen information provided in Suppl. material 2.

Figure 8. 

H. pococki male palp comparison, retrolateral views. NE lineage A Green Mtn (MCH 01_162) B Boone Fork (MCH 01_159); ELK lineage C Elk River (MCH 01_155) D Linville Gorge (MCH 01_165); VA lineage E Cliff Mtn (MCH 04_028) F Guest River (MCH 04_027); CENT lineage G Hickory (MCH 01_144) H Wagon Road Gap (MCH 01_181); WEST lineage I Alarka (MCH 02_168) J Starr Mtn (MCH 02_156) K Backbone Rock (MCH 04_025) L Chunky Gal Mtn (MCH 02_142). Detailed specimen information provided in Suppl. material 2.

Both nuclear and mitochondrial genetic structuring is very prominent in the California region, and our analyses show that this structure generally follows a pattern of relatedness by drainage basin (Figs 4, 5, 9). The observed divergence within H. petrunkevitchi was not surprising as it has previously been noted as having high levels of intraspecific mitochondrial variation (Hedin 2001). As also found for Appalachian samples, both SNAPP and TR2 show a trend of increasing support for models with increasing numbers of species (Tables 1, 3). GMYC analysis of CO1 data similarly delimits all unique geographic locations as distinct species (Fig. 5 inset). We contend that not all local populations can represent unique species, and instead view this as another example of algorithmic over-splitting.

Figure 9. 

Southern Sierra Nevada topography map with genetic and morphological sample locations (see Suppl. material 1 and Suppl. material 2). Geographic gaps and other notable geographic features mentioned in the text are highlighted.

Based on bootstrap values, nuclear data strongly support the hypothesis that H. bernardino is phylogenetically nested within H. petrunkevitchi (bootstrap = 100 for all matrices across all analyses; Figs 2, 3; https://doi.org/10.5061/dryad.g1jwstqsd). Acknowledging that bootstrap values can provide an inflated view of support (Lanfear 2018; Minh et al. 2018), we considered several other concordance factor values for this node (gCF = 34.8, sCF = 48.3, Fig. 3). From a gene (individual UCE locus) perspective, of 155 alignments that could have included the branch of interest (gN = 155), 34.84% or ~ 50 alignments (= gCF) showed the Fig. 3 topology, with very few alignments supporting an alternative topology at high frequency (gDF1 = 3.23, gDF2 = 8.39). Similarly, from a site perspective, of ~ 200 decisive sites for the quartet of interest (sN = 199.42), approximately half support the Fig. 3 topology (sCF = 48.28), with fewer supporting alternative resolutions (sDF1 = 18.33, sDF2 = 29.4). Overall, we view these values (see Lanfear 2018), in concert with bootstrap values, as strongly supporting the paraphyly of H. petrunkevitchi with respect to H. bernardino.

Based on this phylogenomic pattern we see two obvious taxonomic alternatives. The first is to sink H. bernardino into a broadly distributed, highly genetically-structured H. petrunkevitchi. The second is to elevate and formally describe the distinct genetic lineage that is sister to H. bernardino. We prefer and argue for the latter approach, for reasons concisely summarized as follows: 1) H. bernardino is already described based on a diagnostic morphology (Catley 1994, and revised diagnosis below), 2) H. bernardino is geographically-localized, known only from a single mountain range in southern California, highly disjunct from more northern populations (Fig. 9), and 3) H. bernardino is supported both as a distinct nuclear and mitochondrial genetic lineage (Figs 2–5). By these multiple measures of genetic, morphological, and geographic distinctiveness, we view H. bernardino as an evolutionary lineage on a unique and independent trajectory. Catley (1994) described H. bernardino from the southern section of the San Bernardino mountains of southern California. The few known populations are separated from the southern Sierra Nevada by hundreds of kilometers of mostly inappropriate habitat (lower elevations, fewer granite outcrops), where these spiders (or their distinctive webs) have never been found (Fig. 9). This lack of records includes not only our own extensive field work in the intervening region, but also that of the many arachnologists that have conducted research in California, as well as modern-day tools such as iNaturalist. This sort of geographic disjunction has been found in some other taxa, for example, the iconic Ensatina salamander complex, where the geographic disjunction is known as “Bob’s Gaps” (Jackman and Wake 1994). In this instance, separated populations have been described as separate subspecies (E. eschscholtzii klauberi in the Tehachapi mountains distinct from E. eschscholtzii croceater in the San Bernardino mountains; Jackman and Wake 1994), but Ensatina taxonomy is generally regarded as being fairly conservative.

Accepting H. bernardino as an independently evolving lineage, the current taxonomy must be updated to reflect unique lineages discovered within H. petrunkevitchi. Conservatively, we retain the northern lineages that include the type locality (male holotype of H. petrunkevitchi from Cedar Grove, Fresno County = KINGS lineage) as H. petrunkevitchi. This lineage is distributed across the Merced, San Joaquin, Kings and Kaweah River basins (Fig. 9). Again, accepting H. bernardino as a unique species, our data indicate that populations from the Tule River and Cedar Creek drainage basins need to be recognized as a new species, which we formally describe below. Specimens from these drainage basins represent new locality records and the southern-most known observations of Hypochilus spiders in the Sierra Nevada (Fig. 9). More generally, this part of the southern Sierra Nevada is a well-known area of active speciation, with many short-range endemic arthropods and vertebrates (Bond 2012; Jockusch et al. 2012; Papenfuss and Parham 2013; Satler et al. 2013; Leavitt et al. 2015; Emata and Hedin 2016; Starrett et al. 2018; Bennett et al. 2021; Weng et al. 2021). In this regard, discovering a new species in the southern Sierra Nevada is not surprising.

We define species as evolutionary lineages on a “unique and independent trajectory”. Following Davis et al. (2021), we consider species to be cryptic “if they depend on additional sources of data to formulate the delineation hypothesis prior to establishing diagnostic morphological characters”. This definition applies well to the new species description below, which is motivated by a combination of genetic divergence and uniqueness, phylogenetic pattern (paraphyly w.r.t. H. bernardino) and geographic allopatry, which has prompted us to take a closer look at morphology. Additionally, the definition allows for an initial hypothesis of morphological crypsis that does not preclude future downstream studies that in fact find morphological (or other) differences, making the species no longer technically “cryptic”. In many instances, species are morphologically cryptic because of their youth, where morphological divergence has not yet caught up with molecular divergence. But as noted above, the measured mitochondrial differences in Hypochilus are among the highest known in spiders (a clade with ~ 50,000 described species), and we hypothesize relatively ancient divergences in the genus. These data and arguments are inconsistent with recent evolutionary divergence, and we instead favor a model of long-term phylogenetic niche conservatism constraining morphological evolution (Fišer et al. 2018), as argued above.

Overall, we view our new cryptic species hypotheses as conservative, consistent with perspectives that genomic data should be interpreted conservatively when describing new species (Coates et al. 2018). In Hypochilus, this is particularly true since genomic diagnosability is ubiquitous, and extends to the level of localized populations, as reflected in nuclear and mitochondrial algorithmic delimitations. Our hypotheses below do not recognize all genetically divergent lineages as species and allow for varying degrees of population divergence within described species (e.g., H. petrunkevitchi, new species below). In particular, there is evidence that the Yosemite Valley (Merced River) population is on a unique evolutionary trajectory due to its disjunct distribution and because samples from the isolated valley floor routinely fall out as a divergent group in both nuclear and mitochondrial analyses (Figs 2–5). The single sampled population from the San Joaquin drainage is similarly genetically divergent and geographically isolated. This disjunction might be natural, as spiders have never been collected in the apparently suitable granite outcrop habitats between San Joaquin locations and the Yosemite Valley, despite concerted collecting efforts (Fig. 9). In the context of an integrative taxonomic framework, we weigh conservation considerations, extreme geographic allopatry, and well-supported paraphyly as particularly important. Genomic diagnosability is important but not decisive, with morphological diagnosability as least important, again reflecting phylogenetic niche conservatism. Regarding conservation value, we argue that sinking H. bernardino into a broadly distributed, highly genetically-structured H. petrunkevitchi would immediately decrease the conservation value and importance of the former.