Abstract

Species are a fundamental unit for biological studies, yet no uniform guidelines exist for determining species limits in an objective manner. Given the large number of species concepts available, defining species can be both highly subjective and biased. Although morphology has been commonly used to determine species boundaries, the availability and prevalence of genetic data has allowed researchers to use such data to make inferences regarding species limits. Genetic data also have been used in the detection of cryptic species, where other lines of evidence (morphology in particular) may underestimate species diversity. In this study, we investigate species limits in a complex of morphologically conserved trapdoor spiders (Mygalomorphae, Antrodiaetidae, Aliatypus) from California. Multiple approaches were used to determine species boundaries in this highly genetically fragmented group, including both multilocus discovery and validation approaches (plus a chimeric approach). Additionally, we introduce a novel tree-based discovery approach using species trees. Results suggest that this complex includes multiple cryptic species, with two groupings consistently recovered across analyses. Due to incongruence across analyses for the remaining samples, we take a conservative approach and recognize a three species complex, and formally describe two new species (Aliatypus roxxiae, sp. nov. and Aliatypus starretti, sp. nov.). This study helps to clarify species limits in a genetically fragmented group and provides a framework for identifying and defining the cryptic lineage diversity that prevails in many organismal groups. [California; cryptic species; multispecies coalescent model; mygalomorph spiders; species delimitation.]

