Integrative Taxonomy and Species Delimitation in Harvestmen: A Revision of the Western North American Genus Sclerobunus (Opiliones: Laniatores: Travunioidea)

Alpha taxonomy, and specifically the delimitation of species, is becoming increasingly objective and integrative. The use of coalescent-based methods applied to genetic data is providing new tools for the discovery and delimitation of species. Here, we use an integrative approach via a combination of discovery-based multivariate morphological analyses to detect potential new species. These potential species are then used as a priori species in hypothesis-driven validation analyses with genetic data. This research focuses on the harvestmen genus Sclerobunus found throughout the mountainous regions of western North America. Based on our analyses, we conduct a revision of Sclerobunus resulting in synonymy of Cyptobunus with Sclerobunus including transfer of S. cavicolens comb. nov. and elevation of both subspecies of S. ungulatus: S. ungulatus comb. nov. and S. madhousensis comb. nov., stat. nov. The three subspecies of S. robustus are elevated, S. robustus, S. glorietus stat. nov., and S. idahoensis stat. nov. Additionally, five new species of Sclerobunus are described from New Mexico and Colorado, including S. jemez sp. nov., S. klomax sp. nov., S. skywalkeri sp. nov., S. speoventus sp. nov., and S. steinmanni sp. nov. Several of the newly described species are single-cave endemics, and our findings suggest that further exploration of western North American cave habitats will likely yield additional new species.