Species are fundamental units of study in evolutionary biology, ecology, and conservation research. Inaccurate understanding of species diversity may lead to errors in analyses that use species as units (e.g., phylogenetic community structure analyses), and may hinder conservation efforts (Wiens 2007). However, the question of what defines a species is contentious, as evident by the large number of species concepts invoked by evolutionary biologists (summarized in de Queiroz 2007). Although multiple types of data can be used for species delimitation (e.g., morphology, ecology, behavior, etc.), multilocus genetic data are increasingly easy to gather for nonmodel taxa. Additionally, genetic data provide potential evidence of divergence at an early stage of diversification; this evidence is contained within the shared ancestral polymorphisms and allows species to be delimited before individual loci acquire monophyly (Knowles and Carstens 2007). This approach to species delimitation is founded upon coalescent theory (Kingman 1982; Hudson 1991), where probabilistic expectations for the sorting of alleles can be generated, and proceeds under the general lineage concept, which defines species as independently evolving lineages (see de Queiroz 2007). Coalescent-based methods for species delimitation model lineage divergence using species tree inference (e.g., Edwards 2009), and several studies have delimited species using this approach (Carstens and Dewey 2010; Leaché and Fujita 2010; Kubatko et al. 2011; Niemiller et al. 2012; Carstens and Satler 2013). Notably, many of these studies have investigated systems where preexisting taxonomic divisions were present (e.g., described subspecies or allopatric populations). Species delimitation using species trees is more difficult when empirical systems lack clear preexisting divisions (see O'Meara 2010), due to the immensity of tree space with an unspecified number of species.

The distinction between species delimitation with and without a priori groupings is particularly important when investigating lesser-known groups such as arthropods, where the potential lack of morphological differentiation precludes a priori hypotheses of species groupings. In fact, evidence from genetic data has led to the detection of cryptic lineage diversity across the tree of life (Sáez and Lozano 2009; Bickford et al. 2007), with multiple approaches available for objectively defining species boundaries (Fujita et al. 2012). However, in such taxa where a priori knowledge about potential species groupings is not available, species discovery approaches that proceed without the prior partitioning of data are likely to be useful (see Ence and Carstens 2011). Another approach would utilize independent lines of evidence to generate a priori groupings necessary for validation approaches. Although simulation studies have been conducted for a variety of discovery (O'Meara 2010; Rittmeyer and Austin 2012) or validation (Ence and Carstens 2011; Camargo et al. 2012) methods, it is difficult to predict which of these general approaches will be most applicable in a given empirical system. Here we explore this question in a species complex of mygalomorph spiders from California.

The spider infraorder Mygalomorphae contains over 2600 described species (Raven 1985; Hedin and Bond 2006; Bond et al. 2012), and includes such taxa as tarantulas, trapdoor spiders, purseweb spiders and kin. The majority of mygalomorph species are ground-dwelling spiders that live in subterranean burrows; females rarely leave their burrows, whereas males generally only leave burrows in search of a mate after reaching adulthood. Male dispersal is likely limited to short distances, and “ballooning,” a long-range dispersal technique used by many immature “true” spiders (Infraorder Araneomorphae), has rarely been observed in mygalomorph species (Coyle 1983). Due to their burrow fidelity and limited dispersal abilities, mygalomorph populations are typically clustered in isolated aggregations (Bond et al. 2006). As predicted for any species with such life history traits, species-level investigations of mygalomorphs typically reveal extensive population genetic structure (Starrett and Hedin 2007; Bond and Stockman 2008; Hedin et al. 2013). Mygalomorph spiders also tend to be morphologically conserved at shallow phylogenetic levels, leading several authors to recognize previously undiscovered species diversity in this clade (Bond et al. 2001; Arnedo and Ferrandez 2007; Hendrixson and Bond 2007; Starrett and Hedin 2007; Satler et al. 2011; Hedin et al. 2013). Cryptic species in mygalomorphs have generally been discovered through the use of a single gene marker (e.g., mitochondrial DNA), or through the concatenation of a handful of gene regions (Hendrixson and Bond 2005; Starrett and Hedin 2007; Satler et al. 2011; Hendrixson et al. 2013). In addition, molecular data have been combined with ecological niche modeling for species delimitation, identifying lineages that are both genetically and ecological distinct (Stockman and Bond 2007; Bond and Stockman 2008; Hendrixson et al. 2013).

To date, studies assessing species limits in mygalomophs have not adopted species tree-based approaches for species delimitation. As such, many previous investigations into mygalomorphs have relied on gene tree monophyly as a principle criterion for assessing species limits, but this criterion can be problematic. Extensive population-level structure leads to the presence of a high number of “micro-clades” in mtDNA gene trees, which could result in splitting each nominal mygalomorph species into dozens of new species that are both geographically distinct and contain diagnostic haplotypes (but see Bond and Stockman 2008). Whereas our work is focused on delimitation within mygalomorph spiders, similar issues are possible with any dispersal-limited taxon where single locus studies reveal extensive population-level structure. A potential solution is to collect data from multiple autosomal nuclear markers, which should not be biased by sex-biased gene flow (Brito and Edwards 2009), generally have a lower mutation rate, and on average retain ancestral polymorphisms longer than mitochondrial genes.

The California Floristic Providence is a biodiversity hotspot (Myers et al. 2000; Calsbeek et al. 2003), and is home to a high diversity of trapdoor spiders (e.g., Starrett and Hedin 2007; Stockman and Bond 2007; Bond and Stockman 2008, Satler et al. 2011; Bond 2012; Hedin et al. 2013). The genus Aliatypus (Mygalomorphae, Antrodiaetidae) includes 12 described species, 11 of which are endemic to California, with one species (A. isolatus) distributed in Arizona (Coyle 1974; Hedin and Carlson 2011; Satler et al. 2011). Aliatypus thompsoni is distributed throughout the Transverse Ranges of southern California, and includes two distributional “arms,” one extending northwest into the south Coast Ranges, and the other extending northeast into the southern Sierra Nevada Mountains (Fig. 1). Through a thorough statistical analysis of morphological traits, Coyle (1974) recognized A. thompsoni under a morphological species concept, identifying multiple diagnostic characters (including a thoracic groove that is either absent or very shallow and posterior sigilla that are large, faint, and close together) unique to this species. However, it was also noted that morphological variation in some traits was present throughout the species' distribution. Recently, Satler et al. (2011) suggested the possibility of cryptic species within A. thompsoni (based upon a phylogenetic species concept, with the consistent recovery of three genetic clades that are both geographically cohesive and isolated). This hypothesis is consistent with the reputation of the region as a known diversification hotspot reflected in areas of endemism and phylogeographic breaks (Jockusch and Wake 2002; Chatzimanolis and Caterino 2007; Davis et al. 2008; Parham and Papenfuss 2009; Polihronakis and Caterino 2010; Jockusch et al. 2012). Here we analyze data from 6 gene regions to test for cryptic species-level lineage diversity within A. thompsoni. We apply several methods for species delimitation, including both discovery and validation approaches, as well as a chimeric approach following Leaché and Fujita (2010). In addition, we introduce a novel species tree-based discovery approach.

Figure 1

Map of southern California showing A. thompsoni sampling localities, with province designations [following Bailey and Jahns (1954] and Buwalda (1954)) highlighted with letters and regions highlighted with shading. Localities used in multilocus species delimitation analyses highlighted with white circles (see Table 1). Detailed collection information for all localities can be found in Supplemental Table 1.

Figure 1

Map of southern California showing A. thompsoni sampling localities, with province designations [following Bailey and Jahns (1954] and Buwalda (1954)) highlighted with letters and regions highlighted with shading. Localities used in multilocus species delimitation analyses highlighted with white circles (see Table 1). Detailed collection information for all localities can be found in Supplemental Table 1.

MATERIALS AND METHODS

Taxon Sampling

Specimens were collected from 83 localities throughout the known geographic distribution of A. thompsoni. In addition, our sampling includes a substantial range extension to the northwest in the south Coast Ranges (sites 37–42; Fig. 1). A sampling gap is present in the Tehachapi Mountains, as we were unable to collect spiders from the large and privately owned Tejon Ranch (between d and e; Fig. 1). At each locality we attempted to collect between 2 and 4 adult spiders, although some collecting sites included either a single individual or only immature spiders. Dense geographic sampling throughout the species distribution, combined with essentially no genetic variation within sampled localities (see Results), informed our decision to only collect a handful of individuals per sampling site. Adult and immature spiders were identified as A. thompsoni based on somatic morphology following Coyle (1974). Specimens in this study have been assigned a unique specimen number (MY or GMY; see Supplemental Table 1, available from the Dryad data repository, doi:10.5061/dryad.68s40). Following phylogenetic hypotheses of Coyle (1994) and Satler et al. (2011), sequences from out-group taxa were used to root phylogenetic trees.

TABLE 1

Sampling locality information

Site Site acronym Region Province Locality 
EGlenville SN Sierra (a) CA: Kern Co., E Glenville 
PosoFlatRd SN Sierra (a) CA: Kern Co., Poso Flat Rd. 
SierraWay SN Sierra (a) CA: Kern Co., Sierra Way 
LakeIsabella SN Sierra (a) CA: Kern Co., Lake Isabella 
10 Highway178 SN Breckenridge (b) CA: Kern Co., Highway 178 
14 RollingOaksRd SN Piute (c) CA: Kern Co., Rolling Oaks Rd. 
16 CalCreekRd SN Piute (c) CA: Kern Co., Cal Creek Road 
17 WoodfordTehRd SN Tehachapi (d) CA: Kern Co., Woodford Tehachapi Rd. 
19 TehachapiMtnParkCG SN Tehachapi (d) CA: Kern Co, Tehachapi Mtn Park CG 
20 ClearCreekRd SN Tehachapi (d) CA: Kern Co., Clear Creek Rd. 
21 ComanchePointRd SN Tehachapi (d) CA: Kern Co., Comanche Point Rd. 
22 DigierRd TR San Emigdio (e) CA: Kern Co., Digier Rd. 
23 LebecOaksRd TR San Emigdio (e) CA: Kern Co., Lebec Oaks Rd. 
25 Highway5 TR Frazier (f) CA: Kern Co., Highway 5 
26 FrazierPark TR Frazier (f) CA: Kern Co., ~ Frazier Park 
27 LakeviewDr TR Frazier (f) CA: Kern Co., Lakeview Dr. 
29 MtPinosRd TR Frazier (f) CA: Kern Co., Mt. Pinos Rd. 
30 CuddyValleyRd TR Frazier (f) CA: Kern Co., Cuddy Valley Rd. 
36 SantaBarbaraCynRd TR San Rafael (g) CA: Santa Barbara Co., Santa Barbara Canyon Rd. 
44 ParadiseRd TR Santa Barbara (h) CA: Santa Barbara Co., Paradise Rd. 
51 LakeCasitas TR Santa Barbara (h) CA: Ventura Co., Lake Casitas (E side) 
59 GrimesCynRd TR LA Basin (j) CA: Ventura Co., Grimes Canyon Rd. 
74 NESantaClarita TR Sierra Pelona (l) CA: Los Angeles Co., NE Santa Clarita 
75 TujungaCynRd TR San Gabriel (m) CA: Los Angeles Co., Tujunga Canyon Rd. 
78 EatonCyn TR San Gabriel (m) CA: Los Angeles Co., Eaton Cyn. 
37 OldSierraMadreRd CR South Santa Lucia (o) CA: Santa Barbara Co., Old Sierra Madre Rd. 
39 RedHillRd CR North Santa Lucia (p) CA: San Luis Obispo Co., Red Hill Rd. 
Site Site acronym Region Province Locality 
EGlenville SN Sierra (a) CA: Kern Co., E Glenville 
PosoFlatRd SN Sierra (a) CA: Kern Co., Poso Flat Rd. 
SierraWay SN Sierra (a) CA: Kern Co., Sierra Way 
LakeIsabella SN Sierra (a) CA: Kern Co., Lake Isabella 
10 Highway178 SN Breckenridge (b) CA: Kern Co., Highway 178 
14 RollingOaksRd SN Piute (c) CA: Kern Co., Rolling Oaks Rd. 
16 CalCreekRd SN Piute (c) CA: Kern Co., Cal Creek Road 
17 WoodfordTehRd SN Tehachapi (d) CA: Kern Co., Woodford Tehachapi Rd. 
19 TehachapiMtnParkCG SN Tehachapi (d) CA: Kern Co, Tehachapi Mtn Park CG 
20 ClearCreekRd SN Tehachapi (d) CA: Kern Co., Clear Creek Rd. 
21 ComanchePointRd SN Tehachapi (d) CA: Kern Co., Comanche Point Rd. 
22 DigierRd TR San Emigdio (e) CA: Kern Co., Digier Rd. 
23 LebecOaksRd TR San Emigdio (e) CA: Kern Co., Lebec Oaks Rd. 
25 Highway5 TR Frazier (f) CA: Kern Co., Highway 5 
26 FrazierPark TR Frazier (f) CA: Kern Co., ~ Frazier Park 
27 LakeviewDr TR Frazier (f) CA: Kern Co., Lakeview Dr. 
29 MtPinosRd TR Frazier (f) CA: Kern Co., Mt. Pinos Rd. 
30 CuddyValleyRd TR Frazier (f) CA: Kern Co., Cuddy Valley Rd. 
36 SantaBarbaraCynRd TR San Rafael (g) CA: Santa Barbara Co., Santa Barbara Canyon Rd. 
44 ParadiseRd TR Santa Barbara (h) CA: Santa Barbara Co., Paradise Rd. 
51 LakeCasitas TR Santa Barbara (h) CA: Ventura Co., Lake Casitas (E side) 
59 GrimesCynRd TR LA Basin (j) CA: Ventura Co., Grimes Canyon Rd. 
74 NESantaClarita TR Sierra Pelona (l) CA: Los Angeles Co., NE Santa Clarita 
75 TujungaCynRd TR San Gabriel (m) CA: Los Angeles Co., Tujunga Canyon Rd. 
78 EatonCyn TR San Gabriel (m) CA: Los Angeles Co., Eaton Cyn. 
37 OldSierraMadreRd CR South Santa Lucia (o) CA: Santa Barbara Co., Old Sierra Madre Rd. 
39 RedHillRd CR North Santa Lucia (p) CA: San Luis Obispo Co., Red Hill Rd. 

Notes: Information for sampling localities used in multilocus species delimitation analyses. Regions are abbreviated.

Molecular Data

Genomic DNA was extracted from ethanol-preserved leg tissue using a DNeasy kit (Qiagen). Six separate gene fragments were PCR-amplified, including 1 mitochondrial and 5 nuclear genes (cytochrome oxidase I (COI) mtDNA, 28S rRNA, 18S rRNA, EF-1γ nDNA, Mid1 nDNA, Anonymous nDNA). Specimens from all localities were sampled for COI (including multiple specimens from each locality when possible); a subsample of individuals was sequenced for all nuclear genes, with individuals selected to maximize geographical coverage (27 localities; see Fig. 1). PCR and sequencing protocols followed Satler et al. (2011). Marker development details for Mid1 and Anonymous are included in the Supplemental Material. PCR amplicons were sequenced in both directions, with sequences edited and assembled in Sequencher v4.5 (Gene Codes Corporation, MI). Nuclear gene sequences containing multiple heterozygous sites were further analyzed in PHASE v2.1 (Stephens et al. 2001; Stephens and Donnelly 2003) to resolve gametic phase. Allelic combinations phased with 0.90 pp (posterior probability) or higher were retained; standard ambiguity codes were used for heterozygous sites below that threshold. Recombination was tested for with TOPALi v2.5 (Milne et al. 2004; Milne et al. 2009), implementing the DSS (Difference of Sums of Squares) method with default program settings. Matrices for COI mtDNA, EF- nDNA, and Mid1 nDNA were manually aligned in MacClade v4.08a (Maddison and Maddison 2005) based on protein translation; a lack of indels allowed for confidence in the alignments. Length variable data (28S and 18S rRNA) were aligned with the program MAFFT (Katoh et al. 2005) using the G_INS-i alignment algorithm, or by using the MAFFT Geneious Pro plug-in (Drummond et al. 2011) with default settings (Anonymous nDNA). The program Gblocks (Castresana 2000) was used to remove alignment ambiguous regions in the 28S rRNA alignment, reducing the matrix from 1278 to 1052 positions under a “less stringent” criterion (minimum number of sequences for a conserved position and flanking regions: 28; maximum number of contiguous non-conserved positions: 8; minimum length of a block: 5; allowed gap positions: with half).

Individual Gene Tree Analyses

Models of DNA sequence evolution were estimated with jModelTest v0.1.1 (Guindon and Gascuel 2003; Posada 2008) using the Akaike Information Criterion (AIC). Redundant COI sequences from the same locality were combined into haplotypes. Individual gene trees were estimated using maximum likelihood in RAxML v7.2.8 (Stamatakis 2006; Stamatakis et al. 2008). Due to restrictions in model choice in RAxML, gene trees were estimated under a GTRGAMMA model, with COI further partitioned by codon position. Nodal support was assessed with 1000 bootstrap replicates. Gene trees estimated with Bayesian inference (for species tree analysis; see below) and non-RAxML maximum likelihood programs used models estimated with jModelTest.

Tests of Species Limits

Understanding species limits in any given taxonomic group can be both highly complex and subject to researcher bias. Furthermore, identifying cryptic lineage diversity in an objective manner, where morphological differences are either absent or not immediately apparent, can be challenging. Here we provide a framework for the discovery of cryptic lineage diversity, exploring both genetic clustering and tree-based approaches. As the fields of phylogenetics and population genetics merge at this shallow phylogenetic level, both approaches have merit and are expected to provide relevant information. As it is our contention that robust evidence is required to recognize cryptic species, we work under the assumption that if cryptic lineage diversity is present, consistent signal will emerge across multiple species delimitation approaches.

Aliatypus thompsoni was described as a single species based upon a morphological species concept (Coyle 1974), but molecular evidence (see Satler et al. 2011) suggests that this species may be composed of multiple cryptic lineages. Several preliminary analyses were conducted prior to delimitation. First, we conducted an analysis of molecular variance (AMOVA; Excoffier et al. 1992) using the program Geno found within GeneticStudio v131 (Dyer 2009). Genetic data were partitioned in a hierarchical manner, with data partitioned both within three major regions (e.g., south Coast Ranges (CR), Transverse Ranges (TR) and southern Sierra Nevadas (SN); see Fig. 1), and then further partitioned within each geologic province (see below). All 6 loci were used, and null distributions were generated with 10,000 permutations to test for significance. In addition to an AMOVA, we conducted a Mantel Test in Geno on all 6 loci (combined) to test for isolation by distance (IBD). A matrix of the covariate distance among individuals was calculated from locality data (latitude/longitude).

Validation approaches to species delimitation require a priori partitioning of samples. In other studies, these partitions have been defined using existing taxonomy (e.g., Carstens and Dewey 2010), morphology (e.g., Barrett and Freudenstein 2011), and/or geography (e.g., Camargo et al. 2012). Because A. thompsoni lacks taxonomic or morphological divisions, we used geological data to help inform our initial species hypotheses. The mountains of southern to central California constitute a heterogeneous landscape impacted by episodes of mountain orogeny, marine embayments, and the action of multiple fault lines (Hill and Dibblee 1953; Wakabayashi and Sawyer 2001; Hall 2002). The distribution of A. thompsoni partially overlaps three primary mountain ranges, including the south CR, the TR, and the southern SN. Several well-defined topographic and geologic units have been identified within these ranges, defined by both structural features and spatial separation by narrow to broad valleys (Bailey and Jahns 1954; Buwalda 1954; see Fig. 1).

Essentially all population-level studies of mygalomorphs recover high levels of population genetic structure with phylogeographic groupings that are geographically cohesive; these results are attributed to the limited dispersal abilities and strong philopatric tendencies of these spiders (see references above). Based on this information, we hypothesize that spiders found in the same area and within regions with shared geological history are likely to be genetically closely related, and use this criterion to determine our initial species hypotheses. We take a hierarchical approach by first testing the relatedness of a priori groupings recovered within the three primary geological groupings (CR, TR, SN; Table 1), and then test relatedness of the retained lineages among the lineages recovered within these regions. This hierarchical approach was required due to computational limits of the delimitation methods (see below). Briefly, the 27 sampling localities create an excessively large number of permutations (>4.9 × 1019) and thus cannot be exhaustively explored in spedeSTEM (which is limited to 8 putative lineages/4140 permutations) or with Bayesian Phylogenetics & Phylogeography (BPP, which requires ≤ 19 initial lineages).

Two validation approaches were used—spedeSTEM v1.0 (Ence and Carstens 2011) and BPP v2.1 (Rannala and Yang 2003; Yang and Rannala 2010). The former calculates the probability of different models of lineage composition using maximum likelihood and evaluates these models using information theory (Burnham and Anderson 2002). As spedeSTEM requires all terminals to be represented across loci, Anonymous nDNA data were removed (due to missing specimen data) resulting in a 5 locus data set for this analysis. As most individuals were homozygous across all loci, a single allele was sampled per individual, with the removal of redundant alleles improving the calculation of likelihood values (see Supplemental Material for details regarding subsampling). Analyses consisted of gene tree estimation using the tree bisection-reconnection search strategy, DNA substitution models selected in jModelTest, and θ = 4Neμ values estimated using Migrate-n v3.2.7 (Beerli and Felsenstein 1999; Beerli and Felsenstein 2001). Each analysis included 500 replicates.

We also used the Bayesian species delimitation analysis, BPP, which requires an input guide tree representing the species phylogeny with all possible species. BPP evaluates speciation models using a reversal jump Markov Chain Monte Carlo (rjMCMC) algorithm to determine whether to collapse or retain nodes throughout the phylogeny. Since the selection of an accurate guide tree is critical to delimitation accuracy (see Leaché and Fujita 2010), we estimated the input phylogeny using *BEAST v1.6.1 (Heled and Drummond 2010). Multiple analyses with varying priors (θ and τ; see Leaché and Fujita 2010) were run to discern how varying the ancestral population sizes and root ages influenced results. A conservative approach was used for the BPP analysis; we required strong support (pp ≥ 0.95) across all runs to retain a given node (i.e., indicating lineage splitting). Analyses included all 6 loci, with a single allele sampled at random per individual. The settings for BPP were as follows: species delimitation was set to 1, the algorithm was set to 0, and the fine-tune parameter (ϵ) was set to 10. Analyses were run for 500,000 generations (first 10,000 were burn-in), with a sampling interval of 5. Each analysis was repeated twice using different starting seeds to confirm consistency between runs.

The two validation approaches discussed above require operational taxonomic units (OTUs) to be designated a priori, and may mislead delimitation inferences if these groupings are inaccurate. Thus, we also utilized several species discovery approaches, where data are analyzed without assignment into a priori groupings. The first approach involved the use of Structure (Pritchard et al. 2000), a genetic clustering algorithm that attempts to infer population structure through the use of allele frequencies and identifies genetic clusters in both Hardy–Weinberg and linkage equilibrium. Structure runs were conducted assuming between 1 and 27 populations (K = 1 through K = 27), with each K step replicated 5 times. Analyses used an admixture model, were conducted with a burn-in of 1 ×105 steps (with 1 ×106 Markov Chain Monte Carlo (MCMC) steps after burn-in), and allele frequencies were considered independent among populations. The optimal K was selected using the Evanno method (Evanno et al. 2005) estimated in Structure Harvester (Earl and vonHoldt 2011), and the data were summarized with CLUMPP (Jakobsson and Rosenberg 2007) using the FullSearch algorithm and visualized with DISTRUCT (Rosenberg 2004).

Whereas recent simulations suggest that genetic clustering algorithms can produce accurate species delimitations (Rittmeyer and Austin 2012), empirical studies have also observed incongruence between the hierarchical pattern of lineage divergence and the pattern of population clustering across varying levels of K (e.g., Jackson and Austin 2010). Therefore, we utilized additional discovery approaches that co-estimate species limits and relationships. The first approach used was Brownie (O'Meara 2010), a nonparametric method that recovers species groupings by minimizing the amount of incongruence on interspecific branches, while maximizing incongruence on intraspecific branches (i.e., limiting population structure). Brownie is based on the assumption that the most probable gene tree is congruent with the species tree, and aims to quantify the amount of gene tree conflict by minimizing intraspecific genetic structure. Five loci were used for this analysis (excluding Anonymous nDNA due to missing data), and a single allele was randomly selected from individuals representing all 27 localities (see Niemiller et al. 2012). For input, maximum likelihood gene trees (from RAxML) were converted into ultrametric trees using PAUP* v4.0 (Swofford 2002). Analyses were run under default settings using a heuristic search, and replicated 10,000 times. Species limits hypotheses were summarized from the replicates, with probabilities assigned to lineage compositions that were recovered in at least one run (calculated as their frequency over the total number of replicates). Because the parameter space explored by the Brownie heuristic algorithm is much smaller than the number of possible species trees for a given number of tips, it is difficult to assess how well the heuristic search is performing (see O'Meara 2010). Therefore, we developed an additional discovery approach that does not rely on a heuristic search of this parameter space.

Most species delimitation analyses that incorporate species trees compare the probability of trees with fewer or greater numbers of OTUs to identify optimal partitions of the data (e.g., spedeSTEM, BPP). We extended this strategy to its maximum extent by beginning at the level of the individual sample site, and calculated the probability of the phylogeny that treats individual samples as putative lineages. We then collapsed the two samples (or multiple samples if represented as a polytomy) that were most closely related, as measured by multilocus branch length, into a single lineage and recalculated the probability of the new phylogeny. This process of collapse and recalculation was continued until all members of the nominal species were collapsed, and information theory (Burnham and Anderson 2002) was used to identify the optimal model of lineage composition. STEM v2.0 (Kubatko et al. 2009) was used to calculate the probability of the species tree given the gene trees [ST | GT], using the settings and gene trees estimated above. Because STEM computes the maximum likelihood species tree that maximizes the probability of a given set of gene trees under the model (i.e., sample partitioning scheme), our method will necessarily compare the optimal trees for each level of assumed lineage composition. The log likelihoods of the data given the models were evaluated using information theory following Carstens and Dewey (2010). We also repeated this general approach in a Bayesian framework, using the Maximum Clade Credibility (MCC) species tree estimated from *BEAST and all samples, combined with BPP as above.

In addition to the multilocus discovery approaches summarized above, we employed a single locus discovery method using the GMYC model (Pons et al. 2006) applied to the COI data. The GMYC model seeks to identify the point in the phylogeny where coalescent branching events within species transition to those corresponding to species level divergence. However, rather than using the likelihood version of this model, we used a Bayesian implementation (bGMYC; Reid and Carstens 2012) that incorporates gene tree uncertainty by sampling over the posterior distribution of sampled gene trees. For the bGMYC runs, ultrametric gene trees were estimated in BEAST v1.6.1 (Drummond and Rambaut 2007) under a strict clock model using a MCMC run of 5 × 107 generations with sampling every 5 × 103generations. Following a 40% burn-in, the posterior distributions were further trimmed to 100 trees (evenly sampled throughout the posterior), and used as input for bGMYC. Default settings were used, although the parameter py2 was set to 1.2 (recommended by the authors), and the starting number of species was set to half the total number of tips. This analysis was run on the full (all sampled localities) and geographically reduced (localities with all genes sequenced) data sets. For both analyses, bGMYC was run for 20,000 generations, with a burn-in of 10,000 generations, and sampled every 200th generation.

Finally, we followed Leaché and Fujita (2010) in applying a chimeric approach that uses a combination of discovery and validation methods. Rather than using external information to inform our initial species hypotheses (as above with geological information), a discovery approach was first applied to generate an estimate of the optimal genetic partitioning, and then a validation approach was used to assess the support for these groupings. Similar approaches to delimitation have combined Structurama with BPP (Leaché and Fujita (2010)) and Brownie with BPP (Niemiller et al. 2012). We used Structure to partition the individuals based upon recovered genetic clusters (see above), and then further analyzed these partitions with BPP; analyses for each were conducted as summarized above.

Formal Species Descriptions

Specimens were imaged using a Visionary Digital BK plus system (http://www.visionarydigital.com), including a Canon 40D digital camera, Infinity Optics Long Distance Microscope, P-51 camera controller, and FX2 lighting system. Individual images were combined into a composite image using Zerene Stacker V1.04 software, which was then edited using Adobe Photoshop CS3. Seminal receptacles (=spermathecae) were dissected from adult female specimens using fine forceps, immersed for 2–5 min in BioQuip specimen clearing fluid (www.bioquip.com) on a depression slide, then imaged directly in this fluid on slides. Other images were taken with specimens immersed in filtered 70% EtOH, using KY jelly or fine sand to secure specimens. Specimen measurements were taken from digital images using a calibrated ruler tool in Photoshop CS3. All appendage measurements were recorded from the left appendage, unless otherwise indicated. Measurements mirrored those of Coyle (1974), and only brief character descriptions are provided here. The reader is referred to Coyle (1974, figs. 2–7) for more thorough character definitions: CL, CW—carapace length and width, PCL—length of pars cephalica; IFL, ITL, IML, ITarL—lengths of leg I segments (except for patella), viewed retrolaterally; IVFL, IVTL, IVML, IVTarL—lengths of leg IV segments (except for patella), viewed retrolaterally; PFL, PPL, PTL—lengths of male pedipalp segments, viewed retrolaterally; PTT—maximum depth of male pedipalp; PTX—length of male pedipalp at maximum depth; PED—distance of base of embolus to tip of conductor; PCA—distance from PED line to edge of outer conductor sclerite; SL, SW—sternum length and width; PSS—minimum distance between posterior sigilla; PSL—maximum diameter of right posterior sigilla. All measurements are reported in millimeters (mm).

Figure 2

Species delimitations estimated with validation approaches. Results show final analysis with combined data. a) Results from spedeSTEM, represented on species trees estimated with STEM. Scale bar is in units of 2Ne generations. b) Results from BPP, represented on species tree estimated with *BEAST. Numerical values above nodes represent pp values; asterisks below nodes denote presence of split in BPP analyses. Scale bar represents substitutions per site.

Figure 2

Species delimitations estimated with validation approaches. Results show final analysis with combined data. a) Results from spedeSTEM, represented on species trees estimated with STEM. Scale bar is in units of 2Ne generations. b) Results from BPP, represented on species tree estimated with *BEAST. Numerical values above nodes represent pp values; asterisks below nodes denote presence of split in BPP analyses. Scale bar represents substitutions per site.

RESULTS

Molecular Data

GenBank accession numbers for newly generated sequences are provided in Supplemental Table 1, augmenting those published in Satler et al. (2011). Data matrices and trees used for species delimitation analysis have been uploaded to TreeBase (http://purl.org/phylo/treebase/phylows/study/ TB2:S14119). General information for all gene fragments can be found in Table 2 (e.g., alignment length, number of parsimony informative sites, etc.). Results from TOPALi indicate possible recombination within the 28S rRNA gene, but this is likely an artifact of high among-site rate variation (McGuire et al. 1997). The full COI data set includes 103 individuals from 82 localities; reduced data sets include all 6 gene regions representing 27 localities (5 for a small number of locations). Extensive geographic structure characterizes the mitochondrial DNA data, and individuals from essentially all distinct geographic locations possess haplotypes that are not shared among sampling localities (only two haplotypes are shared across localities; sites 22–24 share a haplotype and sites 29–30 share another haplotype). This extreme sampling-level divergence is reflected in COI genetic distances, where the average pairwise genetic distance between locales is 9.46% (GTR model), with a maximum of 15.92% (between sites 3 and 21). Results from the AMOVA confirmed extensive levels of population structure, with data being significantly partitioned among regions and provinces (Table 3). Results from the Mantel test are also significant for isolation by distance on the combined data set (P = 0.008). Maximum likelihood gene trees derived from autosomal loci are less phylogenetically structured than the COI gene tree, but still include many strongly supported nodes (Supplemental Fig. 1).

TABLE 2

Gene information for multilocus analyses

Locus bp Sequences PI Model 
COI mtDNA 969 27 25 289 213 GTR + gamma (pos 1) 
GTR + gamma (pos 2) 
GTR + gamma (pos 3) 
28S rRNA 1052 54 21 278 265 GTR + gamma 
18S rRNA 754 54 55 45 GTR + I 
EF-1γ nDNA 642 54 21 24 22 GTR + I 
Mid1 nDNA 739 54 13 16 12 K80 + I 
Anonymous nDNA 1023 42 24 182 132 GTR + gamma 
Locus bp Sequences PI Model 
COI mtDNA 969 27 25 289 213 GTR + gamma (pos 1) 
GTR + gamma (pos 2) 
GTR + gamma (pos 3) 
28S rRNA 1052 54 21 278 265 GTR + gamma 
18S rRNA 754 54 55 45 GTR + I 
EF-1γ nDNA 642 54 21 24 22 GTR + I 
Mid1 nDNA 739 54 13 16 12 K80 + I 
Anonymous nDNA 1023 42 24 182 132 GTR + gamma 

Note: Locus information includes name of locus, length (in base pairs), number of sequences, number of unique alleles, number of segregating sites, number of parsimony informative sites, and model of DNA sequence evolution.

TABLE 3

AMOVA results

Source df SS MS 
Among region 34.0648 17.0324 
Among province 10 105.0040 10.5004 
Error 14 76.8571 5.4898 
Total 26 215.9259  
Parameter Value P  
ΦRT 0.0915 0.0001  
ΦPR 0.3157 0.0001  
ΦPT 0.3783 0.0001  
Source df SS MS 
Among region 34.0648 17.0324 
Among province 10 105.0040 10.5004 
Error 14 76.8571 5.4898 
Total 26 215.9259  
Parameter Value P  
ΦRT 0.0915 0.0001  
ΦPR 0.3157 0.0001  
ΦPT 0.3783 0.0001  

Notes: AMOVA results based on all six 6 loci (combined). Samples are partitioned into three regions and thirteen provinces (see Table 1 for details). For each level, degrees of freedom (df), sum of squares (SS), and mean of squares (MS) are shown. Null distributions were generated with 10,000 permutations to test for significance.

Tests of Species Limits

Species limits were tested using multiple approaches, including validation, discovery and a chimeric approach. We derived a priori partitions from geological data, and validated the species limits implied by these partitions using spedeSTEM and BPP. In the south CR, both provinces were collapsed into a single lineage with strong support (wi = 0.73 for spedeSTEM, see Table 4; pp = 0.91 for BPP; Supplemental Fig. 2a,b). In the TR, the model with the highest probability (wi = 0.08) in spedeSTEM recovers Frazier (f) as a distinct lineage, and collapses the remaining provinces (e, g, h, j, l, m) into a single lineage (Supplemental Fig. 2c). BPP recovers three lineages (pp = 0.07) in the TR—Frazier (f), San Emigdio (e), and a lineage comprising the remaining provinces (h, j, l, m; Supplemental Fig. 2d). In the southern SN, two lineages were recovered—one consisting of Sierra (a) and a second lineage consisting of the remaining provinces (b–d) (wi = 0.53 for spedeSTEM; pp = 0.69 for BPP; Supplemental Fig. 2e,f). After recovering the within region lineages, spedeSTEM and BPP analyses were extended to the remaining lineages to further evaluate species diversity. One model in spedeSTEM gathered nearly all of the support (wi = 0.93), and included four lineages (see Fig. 2a): Sierra (a), Frazier (f), and the recovered Transverse Range lineage (e, g, h, j, l, m) were retained, while the lineage consisting of provinces b-d (in the Tehachapi and Piute Mountains) was combined with the south CR (o–p). All six lineages recovered with BPP within regions were retained (pp = 0.99; see Fig. 2b) in the combined BPP analysis (Sierra (a), Frazier (f), San Emigdio (e), (b–d), (h, j, l, m), and the south CR (o–p)). Thus, validation approaches with spedeSTEM and BPP suggest 4 and 6 lineages, respectively.

TABLE 4

spedeSTEM validation results

Region Model AICc Δi Model-likelihood wi 
CR o_p 1240.1033 0.73073447 
o, p 1242.1000 1.9967041 0.3684862 0.26926556 
TR f, e_g_h_j_l_m 1298.1117 0.08322176 
f, e_j_m, g_h_l 1298.6167 0.5050049 0.77685434 0.06465119 
SN a, b_c_d 1015.3145 0.5324434 
a, b_c, d 1017.3106 1.9960938 0.36859867 0.19625793 
Total a, f, b_c_d_o_p, e_g_h_j_l_m 2482.9302 0.9329383 
a, f, b_c_d, e_g_h_j_l_m_o_p 2488.8738 5.9436035 0.05121096 0.04777667 
Region Model AICc Δi Model-likelihood wi 
CR o_p 1240.1033 0.73073447 
o, p 1242.1000 1.9967041 0.3684862 0.26926556 
TR f, e_g_h_j_l_m 1298.1117 0.08322176 
f, e_j_m, g_h_l 1298.6167 0.5050049 0.77685434 0.06465119 
SN a, b_c_d 1015.3145 0.5324434 
a, b_c, d 1017.3106 1.9960938 0.36859867 0.19625793 
Total a, f, b_c_d_o_p, e_g_h_j_l_m 2482.9302 0.9329383 
a, f, b_c_d, e_g_h_j_l_m_o_p 2488.8738 5.9436035 0.05121096 0.04777667 

Note: Results from spedeSTEM validation approach show model, number of free parameters (k), AICc score, AICc differences (Δi), likelihood of the model, and probability of the model (wi).

The validation results could be misleading if the history of A. thompsoni does not reflect landscape features, or if our interpretation of the geological literature or the boundaries delineated in the geological literature are inaccurate. Consequently, we utilized several species discovery approaches that do not require a priori groupings. Results from Structure suggest three genetic partitions, with K = 3 containing the largest ΔK (155.539) as estimated with the Evanno method. The K = 3 partitioning scheme combines Sierra (a) and Frazier (f) into a single group, samples from the Tehachapi and Piute mountains (b–d) into a single group, and the remaining samples into a third group (Fig. 3). Brownie analyses resulted in 10,536 trees due to some runs containing multiple trees with equally high scores. Of the trees recovered, 53 unique lineage groupings were found, with four partitions recovered in over 70% of trees (Fig. 3; Supplemental Table 2); all sampling locales were assigned to 1 of the 4 groups. Notably, a group consisting of the Sierra locales (a) occurred at a frequency of 0.99, and a group consisting of the Frazier locales (f) with a frequency of 0.87. The remaining samples are separated into two groups, with three Transverse Range samples (j, l, m) grouped together (frequency = 0.88), and the remaining samples grouped with frequency of 0.70. Overall, the Structure and Brownie results suggest 3 and 4 independent lineages, respectively.

Figure 3

Species delimitations estimated with discovery approaches. Brownie results are represented on a species tree estimated with *BEAST. Values above terminal branches indicate frequency that grouping was recovered in Brownie; nodal support corresponds to pp values estimated with *BEAST. Values on terminal branches correspond to percentage of Brownie runs where grouping was recovered (out of 10,536 final trees; see Supplemental Table 2). Scale bar is in units of substitutions per site. Structure plot shows genetic clustering of individuals, with K = 3 (Sierra and Frazier are highlighted on terminal branches).

Figure 3

Species delimitations estimated with discovery approaches. Brownie results are represented on a species tree estimated with *BEAST. Values above terminal branches indicate frequency that grouping was recovered in Brownie; nodal support corresponds to pp values estimated with *BEAST. Values on terminal branches correspond to percentage of Brownie runs where grouping was recovered (out of 10,536 final trees; see Supplemental Table 2). Scale bar is in units of substitutions per site. Structure plot shows genetic clustering of individuals, with K = 3 (Sierra and Frazier are highlighted on terminal branches).

Using our novel STEM discovery approach, collapsing of nodes tended to increase the probability of the model (although some iterations stayed the same or dropped slightly; see Fig. 4, Supplemental Table 3), until all samples were collapsed into three lineages. This three-lineage model was overwhelmingly favored with the highest support (wi = 0.99), recovers Sierra (a) and Frazier (f) as distinct, and combines the remaining samples into a single lineage. Although three lineages are inferred using this approach, the composition differs from the three partitions found using Structure (see above). Results from the Bayesian implementation of this approach were more difficult to interpret. We found that the maximum number of OTUs that could be analyzed using BPP was 19, so we subdivided the guide tree estimated in *BEAST into three major clades by dividing from the deepest node toward the tips until all samples within a clade could be analyzed. We then adopted a hierarchical approach and analyzed each selected clade independently before conducting a combined analysis. The three clades tested included the Sierra locales (a), the Frazier locales (f), and the remaining locales (see Supplemental Fig. 3). Sierra (a) and Frazier (f) were each collapsed into a single lineage, with the remaining samples partitioned into five lineages as follows—Tehachapi and Piute Mountains (b-d), San Emigdio (e), south CR (o–p), (g–h), and (j, l, m). The combined analysis further collapses the south CR (o–p), (g–h), and (j, l, m) into a single lineage, and retains all remaining lineages, suggesting a 5 lineage complex (Supplemental Fig. 3). It is important to note that lineages that were retained in the initial hierarchical analysis (o–p; g–h; j, l, m) were combined in the final analysis.

Figure 4

Likelihood plots from discovery STEM approach. X-axis is number of lineages, y-axis is log likelihood values. Each iteration is represented by the calculated phylogeny. An asterisk denotes the iteration with the most likely species delimitatation, representing 3 lineages. Out-group data have been removed for visual purposes.

Figure 4

Likelihood plots from discovery STEM approach. X-axis is number of lineages, y-axis is log likelihood values. Each iteration is represented by the calculated phylogeny. An asterisk denotes the iteration with the most likely species delimitatation, representing 3 lineages. Out-group data have been removed for visual purposes.

Results using single locus data (COI) with bGMYC detected many more putative lineages than other approaches. Analysis of the full COI data set suggests an average of 69 coalescent units (mean = 68.98), with a 95% confidence interval (CI) between 61 and 79 coalescent units (Fig. 5). For the recovered units, 48 are found in greater than 95% of the posterior distribution, suggesting strong support for these groupings. Analysis of the geographically reduced data set recovers 12 coalescent units (mean = 12.0108), with a 95% CI between 2 and 23 coalescent units (Supplemental Fig. 4). None of the recovered groupings are found in greater than 95% of the posterior distribution (highest is 0.8464), suggesting that the reduced data set lacks statistical power, perhaps because of extreme COI divergence and/or an overdispersed geographic sample.

Figure 5

Results from bGMYC analysis on complete COI data set. a) COI gene tree. Nodal support corresponds to pp values from BEAST analysis, with an asterisk signifying >0.95 value. Terminal numbers correspond to sample site. Colors next to numbers correspond with pp for recovered coalescent unit. b) Histogram represents recovered species distribution in bGMYC, reflecting range of uncertainty in coalescent unit estimates.