Introduction
In order to accurately identify and describe the diversity of life, species must be delimited operationally [1].With increasing stability of the conceptual definition of a species [2], the operational delimitation of species becomes an easier task [3,4].It is becoming increasingly clear that taxonomy, and species delimitation in particular, must become an integrative field utilizing multiple independent lines of evidence [3,5] as any single dataset may not accurately reflect species limits and relationships (e.g., due to morphological convergence, gene tree/species tree discordance, etc.).Proper integrative species delimitation should incorporate a combination of discovery-based approaches without a priori species hypotheses and, following this, hypothesis-driven approaches, which test those groups that are delimited based on discovery methods [3,4,6].Additionally, in order to decrease the degree of subjectivity that is common in traditional alphataxonomic practices, species delimitation methods are moving towards increasing objectivity [4].For example, recently developed methods accomplish this by relying on statistical thresholds to delimit species in multivariate morphological space [7], or posterior probability distributions to determine support for the presence of a species tree node [8].
The discovery and documentation of biodiversity at the species level is especially important for conservation biology [9], particularly in taxa that are relatively poorly studied (either due to sheer diversity of species and/or proportionally few taxonomists).With the almost ubiquitous incorporation of genetic data in systematics, the discovery and documentation of new and cryptic species is increasing [10], seemingly correlated with the use of coalescent-based species delimitation methods [4,11].Recent studies have shown the potential for these methods in the discovery of new and cryptic species in taxa that are either poorly known, morphologically conserved, show low vagility, and/or have narrow ecological limits [12][13][14][15][16].The order Opiliones (commonly called harvestmen) is a diverse group of arachnids with over 6500 described species [17], distributed on every continent except Antarctica.Despite being a fairly diverse group (i.e., more described species than mammals), harvestmen are relatively poorly studied.The vast majority of harvestmen species show low vagility with high ecological constraints; attributes that have resulted in the discovery of many new, and sometimes cryptic, species when studied using modern molecular approaches [18][19][20][21][22][23][24].
In this study we focus on the harvestmen genera Sclerobunus and Cyptobunus, which are broadly distributed throughout the mountainous regions of western North America.These genera are ecologically limited to moist, dark microhabitats, typically found under logs and rocks in high-elevation forests or in caves.The first species of Sclerobunus was described from Colorado by A. S. Packard in 1877 as Scotolemon robustum [25].Later, Banks [26] described the new genus Sclerobunus, and included two species: Packards' Scotolemon robustum, renamed Sclerobunus robustus (Packard, 1877), and Sclerobunus brunneus.This second species was later transferred to the genus Paranonychus Briggs, 1971.Banks' [26] description of Sclerobunus robustus was based on specimens from Colorado, Utah (originally described by Packard), and Washington.Roewer [27] then described a new species, Sclerobunus parvus, from Vancouver Island, British Columbia.Many years later in a revision of the North American ''Triaenonychidae'', Briggs [28] included two species of Sclerobunus.The first species, Sclerobunus nondimorphicus was described as new from Oregon, Washington, and southern British Columbia (this species includes Banks' [26] specimen of S. robustus from Washington).The second species, S. robustus, was redescribed and recorded with a very broad distribution throughout the American southwest and Rocky Mountains.Within this species, Briggs described three subspecies: S. robustus robustus from Arizona and New Mexico, S. r. glorietus from a single locality in northern New Mexico, and S. r. idahoensis from Idaho.More recently, Shear and Derkarabetian [29], upon examination of the type specimens of S. parvus, synonymized this taxon with Paranonychus brunneus (Banks, 1893).Currently, Sclerobunus includes two species: Sclerobunus nondimorphicus and Sclerobunus robustus, itself including three subspecies.
The taxonomic and phylogenetic relationship between Sclerobunus and its sister genus Cyptobunus remains uncertain.All Cyptobunus are known only from cave habitats and are highly troglomorphic.Banks [30] described the new genus and species, Cyptobunus cavicolus, from Lewis and Clark Caverns, Montana based on a juvenile specimen (as mentioned by Roewer, [31]) and remarked that this taxon ''is but a cavernicolous adaptation of Sclerobunus'' [30].With adult specimens, Crosby and Bishop [32] synonymized Cyptobunus with Sclerobunus and this synonymy was acknowledged by Goodnight and Goodnight [33].Later, Briggs [28] described the new species C. ungulatus, with two subspecies: C. u. ungulatus from Model Cave in Great Basin National Park, Nevada and C. u. madhousensis from North Madhouse Cave near Provo, Utah.Briggs retained separate genera noting differences in several somatic characters, but also remarked on the general similarity in male genitalia.Historically, Sclerobunus are only known from epigean habitats.However, recent fieldwork has uncovered many new populations of troglomorphic Sclerobunus from cave and talus habitats, extending the known habitat preference for this taxon [21].Sclerobunus populations found in caves show varying degrees of troglomorphy including decreased pigmentation, elongation of appendages, and attenuation of spines.Based on this, it is possible that Cyptobunus taxa merely represent highly derived Sclerobunus.A recent genetic analysis by Derkarabetian et al. [21] recovered Cyptobunus phylogenetically nested within Sclerobunus and showed that troglomorphy evolved independently at least three times; however, this was only a twogene study.Although there remains uncertainty in the relationship between Sclerobunus and Cyptobunus, their monophyly relative to other closely related genera is certain and reflected in penis morphology [28] and genetic data [21].
The morphological and genetic analyses of Derkarabetian et al. [21] suggested that there are several new morphologically distinct, genetically divergent species within Sclerobunus.Here, we use discovery-based methods to identify multivariate clusters using morphometric data, which are then treated as putative species in species tree analyses using eight newly developed nuclear loci.Following this, we use a hypothesis-driven approach to test alternative species delimitation hypotheses using the multigenic nuclear dataset.Based on results of our species delimitation analyses, we conduct a revision of Sclerobunus, including Cyptobunus.

Taxon Sampling
Fieldwork was conducted in the summer months of 2006, 2007, and 2008 throughout the mountainous regions of western North America.Samples representing all species and subspecies of Sclerobunus and Cyptobunus were collected for analyses (TABLE S1).Fresh specimens were collected from type localities for all species and subspecies, except S. nondimorphicus, S. r. robustus, and S. r. idahoensis.For these taxa, type localities could not be accessed either due to vague locality records or unsuitable habitat, and fresh samples were collected from as near the type locality as possible.All specimens were collected by hand; those destined for molecular analyses were placed in 100% EtOH and stored at 2 80uC, while specimens to be used in morphological analyses were stored in 80% EtOH.
Samples from caves in Great Basin National Park, Nevada were collected under a Scientific Research Collecting Permit from the U.S. National Park Service permit granted to the authors (GRBA-2007-SCI-0008). Permission to collect in Lewis and Clark Caverns, Montana was given by the Montana State Park staff and samples from Cave of the Winds, Colorado were collected with permission from the management staff.In general, sample sizes from cave populations were limited either purposefully due to conservation concerns, or due to rarity of individuals.To increase the sample size for some cave populations, recently collected specimens from Cave of the Winds and newly discovered Mallory Cave, Colorado specimens from the Denver Museum of Nature and Science were included.In addition, specimens collected from Terrero Cave, New Mexico held in the California Academy of Sciences collection were included in morphological analyses.

Morphometric Analyses
The morphometric dataset included 18 linear measurements taken using an Olympus SZX12 microscope equipped with an ocular micrometer (TABLE S2).The characters chosen included the length and width of the scute, chelicerae, and genital operculum, height and width of the ocularium, height of the pedipalpal femur, and length of the second leg.For some specimens, certain characters could not be measured, for example, missing genital opercula or legs (total of 2.4% missing data).In these cases, missing values were estimated with SPSS v. 20 (IBM Corp.) using the Multiple Imputation function with default settings and 10 iterations.The average of 10 iterations was used as the missing value in the final analyses.Because of sexual dimorphism, analyses were conducted on separate male and female datasets.
We utilize an algorithmic approach described by Ezard et al. [7] that relies on statistical thresholds at multiple steps to objectively define species clusters in multivariate morphospace.This method is used as a discovery-based approach to identify putative species, which we then further tested with genetic data.The Ezard et al. method has several advantages over traditional PCA [7].First, traditional PCA using a covariance matrix assume normally distributed data, which is usually not met in biological datasets that contain multiple species.Our data do not meet this assumption as multiple species are included and 7 of 18 characters are not normally distributed (Shapiro-Wilks test; results not shown).Second, analyses computed via correlation matrices can be negatively influenced by outliers, which can skew orientation of the axes.To circumvent these issues, this method uses robust scale estimators that perform better with non-normal data.Use of the median to scale the data de-emphasizes the effects of any outliers, and is more useful in discriminating incipient species.The Bayesian Information Criteria (BIC) is used to choose the optimal model among possible models that vary in the size, orientation, and number of clusters.Importantly, this method can be used as a null model of morphological homogeneity.Finally, the automation and use of multiple thresholds increases reproducibility and similarity in interpretations across datasets and taxa.We implemented the Ezard et al. method using the default protocol with the k-value set to 18 (total number of characters measured in this study) using R 3.0.2(R Core Team, 2013).In addition, standard principal components (PCA) and discriminant function (DFA) analyses were conducted (see FILE S1 for details); however, we do not use these results in our species delimitation decisions.The PCA and DFA analyses and results are included only as a simple results comparison of the recently developed method of Ezard et al. [7].
We did not include holotype specimens in morphometric analyses for two main reasons: First, all species of Sclerobunus are allopatric in distribution, with one exception at Taos Ski Valley, New Mexico, which is discussed below and at length in Derkarabetian et al. [21].As such, we argue that inclusion of fresh specimens from the type locality is just as valuable as including holotypes.The allopatric distribution of Sclerobunus species eliminates the possibility that type series specimens are heterospecific or that specimens collected from type localities represent multiple species.Second, holotype specimens range from 40-100+ years old and may not be in suitable condition for morphometric analyses, particularly the more fragile cave species.

Molecular Data Collection
Based on a comparison of transcriptomic data (see FILE S2), we developed primers for 61 exon and untranslated (UTR) nuclear gene regions, 23 of which successfully amplified for a small panel of 4 Sclerobunus taxa (2 S. robustus, 1 S. glorietus, and 1 S. nondimorphicus).Of these, 8 gene regions were pursued further based on informativeness and minimal levels of heterozygosity.It has been shown that species limits can be successfully determined with $5 loci and moderate sample sizes (5 specimens per species), even with short node depths and gene flow [34,35].The 8 nuclear loci included six exonic regions and two 39UTR regions, and were sequenced for 43 total samples of Sclerobunus and Cyptobunus (TABLE S3) using the primers and PCR conditions detailed in FILE S2.Bi-directional Sanger reads were assembled into contigs, edited, and unambiguously aligned using Geneious Pro 6 (http:// www.geneious.com).Diversity statistics and tests for recombination and neutrality for each locus were calculated using DnaSP v 5 [36].Haplotypes were determined using Phase 2.1.1 [37].Both of the 39 UTR loci contained some sequences that could not be resolved to haplotypes.For these two loci, unresolved sequences were removed prior to tests of recombination and neutrality.GenBank accession numbers are provided in TABLE S3.Matrices with phased sequences are available from the Dryad Digital Repository: doi:10.5061/dryad.hj6r1.

Species Tree Analyses
The clusters recovered from morphological, discovery-based analyses were used in conjunction with clades recovered in the COI mtDNA analyses of Derkarabetian et al. [22] to define a priori species hypotheses for species trees analyses.In previous analyses, three well-supported, deeply divergent species groups are recovered: a clade from the Pacific Northwest including S. nondimorphicus and S. r. idahoensis, a clade including all Cyptobunus, and a clade including all southwestern Sclerobunus.The monophyly and distinctiveness of these three species groups is also reflected in genitalic morphology [21,28].Morphometric clusters may include multiple species; in these clusters, species group assignments can help resolve cases of clustering due to broad ''morphospace convergence'' between taxa from different species groups, or cases of potential conspecificity of putatively different species within the same species groups (FIGURE 1).In instances where clusters contain specimens from two different species groups, those divergent taxa will be considered separate putative species.Clusters with potentially different taxa from the same species group will rely on validation analyses.
All genetic analyses conducted use only the 8 nuclear genes sequenced here; we do not reanalyze or include previously collected nuclear [21] or mitochondrial data [22] (see Discussion).Different populations of the geographically widespread subspecies S. robustus robustus represent potential cryptic species [22], but because our current nuclear genetic sampling for this taxon is limited, this subspecies was treated as a single taxon in this paper.
To determine rooting and ingroup polarity, we conducted preliminary analyses including only the six exons with Theromaster as the outgroup.The 39 UTR gene regions were either not present in the sequenced Theromaster transcriptome, or could not be unambiguously aligned to the ingroup.Following this, all subsequent species tree analyses were conducted on the complete 8-gene set with Theromaster sequences removed.Models of evolution and optimal partitioning strategy were determined with PartitionFinder [38] using the BIC to choose optimal models.The species tree was reconstructed using the *BEAST algorithm [39].Analyses were run for 200 million generations logging every 1000 generations, implemented in BEAST 1.8 [40] using the CIPRES portal (http://www.phylo.org/sub_sections/portal/).Each resulting run was checked for stationarity and for ESS values above 200 with Tracer 1.5 [41].Species trees were reconstructed with TreeAnnotator (http://beast.bio.ed.ac.uk/TreeAnnotator) from 180,000 trees (10% burnin removal).Individual gene trees were reconstructed using RAxML version 8 [42] with the GTRGAMMA model and 1000 rapid bootstrap replicates.All gene and species trees are available from the Dryad Digital Repository: doi:10.5061/dryad.hj6r1.

Species Delimitation Analyses
We use the general lineage concept (as defined in [2]) as our theoretical concept of a species and utilize multiple analyses to provide supporting criteria with which we can operationally delimit general lineage species.Groups of samples showing concordance across different data sets and analyses (identified as different morphometric clusters, highly supported mtDNA clades, validated nuclear clades, and distinguishable via general morphology) are considered species [43].A general workflow for the species delimitation decision-making process is provided in FIGURE 1.
The Bayesian program BPP [8,44] was used to test specific species delimitation hypotheses.This program conducts multilocus, coalescent-based analyses requiring a guide tree and specification of two priors controlling population size and divergence time.The program incorporates a reversible-jump Markov chain Monte Carlo algorithm (rjMCMC) that explores all possible species delimitation models, ultimately providing an assessment of the probability of a node being present.Although BPP is known to recover a higher number of partitions relative to other programs, particularly for micro-allopatric taxa with deep genetic divergences, we believe that BPP used in conjunction with discovery-based analyses can lead to informed, conservative decisions about species delimitations [16].Specific hypotheses tested included the following: S. nondimorphicus+S.robustus idahoensis, C. ungulatus ungulatus+C.u. madhousensis, Cave of the Winds+Mallory Cave, and within the S. r. glorietus complex: S. r. glorietus (Glorieta Canyon, New Mexico)+S.r. glorietus (Taos Ski Valley), southern S. r. glorietus+S.r. glorietus (Glorieta Canyon), and the syntopic surface and troglomorphic populations of S. r. glorietus from Taos Ski Valley.The species tree recovered from *BEAST analyses based on morphological clusters was used as the input guide tree.In BPP, the rjMCMC species delimitation method [44] was used with algorithm 0 (e = 10), rates were allowed to vary among loci (locusrate = 1), gaps and ambiguous columns were removed (cleandata = 1), and the analysis was set for automatic fine-tune adjustments.Four different combinations of theta (h) and tau (t) priors were used to span a diversity of possible population sizes and divergence times.Each prior combination was run twice to check convergence of runs and proper mixing.Analyses were run for 200,000 generations, sampling every 5 generations with 20,000 burnin.Posterior probabilities greater than 95 for the presence of a node are considered supported species delimitations.
Using the gene trees estimated with RAxML, we calculated the genealogical sorting index (gsi) for each gene and the ensemble genealogical sorting index (gsi T ) [45] using the gsi website (http:// www.genealogicalsorting.org) with 10,000 replicates.The gsi statistic assesses the level of genealogical exclusivity of a group where values range from 0 (a random distribution of haplotypes) to 1 (exclusive ancestry, monophyly), and the ensemble analysis (gsi T ) incorporates uncertainty of relationships and integrates across all gene trees.The groups assessed consisted of the putative species identified via discovery-based analyses (FIGURE 1).
Based on multivariate morphological clustering, there is evidence that the S. r. glorietus subspecies comprises more than one species; however, species limits are uncertain due to sparse geographic sampling.It is important to note that the subspecies S. r. glorietus was described from only a single locality, Glorieta Canyon, New Mexico.As such, and because one of the highly troglomorphic populations is recovered within the S. r. glorietus complex, we simultaneously estimated species limits and species trees for this complex using the Bayes factor delimitation (BFD) method of Grummer et al. [46], testing several different species limit combinations (FILE S3).This method, based on Bayes factors, is beneficial in cases where species limits are uncertain as topological uncertainty is accounted for, and a priori species relationships are not required [46].Additionally, we ran a *BEAST analysis in which each individual S. r. glorietus population was considered a separate species, then tested each species (population) using BPP (FILE S3).

Taxonomy
Digital images were taken using a Visionary Digital system (http://www.visionarydigital.com).Several images were taken at different focal planes and combined using Zerene Stacker (Zerene Systems LLC).Scanning electron microscopy was conducted using methods outlined in Derkarabetian et al. [22].Leg and pedipalp illustrations were done using a drawing tube attached to an Olympus SZX12 microscope with subsequent tracing in Adobe Illustrator CS5.Penis illustrations were traced from digital images.Digital images have been deposited to Morphbank; image ID numbers are included in the figure legends.A Google Earth kmz file with all known localities including sites from this study, previous studies, and geo-referenced museum records for each species is included as FILE S4.All relevant information will be integrated into Encyclopedia of Life webpages (http://www.eol.org) for each species.

Nomenclatural Acts
The electronic edition of this article conforms to the requirements of the amended International Code of Zoological Nomenclature, and hence the new names contained herein are available under that Code from the electronic edition of this article.This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the ICZN.The ZooBank LSIDs (Life Science Identifiers) can be resolved and the associated information viewed through any standard web browser by appending the LSID to the prefix ''http://zoobank.org/''.The LSID for this publication is: urn:lsid:zoobank.org:pub:46DB0A1A-E317-4CD3-A124-53E97418C4B0.The electronic edition of this work was published in a journal with an ISSN, and has been archived and is available from the following digital repositories: PubMed Central, LOCKSS.

Taxon Sampling
All collecting localities from recent fieldwork for this study, including specimens used in morphological and/or genetic analyses, are shown in FIGURES 2 and 3.

Morphometric Analyses
The results of the traditional PCA and DFA are presented in FILE S1.Here, we present the results of the multivariate clustering method of Ezard et al. [7].The analysis of male specimens retained 3 components (FIGURE 4, FILE S5) and selected the EEE model (ellipsoidal, equal volume, shape and orientation) with 9 clusters as optimal.Clusters correspond to ungulatus, Mallory Cave, Cave of the Winds+madhousensis, cavicolens+Terrero Cave, northern glorietus (Taos Ski Valley and Glorieta Canyon), southern glorietus+robustus from Bradford Canyon, New Mexico, and three separate clusters largely corresponding to nondimorphicus, idahoensis, and robustus.The recovery of the Bradford Canyon robustus specimens with the generally smaller-bodied glorietus subspecies is not surprising considering Bradford Canyon specimens have the smallest body size of any robustus population known.The nondimorphicus cluster included 2 specimens of idahoensis, one of which is the largest idahoensis sampled (Goose Creek, Montana).Two specimens of ungulatus were considered significant outliers but were not removed from analyses.
Analysis of female specimens retained 2 components (FILE S5) and the optimal model chosen was the EEE model with 4 clusters corresponding to glorietus (northern and southern), robustus+ idahoensis+nondimorphicus, Cave of the Winds+Mallory Cave, and cavicolens+Terrero Cave+Taos Ski Valley troglomorph.Although we include female morphometric analyses, we mainly rely on results of the male analyses for downstream decisions, as males are more variable and easily distinguished.

Molecular Data Collection
GenBank accession numbers and the final genetic dataset are shown in TABLE S3 and alignment statistics are provided in TABLE 1.In total, original data for 8 nuclear gene regions were gathered from 36 Sclerobunus and 7 Cyptobunus samples.For each locus, no more than 5 samples were missing data (average of 2 sequences missing per gene), and no single specimen had more than one missing gene region (in total, 4.3% missing data).

Species Tree Analyses
Model and partition selection via BIC resulted in single partitions for each gene region (TABLE 1).Individual RAxML gene trees are included in FILE S6.Preliminary rooting analyses recovered northwestern Sclerobunus (nondimorphicus and idahoensis) as sister to all remaining lineages.Nodes were considered highly supported if the posterior probability values exceeded 95.In the species tree (FIGURE 5), all nodes were resolved with 100% posterior probability, with the exception of the node containing the southwestern Sclerobunus and the node between glorietus and klomax.The genus Cyptobunus is phylogenetically nested within Sclerobunus with 100% posterior probability.Therefore, as the generic name Sclerobunus has precedence [26], we synonymize Cyptobunus with Sclerobunus, and all species of Cyptobunus are transferred to Sclerobunus (see below).
The fact that the species tree is well resolved, with the exception of two nodes, confirms the presence of three species groups recovered in previous analyses [21,22] and helps resolve those morphological clusters containing more than one putative species.In this regard, morphometric clusters containing species from different species groups are easily separated.In the male analysis (FIGURE 4) a single specimen of idahoensis was recovered within a cluster comprised of mostly robustus specimens, a cluster included all cavicolens and Terrero Cave specimens, and another cluster included all Cave of the Winds and madhousensis specimens.Similarly in the female analysis, a single cluster included the Taos Ski Valley troglomorph, Terrero Cave, and cavicolens specimens (FILE S5).As these clusters contain taxa from different species groups, they do not represent potential conspecific taxa clusters, but instead represent cases of morphological convergence across distantly related taxa.The two clusters including multiple taxa from the same species group (nondimorphicus/idahoensis and robustus/southern glorietus) were not delimited at this step and relied upon validation methods.Taken together, the morphometric and species trees analyses, coupled with previous genetic analyses, conservatively support at least eight putative species corresponding to nondimorphicus/idahoensis, cavicolens, ungulatus (two subspecies), Cave of the Winds/ Mallory Cave, robustus, southern glorietus, northern glorietus, and Terrero Cave/Taos Ski Valley troglomorph.

Species Delimitation
Results of the BPP species delimitation analyses are shown in FIGURE 5.All species, including all putative new species and morphometric clusters containing multiple putative species from the same species group, were delimited with 100% posterior  probability in all runs under all prior combinations.The two subspecies of ungulatus are an exception.Morphological clustering methods showed a clear distinction between the two subspecies of ungulatus.However, in the delimitation analyses with higher values of theta, h,C(2, 10), indicative of large effective population sizes, the presence of the node was not supported.For cave species, which generally have small populations sizes, the analyses with low values of theta, h,C(2, 1000), may be more appropriate priors for delimitation.Additionally, these isolated cave populations are separated by over 240 kilometers of arid, uninhabitable terrain where dispersal and gene flow is impossible.This extreme isolation clearly supports these two subspecies as separately evolving lineages.As such, we elevate both subspecies to full species: S. ungulatus comb.nov.and S. madhousensis comb.nov., stat.nov.Morphometric clustering methods (FIGURE 4) and genetic delimitation analyses (FIGURE 5) distinguished the new species S. speoventus stat.nov.(Cave of the Winds) from S. steinmanni stat.nov.(Mallory Cave).Although the S. madhousensis specimens clustered with S. speoventus, these two taxa are from different species groups (FIGURE 4) and are qualitatively morphologically different.Results of gsi analyses show exclusive ancestry for all species in at least 5 of 8 loci, with the exception of S. robustus, which was exclusive for only 3 of 8 loci and S. nondimorphicus, S. idahoensis, and S. glorietus, which were nonexclusive for all loci (TABLE 2).The gsi T is greater than 0.5 for all species and greater than 0.8 for 8 of 11 species, demonstrating that a high degree of genealogical sorting has occurred in these species and they are progressing towards exclusive ancestry.
A summary of all species delimitation analyses in the context of an integrative taxonomic framework is shown in FIGURE 6. Results presented here support the presence of five new species and the elevation of all three subspecies of S. robustus: S. robustus, S. glorietus stat.nov., and S. idahoensis stat.nov.The species S. idahoensis, a former subspecies of S. robustus, is not sister to other S. robustus subspecies, and was successfully delimited from S. nondimorphicus (see Taxonomy for details).Within the S. glorietus complex three species were delimited based on morphology (FIGURE 4): a northern clade (S. glorietus stat.nov.), a southern clade (S. skywalkeri sp.nov.), and the highly troglomorphic population from Taos Ski Valley (S. klomax sp.nov.).Although BFD analyses slightly favor including the Glorieta Canyon population with the southern S. glorietus clade (FILE S3), the support is not strong (both BF and posterior probability).Species tree analyses in which each population of S. glorietus are considered putative species resulted in little support for relationships within the complex, except for S. skywalkeri (FILE S3).
We also describe a new troglomorphic species from Terrero Cave, S. jemez sp.nov.Genetic data for this cave population were not available for this study and access to the cave to collect fresh specimens could not be granted, as the cave is culturally significant to the Jemez Pueblo.Although the S. jemez specimens clustered with the S. cavicolens specimens in male and female morphometric analyses (FIGURE 4, FILE S5), these taxa are from different species groups and easily distinguished morphologically.In the female morphometric analyses (FILE S5), the S. jemez specimens additionally clustered with the S. klomax specimens.These two species can be differentiated based on several characters (see below).Geographic evidence also supports this distinction as these caves are separated by over 90 km.
Description.Body length 1.68-3.68,length of scute 1.58-2.69.Body surface structure microgranulate-rivulose, integument color reddish/brown, yellow/orange to yellow in troglomorphic species, with or without black pigment.Anterior margin of scute with 1-4 tubercles on each shoulder; highly troglomorphic species without tubercles.OC recessed from anterior margin.Pedipalps strongly armed with many spine-like tubercles each bearing elongate setae subapically; first 3 proximal ventral spines of femur larger followed by 3 or more smaller spines; patella with 2 medial SBTs at distal margin; tibia with row of 4-6 medial spines and row of 4-7 lateral spines; tarsus with 4 large lateral and medial spines.LI femur with ventral row of 3 or more SBTs.
Sexual Dimorphism.As noted by Briggs [28], the PF is noticeably thicker in males of all Sclerobunus (except S. nondimorphicus).We note two additional sexually dimorphic characters.First, in males the CII lobes found along the ventral midline possess 2-3 pairs of apophyses (FILE S7d), while female CII lobes do not possess apophyses.Second, males possess a small distal process on the ventral surface of the pedipalpal tarsus (FILE S7e).Males of the S. nondimorphicus group may additionally possess many smaller asetose tubercles in 2 rows extending proximally.
Distribution.Mountainous regions of western North America (Arizona, New Mexico, Colorado, Utah, Nevada, Montana, Idaho, Oregon, Washington, and British Columbia).
Distribution.Isolated caves and cave systems in Montana, Nevada, and Utah.
Sclerobunus ungulatus (Briggs, 1971) Diagnosis.S. ungulatus differentiated from S. cavicolens by a complete lack of pigment, lack of prongs on hind claws, pedipalpal tibia with 5 or more medial spines.Diagnosed from S. madhousensis by the longer length of LII (.13 mm).
FEMALES: (N = 11).Body length 2.8-3.68,scute length 2.3-2.69,greatest width of anterior scute 1.41-1.69,greatest width of node (FIGURE 5).Despite this, these two species are nonmonophyletic in all gene trees (FILE S6) and gsi statistics do not support monophyly for any genes (TABLE 2).Multispecies coalescent models, like that implemented in BPP, may still support species status despite the non-monophyly of gene trees, as incongruence among gene trees is attributed to deep coalescence [62,63].Previous analyses based on mitochondrial COI support the reciprocal monophyly of these two species [22] providing additional evidence for distinct species and the deep coalescence of nuclear genes.Biogeography also supports species status as these species occupy regions separated by a well-known biogeographical break seen in numerous species [64].This situation requires further study with denser geographic and genetic sampling, especially to determine if more fine-scale geographical patterns exist within either S. nondimorphicus or S. idahoensis.
Sclerobunus robustus group Diagnosis.The dorsal surface of the ventral plate of the penis has many folds, with a pair of small ventral spines on lateral plates (FIGURES 7c and 8).
Diagnosis.Diagnosed from all other surface species in the S. robustus group by its larger body size (scute length generally .2.2 mm) and longer legs (LII generally .7 mm).Some populations of S. robustus are very similar to S. skywalkeri in morphological characteristics (i.e., Bradford Canyon), but can be distinguished based on geographical distribution.Differentiated from cave-adapted species of the S. robustus group by its lack of troglomorphic features.Several populations of cave-inhabiting S. robustus are known, however they can be diagnosed from all species of the S. robustus group based on LII length (,8 mm in S. robustus, .8mm in troglomorphic species).
Material Distribution.Known from the southern Sangre de Cristo Mountain Range in New Mexico.
Comments.In prior analyses [21] the S. glorietus complex only included 4 samples with only one highly supported branch, a sister relationship between S. klomax and the S. glorietus population from Taos Ski Valley.In previous COI-only analyses [22], S. skywalkeri is monophyletic and highly supported, S. klomax and the S. glorietus from Taos Ski Valley were reciprocally monophyletic with low support, and S. glorietus from Glorieta Canyon was recovered sister to the entire S. robustus species group.Similar difficulties occurred in the present study in which morphological clusters were not supported by gene trees.In most gene trees, the Glorieta Canyon population of S. glorietus was recovered either within or sister to S. skywalkeri populations, however, support for these relationships was low (FILE S6).Despite this uncertainty, and given that BPP analyses are sensitive to guide tree topology [8], we argue that defining three species within the S. glorietus complex based on the morphometric clustering and qualitative morphological similarity is the most conservative and intuitive choice and unlikely to result in future synonymy with increased sampling.Further sampling focused on the geographic gap between the two populations of S. glorietus may eliminate the discrepancy between the morphological and genetic data seen here.
Etymology.The specific epithet is a noun in apposition named in honor of the Jemez Pueblo for whom Terrero Cave is culturally and spiritually significant.
Diagnosis.Diagnosed from all other species in the S. robustus group by a combination of scute length of 1.9-2.0mm and LII length between 8.0-9.0 mm. S. jemez is generally less troglomorphic than S. steinmanni, S. speoventus and S. klomax.The first three proximal ventral spines of the PF are stouter and shorter than other troglomorphic species.S. jemez has more spines on CI (8-12) compared to S. klomax (6)(7).
Etymology.The specific epithet is masculine and is a combination of the Latin words for cave (speo-) and wind (vent-), and refers to the type locality of Cave of the Winds.
Diagnosis.Differentiated from all surface Sclerobunus by a highly troglomorphic appearance.Distinguished from S. steinmanni by a relative lack of black pigment on the dorsum, slightly smaller size, the presence of a single ventral tubercle bearing spine on LI tibia and strongly curved subapical spines on the penis. Description

Species Delimitation
Species delimitation is becoming increasingly objective and taxonomy more integrative, utilizing both discovery and validation approaches employing multiple lines of evidence [3].Here, multiple lines of evidence are used to support all species of Sclerobunus: morphometric clustering, genetic analyses utilizing multiple nuclear loci, previously conducted mitochondrial analyses [21,22], and general morphology.To date, only two studies have used the method of Ezard et al. [7], the original study and a second by Pearson and Ezard [70], both examining morphological differentiation in fossil foraminifera.We implemented this method on extant taxa as a means to develop a priori species hypotheses to be tested using hypothesis-based genetic analyses.Here, the method was used on a dataset consisting of multiple species where it successfully delimited clusters representing multiple putative species.Any clusters containing more than one species were easily differentiated using phylogenetic analyses, as either the species were not sister taxa or validation analyses confirmed genetic distinctiveness.This dataset represents a relatively simple situation in which to apply this method considering that most species of Sclerobunus can be distinguished by eye.However, we believe this method has great utility in modern species delimitation studies and testing its efficacy is essential, especially with regards to species complexes with little obvious morphological differentiation.
Other protocols have been developed that rely on morphological data to identify morphometric clusters that are later used as a priori species hypotheses.For example, Seifert et al. [71] developed a hypothesis-free protocol that uses several morpho-metric clustering methods, which may be used in combination with other lines of evidence, to develop species hypotheses.The ''nest-centroid'' clustering method was specifically developed for and demonstrated with eusocial organisms where, in the case of ants, nest samples are treated as distinct classes in multivariate statistical analyses.The application and utility of this protocol using only morphometric data was shown by the identification and description of a new cryptic ant species [72].Given the increasing number of cryptic species being described [4], using the recently developed morphological-based methods [7,71] as a tool for identifying putative species, particularly when used in conjunction with genetic data, is promising.
The integrative nature of this study provides support for the utility of coalescent-based methods in species delimitation, which are becoming increasingly relevant given the common discovery of cryptic species complexes [4].Although mitochondrial genes have been used with great utility in animal systematics, many concerns have been raised about relying on mtDNA as the sole marker in phylogenetic and species delimitation analyses [73].Similarly, some studies have demonstrated that species trees derived from multilocus nuclear datasets can provide higher levels of resolution and support relative to analyses based solely on mitochondrial genes (e.g., [74]), empirically demonstrating the need for less reliance on mitochondrial data.Under an integrative taxonomic framework, mitochondrial and nuclear datasets represent independent lines of evidence [3] and should be analyzed separately.We did not reanalyze or include previously published COI data for Sclerobunus for several reasons.Thorough phylogenetic analyses of COI have already been conducted, including divergence dating and gsi analyses, documenting deep genetic divergences between species and populations [22].On average, animal mitochondrial genes sort four times faster than nuclear genes and tend to show more extreme levels of genetic divergence, particularly in low vagility, microhabitat specialist taxa.As such, nuclear genes provide a more conservative estimate for species delimitation decisions.We believe that for our system and dataset, including mitochondrial genes in multilocus nuclear validation analyses (BPP) would inflate genetic divergences between putative species and potentially result in a more liberal number of species delimited.For example, 37 of 41 populations sampled showed COI gsi values of 1 (monophyly), indicative of the extreme levels of mitochondrial genetic divergence seen in this genus [22]; however, at the species level, most nuclear genes still show some degree of non-monophyly (TABLE 2).

Delimiting Cave-obligate Species
The importance of genetic data in the systematics of caveadapted species is well documented [75,76]; however, studies using explicit species delimitation methodology in cave taxa are scarce (e.g., [12]).All cave species in this study do not experience gene flow with other cave species given the vast distances between caves.The cave species that are recovered as sister taxa in these analyses (the S. cavicolens group and S. steinmanni/S.speoventus) are not likely the result of speciation after an initial cave-invasion by an ancestral species.These cave species are from distinct, unconnected cave systems and suitable habitat is not present in the intervening regions between caves.Similarly, and in addition to morphological differences, we believe this same reasoning allows us to elevate the two subspecies of S. ungulatus to separate species.It is probable that these cave populations are each independently derived from different, now-extinct surface populations.Even with the discovery and increasing recognition of the importance of noncave subterranean habitats (superficial subterranean habitats) [77], post-invasion speciation is still unlikely given the vast distances between caves and likely impossible as the desert basin habitats involved do not possess superficial subterranean habitats.Gene flow with surface forms is possible in some cases, but most focal caves are surrounded by unsuitable surface habitats.

Cave Conservation and Exploration
This study highlights the need for further conservation efforts and exploration in western North American caves.For example, S. klomax and S. steinmanni are only known from a few specimens despite multiple collecting attempts.Most other cave species of Sclerobunus, and many North American arachnids in general, are rare in collections.However, despite the general rarity and numerous destructive threats facing caves and the organisms inhabiting them, only a handful of taxa are listed as endangered [78].Notably, all but one of the 12 arachnid species listed as endangered in the United States is a cave-obligate species (http:// ecos.fws.gov/tess_public/SpeciesReport.do?groups = J&listingType = L&mapstatus = 1).
Given the large number of caves in western North America not yet explored from a biological perspective, the potential for discovery of new cave populations of Sclerobunus is high.Records exist of cave-dwelling Sclerobunus robustus from several caves in isolated mountain ranges in Arizona [79], but these likely represent recent, local cave invasions within S. robustus.More importantly, several unsampled records represent potentially new species, in particular, a juvenile Sclerobunus collected from Crystal Falls Cave, Idaho [80].Additional research would help clarify the interesting situation uncovered at Taos Ski Valley.To date, only 3 female specimens of S. klomax have been found from rock piles within a very small area (,100 m 2 ) of south-facing forested-talus slope.Further fieldwork is needed to determine if this species has a broader geographic distribution and whether it only exists in superficial subterranean habitats or can be found in any nearby caves.Taken together, the description of four new, highly troglomorphic Sclerobunus from Colorado and northern New Mexico, and the possibility of other new cave-adapted species and populations, portray a picture of active troglomorphic evolution in a region that is not typically considered a hotspot for cave biodiversity.

Figure 1 .
Figure 1.Species delimitation decision-making workflow.Symbols to the left of each step represent the different classes of data used at that particular step: star = morphological data; circle = mitochondrial data (COI); square = nuclear data.doi:10.1371/journal.pone.0104982.g001

Figure 5 .
Figure 5. *BEAST species tree with BPP species delimitations.Numbers below nodes on species tree correspond to posterior probability support.Earliest diverging braches shortened for graphical purposes.At bottom left, BPP results with 4 different prior combinations.Numbers correspond to results of the two separate BPP runs.Red letters correspond to those nodes on the species tree.doi:10.1371/journal.pone.0104982.g005

Table 2 .
Results of gsi analyses.
Penis with dorsal surface of ventral plate with many folds (FIGURES 7c and 8); southwestern North America (AZ, NM, UT, CO)… S. robustus group Sclerobunus cavicolens group Diagnosis.Penis with lateral plates extending dorsally into fan-shaped projections and extending ventrally into bifurcate Scute length 1.7-2 mm; LII length under 12 mm; caves near Provo, UT… S. madhousensis comb.nov., stat.nov.Scute length over 2 mm; LII length over 13 mm; caves in Great Basin National Park, NV… S. ungulatus comb.nov.