Figure 5

Results from bGMYC analysis on complete COI data set. a) COI gene tree. Nodal support corresponds to pp values from BEAST analysis, with an asterisk signifying >0.95 value. Terminal numbers correspond to sample site. Colors next to numbers correspond with pp for recovered coalescent unit. b) Histogram represents recovered species distribution in bGMYC, reflecting range of uncertainty in coalescent unit estimates.

For the chimeric approach, the three genetic clusters from Structure were further analyzed for lineage distinctness with BPP. Analyses across all runs recovered all three lineages as distinct, with pp values of 1.0 for all nodes. Thus, the chimeric approach suggests a complex consisting of 3 lineages.

Survey of Morphology

Morphological character states that are most useful in distinguishing different Aliatypus species elsewhere in the genus (e.g., configuration of female sperm storage structures, shape of male palpal bulb [used for sperm transfer], spination of male first leg [used during copulation]) are conserved in the A. thompsoni species group (see Appendix, Supplemental Figs. 5–10). Although Coyle (1974) noted geographic morphological variation in A. thompsoni (e.g., variation in female spermathecal morphology, number of trichobothria on metatarsus IV), his sample in fact did not include spiders within the range of either the Sierra (a) or Frazier (f) lineages. However, although subtle variation is detected (e.g., see Supplemental Figs. 7 and 8), no discrete groupings are recovered, with our qualitative analysis of morphology indicating conservative evolution in the A. thompsoni species group, consistent with the hypothesis of cryptic speciation.

DISCUSSION

Variation in Species Delimitation Across Methods

Accuracy in species delimitation is a function of both the data collected and the appropriateness of choices made by the researcher during the course of analysis. Because it is not clear which approach is most appropriate in the A. thompsoni complex, we utilized multiple approaches for species delimitation, including validation, discovery, and a chimeric approach. Our results consistently indicate that there are multiple cryptic lineages within A. thompsoni, but differ in how such lineages are delimited across approaches. Some of the approaches used are likely inappropriate given our data; for example, we view single locus approaches as fundamentally unsuitable in this system (see below).

Validation approaches appear promising for species delimitation but require samples to be partitioned prior to analysis, a constraint that can be problematic. In order to satisfy this requirement, our initial partitions were developed by dividing samples based on geologic criteria. These validation approaches recovered between 4 and 6 lineages, and for the most part, our results seem biologically plausible. One grouping that is not recovered in any other analysis is the joining of spiders from the Tehachapi and Piute Mountains (b–d) with spiders from the south CR (o–p) in the spedeSTEM analysis. However, the validation analyses make several assumptions, including that species limits in A. thompsoni reflect landscape features and that the three regions can be tested in a hierarchical manner (due to too many models with the full data set; see Materials and Methods). We also assume that the geological data and our interpretation of these data are accurate, but this assumption is difficult for spider biologists and geneticists to verify. In addition, the geologic units are adjacent and generally lack obvious extrinsic barriers, such that samples from the edges of defined units might be potentially misplaced. For example, data from mtDNA (as well as 18S rRNA; see Satler 2011) supports locality 28 to be grouped with Frazier (f), yet this location was grouped with San Emigdio (e) based on geologic criteria. In short, validation approaches are computationally more efficient but potentially misleading, and should be conducted in concert with discovery approaches. It is also important to note that defining a priori partitions based on other types of data, such as morphology or taxonomy, can suffer from the same issues, where a priori partitioning of the samples could be incorrect (e.g., via convergence in morphology, inaccurate taxonomic subdivisions, etc.).

Given the lack of preexisting divisions in A. thompsoni, discovery approaches seem essential. The four multilocus discovery approaches used in this study recovered between 3 and 5 lineages; two of the groupings are recovered in three of the approaches, with incongruence seen regarding the partitioning of remaining samples. We attribute these differences to characteristics of the individual methods. For example, both Brownie and STEM use single gene trees as input, and thus do not integrate over the uncertainty in gene tree space as does the Bayesian method BPP. However, the latter approach does not account for rate heterogeneity within or across loci, and does not allow for more parameter-rich DNA substitution models (the user is restricted to a Jukes–Cantor model), thus representing a less accurate representation of the gene tree distribution. In the case of the A. thompsoni data, gene tree estimates appear to be reasonably well supported (based upon nodal support values; see Supplemental Fig. 1). A more important factor may involve the exploration of tree space with varying species limits. Brownie uses a heuristic search because tree space incorporating uncertainty in species limits (i.e., varying OTU compositions) is exceedingly large. Nevertheless, a majority of the sampled Brownie trees recover four groupings (Fig. 3, Supplemental Table 2), and these groupings are generally consistent with the geographical divisions of the samples (Fig. 1).

We observed that BPP consistently recovers the highest number of partitions as compared with other methods in both the discovery and validation approaches. One explanation for this result might be that BPP does not integrate over the species tree parameter space. Whereas a species tree is required as an input guide tree, previous work has shown that an incorrectly specified species tree can produce false positives and thus result in oversplitting (see Leaché and Fujita 2010). Although simulations demonstrate the efficacy of this approach (Yang and Rannala 2010; Zhang et al. 2011; Camargo et al. 2012) given an accurate guide tree, the simulated scenarios are much simpler than the complex scenario we face in this system. More generally, species trees can never be known with certainty in empirical systems. For our analysis, we assume that the accuracy of our estimates of individual gene trees are easier to verify than our species tree estimates, because the latter necessarily include species delimitations.

We also used a genetic clustering algorithm implemented in Structure to partition the data; simulation studies show that these algorithms may work well for species delimitation (Rittmeyer and Austin 2012). The Structure results recover three groups as do several other approaches, but Structure uniquely groups Sierra (a) and Frazier (f) into a single lineage. These two groups are also clustered in the second ranked model (K = 2; ΔK = 111.037), but are recovered as distinct in the third ranked model (K = 4; ΔK = 32.869). Our determination of three partitions was based on the popular Evanno method, but visual assessment of the Structure plot (Fig. 3) suggests differing genetic compositions between the two groups, highlighting the potential difficulty in summarizing the optimal K value from a Structure analysis. Although Structure has been widely used for investigations of population structure, the isolation by distance signal in our data suggest that it may not be well suited to our system, and may be a contributing factor (among potentially others) driving this possibly spurious grouping (Pritchard et al. 2000; Pritchard et al. 2010). Given that the grouping of Sierra (a) and Frazier (f) together was not supported by any other analysis, this can certainly be problematic if these partitions are used for downstream analyses. For example, this incompatibility may invalidate the chimeric approach (Structure/BPP) that has been utilized in other investigations. Each of the partitions recovered using Structure were validated with BPP (pp = 1.0 for all nodes across all runs). It seems possible that such chimeric approaches could validate incorrect partitions (given that validation approaches cannot split initial groupings) and mislead the researcher regarding species limits, either directly or via errors in species tree estimation, when a method that requires a guide tree is used (e.g., BPP).

Our lack of satisfaction with existing approaches led us to develop a novel approach to species discovery using species trees. This STEM approach is implemented in the same manner as spedeSTEM, but operates by working backwards from the individual sample (in our case, sample sites) rather than by validating a priori partitions. Species divergence is estimated during the process of calculating the P (ST | GT) for a given partition, and the only assumption is that members of a particular lineage (i.e., either a single or group of samples) will be more closely related genetically than they will to members of other lineages. Starting with a species tree that treats each individual sample as an OTU, our novel discovery method works backwards by collapsing the tips of the species tree that are most similar. At each level of clustering STEM is used to calculate the ML tree given the data. Both Sierra (a) and Frazier (f) are recovered as distinct lineages, a result which is found across essentially all other methods. All other sampling localities are collapsed into a single lineage, as opposed to dividing these into various additional groups as seen in the other discovery, validation and chimeric approaches (see Fig. 6).

Figure 6

Summarized results from all species delimitation analyses, represented on *BEAST tree resulting from analysis of all sampling locales. Colors correspond to recovered groupings, with each partition represented by a unique color. Colored bars above phylogeny represent hypothesized species groupings based on multiple analyses. Insert image of an adult male Aliatypus starretti, sp. nov. (Kern Co., Poso Flat Road).

Figure 6

Summarized results from all species delimitation analyses, represented on *BEAST tree resulting from analysis of all sampling locales. Colors correspond to recovered groupings, with each partition represented by a unique color. Colored bars above phylogeny represent hypothesized species groupings based on multiple analyses. Insert image of an adult male Aliatypus starretti, sp. nov. (Kern Co., Poso Flat Road).

Challenges of Species Delimitation in Genetically Fragmented Systems

Mygalomorph spiders, and genetically fragmented systems in general, are challenging groups for species delimitation. Previous studies investigating these spiders with genetic data find extensive population structure, with the presence of numerous “micro-allopatric” populations within species (e.g., Bond et al. 2001; Starrett and Hedin 2007; Bond and Stockman 2008; Hedin et al. 2013). The potential to oversplit diversity is high in such systems. For example, we used an implementation of the GMYC model that takes into account gene tree uncertainty (i.e., bGMYC) in a single locus data set. Although the GMYC model has previously been used to delimit arthropod species (Papadopoulou et al. 2008; Hamilton et al. 2011; Pons et al. 2011; Keith and Hedin 2012), it relies upon a single locus and is subject to the potential errors associated with such data (e.g., selection, incomplete lineage sorting, gene tree estimation error). The GMYC model appears to be most accurate if a “barcode gap” is present (Pons et al. 2006), because coalescent events within the gene tree are expected to reside in units with little to no genetic structure (i.e., recovered coalescent units are assumed to be panmictic). However, the pronounced population structure within our system (and most mygalomorph spiders) is likely to bias the GMYC model by shifting the transition of coalescent to divergence processes toward the present, therefore potentially recognizing a very high number of species (Lohse 2009; Hamilton et al. 2011). This highlights a potential major flaw with the GMYC model, which is the confounding of population-level structure with species boundaries in taxa containing limited dispersal abilities. We suspect that our results reflect this bias (see Figs.5 and 6), and doubt that dozens of cryptic species are present within A. thompsoni. Rather, what is being detected is the signature of female life history traits, where individual females display minimal dispersal abilities and exist in isolated “micro-allopatric” populations. Because only males are expected to display appreciable dispersal behavior (and thereby move alleles across the landscape), it is clear that the inclusion of multiple nuclear markers is necessary for understanding species limits within this group, and the many other organismal groups with similar natural histories.

Another important consideration in genetically fragmented groups is geographic sampling, and the problem of sampling gaps. In the A. thompsoni complex, the region between the central TR and Tehachapi Mountains is the large, privately owned Tejon Ranch, which we were unable to access for this study (see Fig. 1). Are Tehachapi and Piute Mountains spiders genetically distinct from spiders in the central TR (see Fig. 6) because of true genealogical breaks, or is this distinctiveness an artifact of a sampling gap? Because incomplete sampling can lead to perceived phylogenetic breaks (Irwin 2002), we cannot be certain that this is not the case here. If a pattern of isolation by distance is present within this group, we might expect that the inclusion of samples from the Tejon sampling gap could potentially collapse the lineages detected within the remaining samples. Alternatively, genetic breaks detected by these species delimitation methods may be evolutionary breaks, and the inclusion of more samples would only reinforce these findings and provide greater detail on where the phylogenetic breaks occur. Given the results from the Mantel Test, a signature of isolation by distance is present in this system, suggesting that this sampling gap may be a driving force behind analyses splitting the samples (see Jackman and Wake 1994). This uncertainty highlights the need for dense geographic sampling when trying to understand species limits in genetically fragmented systems (e.g., Weisrock et al. 2010; Jockusch et al. 2012; Niemiller et al. 2012; Hedin et al. 2013).

The Difficult but Importance Task of Formally Describing Cryptic Species

We hypothesize that A. thompsoni is a species complex, and recommend the elevation of Sierra (a) and Frazier (f) lineages to species status under criteria of the general lineage concept (de Queiroz 2007). Below (see Appendix) we provide formal species descriptions for the Sierra (a) and Frazier (f) lineages. In light of the incongruence among methods regarding species limits for the remaining samples in the A. thompsoni complex, we advocate a conservative approach and consider these samples as belonging to a single species. Because these remaining samples include the type locality, they retain the name A. thompsoni. In addition, although the multigenic approach adopted in our research does not require reciprocal monophyly for any single gene region for any putative group, both the Sierra (a) and Frazier (f) lineages in fact meet these criteria (see Supplemental Fig. 1). The Frazier (f) lineage is recovered as a clade on all sampled gene trees, typically with strong support. The Sierra (a) lineage is recovered as a clade on 4 of 5 sampled gene trees (only a single population was sampled for Anonymous). We calculated the probability of such a high degree of genealogical congruence occurring by chance using the multigenic Yule process model of Rosenberg (2007, equation 6). For both lineages, very low probabilities (Frazier = 1.23 × 10−28, Sierra = 2.80 × 10−15; see Supplemental Material) provide strong additional support for the evolutionary and genetic distinctiveness of these groups.

Our survey of morphology suggests that both male and female spiders in the Sierra (a) and Frazier (f) lineages retain conserved character states shared across the entire A. thompsoni species group. However, we are not arguing that morphological differences do not exist in this complex—different methods of study or quantitative analyses may indeed uncover such differences in future studies. We also emphasize that although male secondary reproductive structures are commonly used in mygalomorph spider systematics, the few adult male specimens available for study (one male from each new species) precludes a thorough statistical analysis at this time. This rarity of specimens partially reflects the difficulty in collecting adult male Aliatypus, which after molting to adulthood, seek mates for brief periods of time prior to death (Coyle and Icenogle 1994).

Morphologically cryptic species pose problems for traditional animal taxonomy, which is founded on the discovery of morphologically diagnosable groups for both species delimitation and description. The infusion of molecular data in systematics has seen a proliferation of the “discovery” of morphologically cryptic species, but most (a majority?) of these putative taxa await formal description. There are many possible reasons for this mismatch between species discovery and description, including less taxonomic training in systematics, conservative researchers waiting for the accumulation of additional data, and the difficulty of describing species with nontraditional characters under the traditional paradigm (Fujita and Leaché 2010). This last problem is particularly acute for morphologically cryptic species, where diagnostic morphological features are few to absent, thus forcing a reliance on diagnostic DNA characters. Such diagnostic DNA characters are numerous in the A. thompsoni complex (see Appendix), but will not always be available, and their recognition is perhaps inconsistent with a multispecies coalescent perspective (Fujita and Leaché 2010). Despite these difficulties, we view the disconnect between “discovery” and formal description as fundamentally problematic. The process of species delimitation will continue to become a more data-rich process, and these data need to be applied to both species delimitation and description. Species are fundamental units of evolutionary and conservation biology, and formal species description must be an important tenet of systematic biology.

CONCLUSIONS

This study provides a framework for the discovery of cryptic lineage diversity that prevails in many organismal groups. Here we show the importance of taking a combined approach by exploring the performance of both discovery and validation methods, where groupings recovered across various approaches provide a strong signal for the presence of distinct species. In addition to currently available methods, our novel species tree-based discovery approach with STEM recovered a consensus three species model with overwhelming support, suggesting that this method can be effective in recovering species limits in genetically fragmented systems. Our results also suggest that single locus data sets may be unsuitable in such systems, where available methods (e.g., GMYC model) are likely to recover much higher levels of species diversity than is actually present. This study further adds to our specific knowledge of California mygalomorph diversity, and generally increases our taxonomic knowledge in this biodiversity hotspot. New tools for data collection and analysis are allowing systematists to discover biodiversity like never before; we argue that this discovery process must become more directly linked to species descriptions and formal taxonomy.

SUPPLEMENTARY MATERIAL

Supplementary figures, tables and material can be found at http://datadryad.org in the Dryad data repository (doi:10.5061/dryad.68s40).

FUNDING

This work was supported by the American Arachnological Society (Vincent Roth and Arachnological Research Funds to J.D.S.) and San Diego State University (Theodore Cohn Evolutionary Biology Graduate Scholarship, Harry E. Hamber Memorial Scholarship to J.D.S.).

ACKNOWLEDGMENTS

The authors would like to thank Jim Starrett, Jason Bond, David Carlson, Dean Leavitt, Casey Richart, Shahan Derkarabetian, Axel Schönhofer, Jared Grummer, Jenny Nash, Pieter Van Niekerk, Bob Keith, Lars Hedin, Ian Ballard, Steve Lew, Pierre Paquin, and Steven Thomas for help with the collection of spiders. We would like to sincerely thank Jason Bond for providing Antrodiaetus transcriptomic data used to develop molecular markers for this study, Helen Fan for help with calculations of taxonomic distinctiveness, and Noah Reid for assistance with bGMYC. We thank Sarah Hird, John McVay, Tara Pelletier, Jessica Perez, and Noah Reid for helpful discussion and comments on earlier versions of this article. We would also like to thank Robb Brumfield, Matthew Niemiller, Jason Bond, and one anonymous reviewer for valuable comments and suggestions that helped to improve this article.

APPENDIX—SPECIES DESCRIPTIONS

Systematics

Family Antrodiaetidae Gertsch, in Comstock, 1940: 236 [urn:lsid:amnh.org:spiderfam:0002]

Genus Aliatypus Smith, 1908 [urn:lsid:amnh.org: spidergen:00006]

Aliatypus thompsoni Species Group

Group of three closely related taxa, sharing the following morphological characteristics which (in combination) distinguish this group from all other described Aliatypus species: Males possess large, closely spaced posterior sternal sigilla, a relatively long leg I metatarsus, short ventral ensiform macrosetae and sparse dorsal macrosetae on leg I tibia and metatarsus, and an inconspicuous thoracic groove. Females possess an inconspicuous thoracic groove, closely spaced posterior sternal sigilla, few- to many-looped seminal receptacles with medium-sized terminal bulbs, and a large number of ensiform macrosetae on the retrolateral surface of the palpal tarsus.

Aliatypus thompsoni Coyle 1974

Aliatypus thompsoni Coyle 1974: figs. 27–33, 40–44, 50–53, 60, 71–73, 86, 93, 111–113, 155–171; map 4.

For diagnosis and description see Coyle (1974). None of the specimens examined by Coyle occur within the geographic distribution of the two new species described below, and as such, we attribute all such specimens to A. thompsoni. A single possible exception is the female examined by Coyle from the “Piute Mountains, S of Kernville.” New records for A. thompsoni include all populations examined for this study (see Supplemental Table 1), except those from collecting sites 1–6 and 25–30 (see Fig. 1).

Aliatypus roxxiae Satler and Hedin, new species.

Figures 1 and 6, Supplemental Figures 5–8

Type material—male holotype from California, Kern County, Mt Pinos Road, south of Mil Portrero Hwy., 34.8276, −119.0872, collected 23 August, 2009, J. Starrett & E. Floyd (voucher specimen MY4040). GenBank number: KC787766. Female paratype from California, Kern County, Cuddy Valley Road, 34.8168, −119.0765, collected 21 May 2009, J. Satler (MY1742). GenBank numbers: KC787767, KC787890, KC787731, KC787704, KC787944, KC787916. Specimens deposited at the California Academy of Sciences (CAS), San Francisco, California.

Other Records—California, Kern County: near Antelope Canyon, northwest corner of Antelope Valley, 34.8504, −118.6905, collected 26 July 2010, M. Hedin, J. Satler, D. Carlson & L. Hedin (two females (F), one immature (I)); east of Frazier Park, near I-5, 34.8138, −118.8897, collected 17 January, 2003, M. Hedin & J. Starrett (F, 2I); Lakeview Drive, Lake of the Woods, 34.8133, −119.0027, collected 20 May 2009, J. Satler (2I); Mt Pinos Road, 34.8276, −119.0872, collected 23 August, 2009, J. Starrett & E. Floyd (2F, I); Cuddy Valley Road, 34.8168, −119.0765, collected 21 May, 2009, J. Satler (F, 2I). Inclusion of specimens from site 28 (Kern County: west of Frazier Park, 34.8254, −119.0115, collected 31 March 2010, J. Satler (2I)), based on COI gene tree data, plus unpublished 18S data (Satler 2011).

Diagnosis—Corresponds to Frazier lineage (e.g., Fig. 6). Difficult to distinguish morphologically from other members of the A. thompsoni species group, sharing species group features as defined above. Some, but not all populations of A. thompsoni include more than one distal trichobothria on metatarsus IV (both sexes), and fewer seminal receptacle loops (females), but these character differences are not strictly diagnostic because of geographic variation within A. thompsoni (see Coyle 1974). Morphologically very similar to A. starretti, sp. nov.

Best distinguished from other members of the A. thompsoni species group based on allopatric geographic distribution (from collecting sites 25–30, Fig. 1), plus diagnostic nucleotide changes for multiple gene regions. The following subset of site changes are diagnostic for the easily sequenced 18S gene region (positions correspond to 18S alignments available on TreeBase): site 204 T (vs. C in remaining members of complex), site 633 G (vs. T), site 677 G (vs. T).

Holotype male (MY4040)—Carapace 5.0 long, 3.8 wide, pars cephalica 2.9 long. Thoracic groove a very slight double longitudinal depression. Postocular setae forming short narrow longitudinal row. Leg I IFL 5.4, ITL 3.2, IML 3.6, ITarL 1.5. Tibia and metatarsus with conspicuous ventral ensiform macrosetae; a few scattered attenuate macrosetae on dorsal tibia, but most setae appressed. Leg IV IVFL 4.7, IVTL 3.0, IVML 4.5, IVTarL 2.2. One trichobothrium dorsally on distal end of metatarsus IV. Pedipalp PFL 5.3, PPL 3.1, PTL 3.5, PTX 2.7, PTT 0.9, PED 1.1, PCA 0.5. Palpal tibia strongly swollen ventrally near distal end, hirsute on ventral surface of swelling. Conductor tip broad, with distal nipple-like tip. Sternum SL 2.8, SW 2.3, posterior sternal sigilla close (PSS 0.2), large (PSL 0.6), suboval. Abdomen Tergites I and III much smaller than tergite II.

Paratype female (MY1742)—Carapace 6.2 long, 4.5 wide, pars cephalica 4.0 long. Thoracic groove essentially absent. Postocular setae forming narrow longitudinal row. Chelicerae with row of 4 retrolateral macroteeth, 7 prolateral macroteeth, 6 intermediate microteeth. Leg I IFL 4.1, ITL 2.3, IML 1.8, ITarL 1.0. 14 ventral ensiform macrosetae on metatarsus. Leg IV IVFL 3.8, IVTL 2.1, IVML 3.0, IVTarL 1.2. One trichobothrium dorsally on distal end of metatarsus IV. Pedipalp tarsus with 6 ensiform macrosetae on both prolateral and retrolateral surfaces. Sternum SL 4.0, SW 3.3, posterior sternal sigilla close (PSS = 0.3), large (PSL = 0.7), faint, oblong. Seminal receptacles stalks long, many-looped (5–7 bends), with medium-sized terminal bulbs (relative to stalk diameter).

Female variation—All other adult females available (see records above) possess an inconspicuous to absent thoracic groove, and spermathecal receptacles consistent with the paratype female but with slightly fewer (4–6) bends. With a single exception, all females possess a single trichobothrium distally on metatarsus IV (one specimen has two).

Etymology—Named after Roxxi, a favorite family dog of the Satler family.

Distribution and natural history—Most populations known from oak/pine woodlands on north side of Frazier Mountain and Mt. Pinos, south of the San Andreas Rift Zone. Easternmost known population (lower Antelope Canyon) from northwest corner of Antelope Valley, east of Tejon Pass. All collections are from mesic microhabitats, particularly N-facing, stabilized soil banks from beneath trees. Adult females have been collected from burrows throughout the year; an adult male was extracted from a burrow in late August.

Aliatypus starretti Satler and Hedin, new species.

Figures 1 and 6, Supplemental Figures 7–10

Type material—Male holotype from California, Kern County, Poso Flat Road, south of Glenville, 35.6215, −118.7077, collected 31 October 2009, J. Satler (MY4061). Female paratype from California, Kern County, Hwy 155, east of Glenville, 35.7215, −118.6946, collected 30 March 2003, M. Hedin, J. Starrett, P. Paquin, S. Lew (MY485). GenBank numbers: JN540748, JN540564, JN540537, KC787700, KC787939, KC787911, KC787912. Specimens deposited at the California Academy of Sciences (CAS), San Francisco, California.

Other Records—California, Kern County: Hwy 155, east of Glenville, 35.7215, −118.6946, collected 30 March 2003, M. Hedin, J. Starrett, P. Paquin, S. Lew (2F, I); Poso Flat Road, 35.6215, −118.7077, collected 19 May, 2009, J. Satler (2F, I); off Sierra Way, northeast of Kernville, 35.7605, −118.412, collected 19 May, 2009, J. Satler (I); near Lake Isabella dam, off Keyesville Road, 35.6358, −118.4963, collected 30 March 2003, M. Hedin, J. Starrett, P. Paquin, S. Lew (F, 3I). Inclusion of specimens from following sites based on COI gene tree data: California, Kern County: Cook Peak Lookout Road, south of Mountain Mesa, 35.6205, −118.4219, collected 30 March, 2010, J. Satler et al. (F, I); Erskine Creek Road, southeast of Lake Isabella, 35.5869, −118.4383, collected 30 March 2010, J. Satler et al. (2I).

Diagnosis—Corresponds to Sierra lineage (Fig. 6). Difficult to distinguish morphologically from other members of the A. thompsoni group, sharing species group features as defined above. Some, but not all populations of A. thompsoni include more than one distal trichobothria on metatarsus IV (both sexes), and more seminal receptacle loops (females), but these character differences are not strictly diagnostic because of geographic variation within A. thompsoni (see Coyle 1974). Morphologically very similar to A. roxxiae.

Best distinguished from other members of the A. thompsoni group based on allopatric geographic distribution (from collecting sites 1–6, Fig. 1), plus diagnostic nucleotide changes for multiple gene regions. The following subset of site changes are diagnostic for the easily-sequenced 18S gene region (positions correspond to 18S alignment available on TreeBase): site 11 A (vs. T in remaining members of complex), site 18 G (vs. A), site 133 T (vs. C), site 165 C (vs. A).

Holotype male (MY4061)—Carapace 4.6 long, 3.7 wide, pars cephalica 2.8 long. Thoracic groove a very slight double longitudinal depression. Postocular setae forming short narrow longitudinal row. Leg I IFL 5.5, ITL 3.4, IML 3.7, ITarL 1.7. Tibia and metatarsus with conspicuous ventral ensiform macrosetae; a few scattered attenuate macrosetae on dorsal tibia, but most setae appressed. Leg IV IVFL 4.8, IVTL 3.2, IVML 5.3, IVTarL 2.4. One trichobothrium dorsally on distal end of metatarsus IV. Pedipalp PFL 5.4, PPL 3.5, PTL 3.8, PTX 3.0, PTT 0.8, PED 0.9, PCA 0.5. Palpal tibia strongly swollen ventrally near distal end, hirsute on ventral surface of swelling. Conductor tip sharp, like knife blade. Sternum SL 2.7, SW 2.2, posterior sternal sigilla close (PSS 0.2), large (PSL 0.6), suboval, mottled in appearance. Abdomen Tergites I and III much smaller than tergite II.

Paratype female (MY485)—Carapace 5.3 long, 3.7 wide, pars cephalica 3.4 long. Thoracic groove essentially absent. Postocular setae forming narrow longitudinal row. Chelicerae with row of 4 retrolateral macroteeth, 7 prolateral macroteeth, 8 intermediate microteeth. Leg I IFL 3.1, ITL 1.9, IML 1.5, ITarL 0.8. 13 ensiform macrosetae on metatarsus. Leg IV (right) IVFL 3.3, IVTL 2.1, IVML 2.7, IVTarL 1.0. One trichobothrium dorsally on distal end of metatarsus IV. Pedipalp tarsus with 7 ensiform macrosetae on both prolateral and retrolateral surfaces. Sternum SL 3.1, SW 2.5, posterior sternal sigilla close (PSS = 0.2), large (PSL = 0.8), faint, oblong. Seminal receptacles stalks long, several loops (3–4 bends), with medium-sized terminal bulbs (relative to stalk diameter).

Variation—All other adult females available (see records above) possess an inconspicuous to absent thoracic groove, spermathecal receptacles consistent with the paratype female (3–5 bends), and a single trichobothrium distally on metatarsus IV.

Etymology—Named after Dr. James Starrett, outstanding friend, student of mygalomorph spiders, and collector of many Aliatypus specimens used in this study.

Distribution and natural history—Most populations from upper Kern River valley, in woodland habitats surrounding Lake Isabella. Western populations occur in intermediate elevation oak woodland habitats on west side of Greenhorn Mountains near and south of Glenville (see Fig. 1). Collections from N-facing mesic microhabitats, including stabilized soil banks beneath trees, shaded roadcuts, and from the base of large boulders. Adult females have been collected from March to May (but must be in burrows year-round after molting to adulthood); an adult male was extracted from a burrow in late October.

REFERENCES

Arnedo
M.A.
Ferrandez
M.
Mitochondrial markers reveal deep population subdivision in the European protected spider Macrothele calpeiana (Walckenaer, 1805) (Araneae, Hexathelidae)
Conserv. Genet.
 , 
2007
, vol. 
8
 (pg. 
1147
-
1162
)
Bailey
T.H.
Jahns
R.H.
Geology of the Transverse Range province, southern California
California Div. Mines and Geology Bull.
 , 
1954
, vol. 
170
 (pg. 
83
-
106
)
Barrett
C.F.
Freudenstein
J.V.
An integrative approach to delimiting species in a rare but widespread mycoheterotrophic orchid
Mol. Ecol.
 , 
2011
, vol. 
20
 (pg. 
2771
-
2786
)
Beerli
P.
Felsenstein
J.
Maximum-likelihood estimation of migration rates and effective population numbers in two populations using a coalescent approach
Genetics.
 , 
1999
, vol. 
152
 (pg. 
763
-
773
)
Beerli
P.
Felsenstein
J.
Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach
Proc. Natl. Acad. Sci. USA.
 , 
2001
, vol. 
98
 (pg. 
4563
-
4568
)
Bickford
D.
Lohman
D.J.
Sodhi
N.S.
Ng
P.K.L.
Meier
R.
Winker
K.
Ingram
K.K.
Das
I.
Cryptic species as a window on diversity and conservation
Trends Ecol. Evol.
 , 
2007
, vol. 
22
 (pg. 
148
-
155
)
Bond
J.E.
Phylogenetic treatment and taxonomic revision of the trapdoor spider genus Aptostichus Simon (Araneae: Mygalomorphae: Euctenizidae)
Zookeys.
 , 
2012
, vol. 
252
 (pg. 
1
-
209
)
Bond
J.E.
Beamer
D.A.
Lamb
T.
Hedin
M.
Combining genetic and geospatial analyses to infer population extinction in mygalomorph spiders endemic to the Los Angeles region
Animal Conservation.
 , 
2006
, vol. 
9
 (pg. 
145
-
157
)
Bond
J.E.
Hedin
M.C.
Ramirez
M.G.
Opell
B.D.
Deep molecular divergence in the absence of morphological and ecological change in the Californian coastal dune endemic trapdoor spider Aptostichus simus
Mol. Ecol.
 , 
2001
, vol. 
10
 (pg. 
899
-
910
)
Bond
J.E.
Hendrixson
B.E.
Hamilton
C.A.
Hedin
M.
A reconsideration of the classification of the spider infraorder mygalomorphae (Arachnida: Araneae) based on three nuclear genes and morphology
PLoS ONE.
 , 
2012
, vol. 
7
 
6
pg. 
e38753
 
Bond
J.E.
Stockman
A.K.
An integrative method for delimiting cohesion species: finding the population-species interface in a group of Californian trapdoor spiders with extreme genetic divergence and geographic structuring
Syst. Biol.
 , 
2008
, vol. 
57
 (pg. 
628
-
646
)
Brito
P.H.
Edwards
S.V.
Multilocus phylogeography and phylogenetics using sequence-based markers
Genetica.
 , 
2009
, vol. 
135
 (pg. 
439
-
455
)
Burnham
K.P.
Anderson
D.A.
Model selection and multimodel inference: a practical information-theoretic approach
2002
2nd ed
New York
Springer
pg. 
488
 
Buwalda
J.P.
Geology of the Tehachapi Mountains, California
California Div. Mines and Geology Bull.
 , 
1954
, vol. 
170
 (pg. 
131
-
142
)
Calsbeek
R.
Thompson
J.N.
Richardson
J.E.
Patterns of molecular evolution and diversification in a biodiversity hotspot: the California Floristic Province
Mol. Ecol.
 , 
2003
, vol. 
12
 (pg. 
1021
-
1029
)
Camargo
A.
Morando
M.
Avila
L.J.
Sites
J.W.
Species delimitation with ABC and other coalescent-based methods: a test of accuracy with simulations and an empirical example with lizards of the Liolaemus darwinii complex (Squamata: Liolaemidae)
Evolution.
 , 
2012
, vol. 
66
 (pg. 
2834
-
2849
)
Carstens
B.C.
Dewey
T.A.
Species delimitation using a combined coalescent and information-theoretic approach: an example from North American Myotis bats
Syst. Biol.
 , 
2010
, vol. 
59
 (pg. 
400
-
414
)
Carstens
B.C.
Satler
J.D.
The carnivorous plant described as Sarracenia alata contains two cryptic species
Biol. J. Linn. Soc.
 , 
2013
 
doi: 10.1111/bij.12093
Castresana
J.
Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis
Mol. Biol. Evol.
 , 
2000
, vol. 
17
 (pg. 
540
-
552
)
Chatzimanolis
S.
Caterino
M.S.
Toward a better understanding of the “transverse range break”: lineage diversification in southern California
Evolution.
 , 
2007
, vol. 
61
 (pg. 
2127
-
2141
)
Coyle
F.A.
Systematics of the trapdoor spider genus Aliatypus (Araneae: Antrodiaetidae)
Psyche.
 , 
1974
, vol. 
81
 (pg. 
431
-
500
)
Coyle
F.
Aerial dispersal by mygalomorph spiderlings (Araneae, Mygalomorphae)
J. Arachnology.
 , 
1983
, vol. 
11
 (pg. 
283
-
286
)
Coyle
F.A.
Cladistic analysis of the species of the trapdoor spider genus Aliatypus (Araneae, Antrodiaetidae)
J. Arachnology.
 , 
1994
, vol. 
22
 (pg. 
218
-
224
)
Coyle
F.A.
Icenogle
W.R.
Natural history of the Californian trapdoor spider genus Aliatypus (Araneae, Antrodiaetidae)
J. Arachnology.
 , 
1994
, vol. 
22
 (pg. 
224
-
255
)
Davis
E.B.
Koo
M.S.
Conroy
C.
Patton
J.L.
Moritz
C.
The California hotspots project: identifying regions of rapid diversification of mammals
Mol. Ecol.
 , 
2008
, vol. 
17
 (pg. 
120
-
138
)
de Queiroz
K.
Species concepts and species delimitation
Syst. Biol.
 , 
2007
, vol. 
56
 (pg. 
879
-
886
)
Drummond
A.J.
Ashton
B.
Buxton
S.
Cheung
M.
Cooper
A.
Duran
C.
Field
M.
Heled
J.
Kearse
M.
Markowitz
S.
Moir
R.
Stones-Havas
S.
Sturrock
S.
Thierer
T.
Wilson
A.
Geneious v5.4
2011
 
Available from: http://www.geneious.com (last accessed November 15, 2012)
Drummond
A.J.
Rambaut
A.
BEAST: Bayesian evolutionary analysis by sampling trees
BMC Evol. Biol.
 , 
2007
, vol. 
7
 pg. 
214
 
Dyer
R.J.
GeneticStudio: a suite of programs for spatial analysis of genetic-marker data
Mol. Ecol. Resour.
 , 
2009
, vol. 
9
 (pg. 
110
-
113
)
Earl
D.A.
vonHoldt
B.M.
STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method
Conserv. Genet. Resour.
 , 
2011
, vol. 
4
 (pg. 
359
-
361
)
Edwards
S.V.
Is a new and general theory of molecular systematics emerging?
Evolution.
 , 
2009
, vol. 
63
 (pg. 
1
-
19
)
Ence
D.D.
Carstens
B.C.
SpedeSTEM: a rapid and accurate method for species delimitation
Mol. Ecol. Resour.
 , 
2011
, vol. 
11
 (pg. 
473
-
480
)
Evanno
G.
Regnaut
S.
Goudet
J.
Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study
Mol. Ecol.
 , 
2005
, vol. 
14
 (pg. 
2611
-
2620
)
Excoffier
L.
Smouse
P.E.
Quattro
J.M.
Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data
Genetics.
 , 
1992
, vol. 
131
 (pg. 
479
-
491
)
Fujita
M.K.
Leaché
A.D.
A coalescent perspective on delimiting and naming species: a reply to Bauer et al
Proc. Royal Society B.
 , 
2011
, vol. 
278
 (pg. 
493
-
495
)
Fujita
M.K.
Leaché
A.D.
Burbrink
F.T.
McGuire
J.A.
Moritz
C.
Coalescent-based species delimitation in an integrative taxonomy
Trends Ecol. Evol.
 , 
2012
, vol. 
27
 (pg. 
480
-
488
)
Guindon
S.
Gascuel
O.
A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood
Syst. Biol.
 , 
2003
, vol. 
52
 (pg. 
696
-
704
)
Hall
C.A.J.
Nearshore marine paleoclimate regions, increasing zoogeographic provinciality, molluscan extinctions, and paleoshorelines, California: late Oligocene (27 Ma) to late Pliocene (2.5 Ma)
Geol. Soc. Am., Special Paper
 , 
2002
, vol. 
357
 pg. 
v-489
 
Hamilton
C.A.
Formanowicz
D.R.
Bond
J.E.
Species delimitation and phylogeography of Aphonopelma hentzi (Araneae, Mygalomorphae, Theraphosidae): cryptic diversity in North American tarantulas
PLoS ONE.
 , 
2011
, vol. 
6
 
10
pg. 
e26207
 
Hedin
M.
Bond
J.E.
Molecular phylogenetics of the spider infraorder Mygalomorphae using nuclear rRNA genes (18S and 28S): conflict and agreement with the current system of classification
Mol. Phylogenet. Evol.
 , 
2006
, vol. 
41
 (pg. 
454
-
471
)
Hedin
M.
Carlson
D.
A new trapdoor spider species from the southern Coast Ranges of California (Mygalomorphae, Antrodiaetidae, Aliatypus coylei, sp. nov,), including consideration of mitochondrial phylogeographic structuring
Zootaxa.
 , 
2011
, vol. 
2963
 (pg. 
55
-
68
)
Hedin
M
Starrett
J
Hayashi
C.
Crossing the uncrossable: novel trans-valley biogeographic patterns revealed in the genetic history of low-dispersal mygalomorph spiders (Antrodiaetidae, Antrodiaetus) from California
Mol. Ecol.
 , 
2013
, vol. 
22
 (pg. 
508
-
526
)
Heled
J.
Drummond
A.J.
Bayesian inference of species trees from multilocus data
Mol. Biol. Evol.
 , 
2010
, vol. 
27
 (pg. 
570
-
580
)
Hendrixson
B.E.
Bond
J.E.
Testing species boundaries in the Antrodiaetus unicolor complex (Araneae: Mygalomorphae: Antrodiaetidae): “paraphyly” and cryptic diversity
Mol. Phylogenet. Evol.
 , 
2005
, vol. 
36
 (pg. 
405
-
416
)
Hendrixson
B.E.
Bond
J.E.
Molecular phylogeny and biogeography of an ancient Holarctic lineage of mygalomorph spiders (Araneae: Antrodiaetidae: Antrodiaetus)
Mol. Phylogenet. Evol.
 , 
2007
, vol. 
42
 (pg. 
738
-
755
)
Hendrixson
B.E.
DeRussy
B.M.
Hamilton
C.A.
Bond
J.E.
An exploration of species boundaries in turret-building tarantulas of the Mojave Desert (Araneae, Mygalomorphae, Theraphosidae, Aphonopelma)
Mol. Phylogenet. Evol.
 , 
2013
, vol. 
66
 (pg. 
327
-
340
)
Hill
M.L.
Dibblee
T.W.
San Andreas, Garlock, and Big Pine faults, California: a study of the character, history, and tectonic significance of their displacements
Geological Soc. America Bulletin.
 , 
1953
, vol. 
64
 (pg. 
443
-
458
)
Hudson
R.R.
Futuyma
D.
Antonovics
J.
Gene genealogies and the coalescent process
Oxford surveys in evolutionary biology
 , 
1991
New York
Oxford University Press
(pg. 
1
-
44
)
Irwin
D.E.
Phylogeographic breaks without geographic barriers to gene flow
Evolution.
 , 
2002
, vol. 
56
 (pg. 
2383
-
2394
)
Jackman
T.R.
Wake
D.B.
Evolutionary and historical analysis of protein variation in the blotched forms of salamanders of the Ensatina complex (Amphibia: Plethodontidae)
Evolution.
 , 
1994
, vol. 
48
 (pg. 
876
-
897
)
Jackson
N.D.
Austin
C.C.
The combined effects of rivers and refugia generate extreme cryptic fragmentation within the common ground skink (Scincella lateralis)
Evolution.
 , 
2010
, vol. 
64
 (pg. 
409
-
428
)
Jakobsson
M.
Rosenberg
N.A.
CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure
Bioinformatics.
 , 
2007
, vol. 
23
 (pg. 
1801
-
1806
)
Jockusch
E.L.
Martínez-Solano
I.
Hansen
R.W.
Wake
D.B.
Morphological and molecular diversification of slender salamanders (Caudata: Plethodontidae: Batrachoseps) in the southern Sierra Nevada of California with descriptions of two new species
Zootaxa.
 , 
2012
, vol. 
3190
 (pg. 
1
-
30
)
Jockusch
E.L.
Wake
D.B.
Falling apart and merging: diversification of slender salamanders (Plethodontidae: Batrachoseps) in the American West
Biol. J. Linn. Soc.
 , 
2002
, vol. 
76
 (pg. 
361
-
391
)
Katoh
K.
Kuma
K.
Toh
H.
Miyata
T.
MAFFT version 5: improvement in accuracy of multiple sequence alignment
Nucleic Acids Res.
 , 
2005
, vol. 
33
 (pg. 
511
-
518
)
Keith
R.
Hedin
M.
Extreme mitochondrial population subdivision in southern Appalachian paleoendemic spiders (Araneae: Hypochilidae: Hypochilus), with implications for species delimitation
J. Arachnology.
 , 
2012
, vol. 
40
 (pg. 
167
-
181
)
Kingman
J.F.C.
The coalescent
Stochastic Processes and their Applications.
 , 
1982
, vol. 
13
 (pg. 
235
-
248
)
Knowles
L.L.
Carstens
B.C.
Delimiting species without monophyletic gene trees
Syst. Biol.
 , 
2007
, vol. 
56
 (pg. 
887
-
895
)
Kubatko
L.S.
Carstens
B.C.
Knowles
L.L.
STEM: species tree estimation using maximum likelihood for gene trees under coalescence
Bioinformatics.
 , 
2009
, vol. 
25
 (pg. 
971
-
973
)
Kubatko
L.
Gibbs
H.
Bloomquist
E.
Inferring species-level phylogenies and taxonomic distinctiveness using multilocus data in Sistrurus rattlesnakes
Syst. Biol.
 , 
2011
, vol. 
60
 (pg. 
393
-
409
)
Leaché
A.D.
Fujita
M.K.
Bayesian species delimitation in West African forest geckos (Hemidactylus fasciatus)
Proc. Royal Soc. B.
 , 
2010
, vol. 
277
 (pg. 
3071
-
3077
)
Lohse
K.
Can mtDNA barcodes be used to delimit species? A response to Pons et al. (2006)
Syst. Biol.
 , 
2009
, vol. 
58
 (pg. 
439
-
442
)
Maddison
D.R.
Maddison
W.P.
MacClade 4: Analysis of phylogeny and character evolution. Version 4.08a
2005
 
Available from: http://macclade.org
McGuire
G.
Wright
G.
Prentice
M.J.
A graphical method for detecting recombination in phylogenetic data sets
Mol. Biol. Evol.
 , 
1997
, vol. 
14
 (pg. 
1125
-
1131
)
Milne
I.
Lindner
D.
Bayer
M.
Husmeier
D.
McGuire
G.
Marshall
D.F.
Wright
F.
TOPALi v2: a rich graphical interface for evolutionary analyses of multiple alignments on HPC clusters and multi-core desktops
Bioinformatics.
 , 
2009
, vol. 
25
 (pg. 
126
-
127
)
Milne
I.
Wright
F.
Rowe
G.
Marshall
D.F.
Husmeier
D.
McGuire
G.
TOPALi: software for automatic identification of recombinant sequences within DNA multiple alignments
Bioinformatics.
 , 
2004
, vol. 
20
 (pg. 
1806
-
1807
)
Myers
N.
Mittermeier
R.A.
Mittermeier
C.G.
da Fonseca
G.A.B.
Kent
J.
Biodiversity hotspots for conservation priorities
Nature.
 , 
2000
, vol. 
403
 (pg. 
853
-
858
)
Niemiller
M.L.
Near
T.J.
Fitzpatrick
B.M.
Delimiting species using multilocus data: diagnosing cryptic diversity in the southern cavefish, Typhlichthys subterraneus (Teleostei: Amblyopsidae)
Evolution.
 , 
2012
, vol. 
66
 (pg. 
846
-
866
)
O'Meara
B.C.
New heuristic methods for joint species delimitation and species tree inference
Syst. Biol.
 , 
2010
, vol. 
59
 (pg. 
59
-
73
)
Papadopoulou
A.
Bergsten
J.
Fujisawa
T.
Monaghan
M.T.
Barraclough
T.G.
Vogler
A.P.
Speciation and DNA barcodes: testing the effects of dispersal on the formation of discrete sequence clusters
Phil. Trans. R. Soc. B.
 , 
2008
, vol. 
363
 (pg. 
2987
-
2996
)
Parham
J.F.
Papenfuss
T.J.
High genetic diversity among fossorial lizard populations (Anniella pulchra) in a rapidly developing landscape (Central California)
Conserv. Genet.
 , 
2009
, vol. 
10
 (pg. 
169
-
176
)
Polihronakis
M.
Caterino
M.S.
Contrasting patterns of phylogeographic relationships in sympatric sister species of ironclad beetles (Zopheridae: Phloeodes spp.) in California's Transverse Ranges
BMC Evol. Biol.
 , 
2010
, vol. 
10
 pg. 
195
 
Pons
J.
Barraclough
T.G.
Gomez-Zurita
J.
Cardoso
A.
Duran
D.P.
Hazell
S.
Kamoun
S.
Sumlin
W.D.
Vogler
A.P.
Sequence-based species delimitation for the DNA taxonomy of undescribed insects
Syst. Biol.
 , 
2006
, vol. 
55
 (pg. 
595
-
609
)
Pons
J.
Fujisawa
T.
Claridge
E.M.
Savill
R.A.
Barraclough
T.G.
Vogler
A.P.
Deep mtDNA subdivision within Linnean species in an endemic radiation of tiger beetles from New Zealand (genus Neocicindela)
Mol. Phylogenet. Evol.
 , 
2011
, vol. 
59
 (pg. 
251
-
262
)
Posada
D.
jModelTest: Phylogenetic Model Averaging
Mol. Biol. Evol.
 , 
2008
, vol. 
25
 (pg. 
1253
-
1256
)
Pritchard
J.K.
Stephens
M.
Donnelly
P.
Inference of population structure using multilocus genotype data
Genetics.
 , 
2000
, vol. 
155
 (pg. 
945
-
959
)
Pritchard
J.K.
Wen
X.
Falush
D.
Documentation for structure software: Version 2.3
2010
 
Available from: http://pritch.bsd.uchicago.edu/structure.html.1-38 (last accessed October 8, 2012)
Rannala
B.
Yang
Z.
Bayes estimation of species divergence times and ancestral population sizes using DNA sequences from multiple loci
Genetics.
 , 
2003
, vol. 
164
 (pg. 
1645
-
1656
)
Raven
R.J.
The spider infraorder Mygalomorphae (Araneae): cladistics and systematics
Bull. Am. Mus. Nat. Hist.
 , 
1985
, vol. 
182
 (pg. 
1
-
175
)
Reid
N.M.
Carstens
B.C.
Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yule-coalescent model
BMC Evol. Biol.
 , 
2012
, vol. 
12
 pg. 
196
 
Rittmeyer
E.N.
Austin
C.C.
The effects of sampling on delimiting species from multi-locus sequence data
Mol. Phylogenet. Evol.
 , 
2012
, vol. 
65
 (pg. 
451
-
463
)
Rosenberg
N.A.
Distruct: a program for the graphical display of population structure
Mol. Ecol. Notes.
 , 
2004
, vol. 
4
 (pg. 
137
-
138
)
Rosenberg
N.A.
Statistical tests for taxonomic distinctiveness from observations of monophyly
Evolution.
 , 
2007
, vol. 
61
 (pg. 
317
-
323
)
Sáez
A.G.
Lozano
E.
Body doubles
Nature.
 , 
2009
, vol. 
433
 pg. 
111
 
Satler
J.D.
Species trees and species delimitation in the California trapdoor spider genus Aliatypus (Mygalomorphae, Antrodiaetidae)
2011
 
Master's thesis, San Diego State University
Satler
J.D.
Starrett
J.
Hayashi
C.Y.
Hedin
M.
Inferring species trees from gene trees in a radiation of California trapdoor spiders (Araneae, Antrodiaetidae, Aliatypus)
PLoS ONE.
 , 
2011
, vol. 
6
 
9
pg. 
e25355
 
Stamatakis
A.
RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models
Bioinformatics.
 , 
2006
, vol. 
22
 (pg. 
2688
-
2690
)
Stamatakis
A.
Hoover
P.
Rougemont
J.
A rapid bootstrap algorithm for the RAxML Web Servers
Syst. Biol.
 , 
2008
, vol. 
57
 (pg. 
758
-
771
)
Starrett
J.
Hedin
M.
Multilocus genealogies reveal multiple cryptic species and biogeographical complexity in the California turret spider Antrodiaetus riversi (Mygalomorphae, Antrodiaetidae)
Mol. Ecol.
 , 
2007
, vol. 
16
 (pg. 
583
-
604
)
Stephens
M.
Donnelly
P.
A comparison of Bayesian methods for haplotype reconstruction from population genotype data
Am. J. Hum. Genet.
 , 
2003
, vol. 
73
 (pg. 
1162
-
1169
)
Stephens
M.
Smith
N.
Donnelly
P.
A new statistical method for haplotype reconstruction from population data
Am. J. Hum. Genet.
 , 
2001
, vol. 
68
 (pg. 
978
-
989
)
Stockman
A.K.
Bond
J.E.
Delimiting cohesion species: extreme population structuring and the role of ecological interchangeability
Mol. Ecol.
 , 
2007
, vol. 
16
 (pg. 
3374
-
3392
)
Swofford
D.L.
PAUP*: Phylogenetic analysis using parsimony (*and Other Methods). Version 4
2002
Sunderland (MA)
Sinauer Associates
Wakabayashi
J.
Sawyer
T.L.
Stream incision, tectonics, uplift, and evolution of topography of the Sierra Nevada, California
J. Geology.
 , 
2001
, vol. 
109
 (pg. 
539
-
562
)
Weisrock
D.W.
Rasoloarison
R.M.
Fiorentino
I.
Ralison
J.M.
Goodman
S.M.
Kappeler
P.M.
Yoder
A.D.
Delimiting species without nuclear monophyly in Madagascar's mouse lemurs
PLoS ONE.
 , 
2010
, vol. 
5
 
3
pg. 
e9883
 
Wiens
J.J.
Species delimitation: new approaches for discovering diversity
Syst. Biol.
 , 
2007
, vol. 
56
 (pg. 
875
-
878
)
Yang
Z.
Rannala
B.
Bayesian species delimitation using multilocus sequence data
Proc. Natl. Acad. Sci. USA.
 , 
2010
, vol. 
107
 (pg. 
9264
-
9269
)
Zhang
C.
Zhang
D.X.
Zhu
T.
Yang
Z.
Evaluation of a Bayesian coalescent method of species delimitation
Syst. Biol.
 , 
2011
, vol. 
60
 (pg. 
747
-
761
)

Author notes

Associate Editor: Robb Brumfield