New species in the Sitalcina sura species group (Opiliones, Laniatores, Phalangodidae), with evidence for a biogeographic link between California desert canyons and Arizona sky islands

Abstract The western United States is home to numerous narrowly endemic harvestman taxa (Arachnida, Opiliones), including members of the genus Sitalcina Banks, 1911. Sitalcina is comprised of three species groups, including the monospecific Sitalcina californica and Sitalcina lobata groups, and the Sitalcina sura group with eight described species. All species in the Sitalcina sura group have very small geographic distributions, with group members distributed like disjunct “beads on a string” from Monterey south to southern California and southeast to the sky-island mountain ranges of southern Arizona. Here, molecular phylogenetic and species delimitation analyses were conducted for all described species in the Sitalcina sura group, plus several newly discovered populations. Species trees were reconstructed using multispecies coalescent methods implemented in *BEAST, and species delimitation was accomplished using Bayes Factor Delimitation (BFD). Based on quantitative species delimitation results supported by consideration of morphological characters, two new species (Sitalcina oasiensis sp. n., Sitalcina ubicki sp. n.) are described. We also provide a description of the previously unknown male of Sitalcina borregoensis Briggs, 1968. Molecular phylogenetic evidence strongly supports distinctive desert versus coastal clades, with desert canyon taxa from southern California more closely related to Arizona taxa than to geographically proximate California coastal taxa. We hypothesize that southern ancestry and plate tectonics have played a role in the diversification history of this animal lineage, similar to sclerophyllous plant taxa of the Madro-Tertiary Geoflora. Molecular clock analyses for the Sitalcina sura group are generally consistent with these hypotheses. We also propose that additional Sitalcina species await discovery in the desert canyons of southern California and northern Baja, and the mountains of northwestern mainland Mexico.


Introduction
Laniatorean harvestmen comprise the majority of Opiliones diversity, with more than 4100 described species (Kury 2013). North temperate Laniatores are typically small (body length usually less than 3 mm), short-legged animals, most often found in sheltered microhabitats such as under rocks and logs. The western United States is home to a diverse laniatorean fauna, with more genera and species than any other region in the Nearctic, making this area a harvestman "hotspot" (e.g., Briggs 1971, Ubick and Briggs 1989, Ubick and Briggs 2008, Derkarabetian and Hedin 2014. Laniatores of the family Phalangodidae are mostly Holarctic in distribution, but particularly diverse in the Nearctic (Ubick 2007). In California there are about 70 described phalangodid species, and this number has been growing with further sampling and use of SEM studies of morphology. Most of the described species have very small geographic distributions, and deserve more evolutionary, biogeographic, and conservation interest than currently afforded (Ubick andBriggs 2008, Emata and. The phalangodid genus Sitalcina is comprised of three species groups (Ubick and Briggs 2008), including the monospecific S. californica and S. lobata groups, and the S. sura group which includes eight described species: S. sura Briggs, 1968, S. seca Ubick & Briggs, 2008, S. chalona (Briggs, 1968), S. flava Briggs, 1968, S. borregoensis, S. rothi Ubick & Briggs, 2008, S. catalina Ubick & Briggs, 2008and S. peacheyi Ubick & Briggs, 2008. With a known distribution from Monterey south to southern California and southeast to the sky-island mountain ranges of southern Arizona (Figure 1), species of the S. sura group occupy a variety of upland habitats. Most of these habitats are dominated by sclerophyllous woody plant taxa (e.g., Quercus, Pinyon pine, Arctostaphylos, Ceanothus, etc.), part of the Madro-Tertiary Geoflora (MTG, Raven and Axelrod 1978, Lancaster and Kay 2013, Baldwin 2014. Exceptions include S. sura from redwood forests, and S. borregoensis from desert canyons, below the elevational level of desert chaparral. Species from the sky-island mountain ranges of southern Arizona occur in mid-elevation Madro-Tertiary habitats, largely below higher elevation conifer forest, but above low desert habitats. In all habitats S. sura group members occupy seemingly similar microhabitats, typically under rocks on shaded north-facing slopes (Ubick and Briggs 2008, Figure 2).
Because of apparent limited dispersal abilities and microhabitat specificity, extreme population genetic structuring and divergence can be expected in Sitalcina, as is observed in other low vagility Laniatores (e.g., Thomas and Hedin 2008, Hedin and Thomas 2010, Derkarabetian et al. 2011, Emata and Hedin 2016. As such, some currently described species with naturally fragmented distributions (e.g., S. peacheyi on isolated mountaintops and in caves) may actually comprise lobata (1 species), and S. sura (8 described species, 4 geographically novel populations). Sampled populations indicated by circles. General species distributions follow Ubick and Briggs (2008). multiple cryptic species. Also, the geographic sampling of previous studies may have missed unique species. Finally, members of the S. sura group are very similar in somatic morphology, likely due to niche conservatism (Wiens andGraham 2005, Keith and. Overall, the S. sura group has the potential to include both newly discovered and cryptic species.
A variety of objective species delimitation methods have been developed in recent years (e.g., Yang and Rannala 2010, Rannala and Yang 2013, Grummer et al. 2014. Several of these species delimitation methods are founded on the multispecies coalescent model, utilizing multilocus genetic data (summarized in Fujita et al. 2012, Carstens et al. 2013. The probable existence of morphologically cryptic species in Sitalcina makes the use of such genetic methods attractive. Also, because Sitalcina populations are allopatric and mostly separated by unsuitable habitat, it can be assumed that interspecific gene flow is minimal. As such, incongruence in gene tree topologies should largely reflect incomplete lineage sorting, consistent with the assumptions of most multispecies coalescent methods (Fujita et al. 2012). Conversely, extreme population genetic structuring across naturally fragmented habitats may represent a deviation from model assumptions (Niemiller et al. 2012, Hedin et al. 2015.
In this research we first use species discovery approaches to formulate alternative species delimitation hypotheses for the S. sura group. Based on genetic species delimitation results, supported by consideration of morphology, two new species and the previously unknown male of S. borregoensis are described. A time-calibrated multilocus species tree is used as a framework to interpret the biogeographic history of the S. sura group. We hypothesize that this history is linked to both plate tectonics and southern ancestry, similar to elements of the MTG.

Taxon sampling
Fieldwork was conducted in the winter and spring months when surface microhabitats were most suitable for successful collections. Voucher specimens used for morphological study and species descriptions were preserved in 80% ETOH, while those used in genetic analyses were preserved in 100% ETOH at -80 °C. All described species of the S. sura group were collected from at or near type localities (Ubick and Briggs 2008). This sampling included previously discussed (Ubick and Briggs 2008), but undescribed, specimens from the Santa Ynez Mountains, California. In addition, several new populations from California and Arizona were discovered in appropriate habitats. Two specimens per locality were used for phylogenetic analyses when available, with the final sample including 29 ingroup specimens from 16 localities (Figure 1, Suppl. material 1: Table S1). Locality data for all specimens are also available on the Symbiota Collections of Arthropods Network (http://symbiota4.acis.ufl.edu/scan/portal/index. php). Specimens are housed in the San Diego State Terrestrial Arthropod Collection (with SDSU_TAC or SDSU_OP catalog numbers); type specimens are deposited at the California Academy of Sciences (CAS). Sitalcina californica (Banks, 1893) and S. lobata Goodnight & Goodnight, 1942 were used as outgroup taxa in all analyses. Here we accept the hypothesis that Sitalcina is monophyletic, as argued on morphological grounds by Ubick and Briggs (2008).

DNA isolation, amplification, and sequencing
Genomic DNA was extracted from leg tissue (2-3 legs) using the Qiagen DNeasy Kit (Qiagen, Valencia, CA). DNA fragments for mitochondrial cytochrome oxidase I (COI) and nuclear 28S rRNA, plus five additional protein-coding nuclear genes, were amplified using PCR (Table 1). Primers and cycling conditions are reported in Suppl. material 1: Table S2. Primers for protein-coding nuclear genes targeted terminal exon plus associated 3' untranslated regions (3'-UTRs), and were developed by comparing transcriptomes of S. lobata  to the phalangodid taxon Texella bifurcata (Briggs, 1968) (SRA numbers reported in Emata and Hedin 2016). Both transcriptomes were generated using Illumina short-read technology, and assembled using Trinity software (Grabherr et al. 2011). PCR products were purified using Millipore plates, Sanger sequenced in both directions at Macrogen USA (Rockville, MD), and edited using Geneious Pro 6 ( Kearse et al. 2012). Sequences were aligned using MUSCLE (Edgar 2004), with the exception of 28S data, which were aligned using MAFFT (Katoh and Standley 2013) and the G_INS-I alignment algorithm (Wilm 2006). Alleles from heterozygous nuclear sequences were inferred using PHASE 2.1.1 (Stephens et al. 2001, Stephens andScheet 2005), with each analysis repeated twice to ensure consistent results.

Gene trees and genetic clustering
Models of DNA sequence evolution were chosen using jModeltest2 (Guindon andGascuel 2003, Darriba et al. 2012), with the Akaike Information Criterion (AIC) used to select models (Table 1). Mitochondrial COI gene tree analyses were conducted using both codon partitioned and un-partitioned models. Gene trees for individual loci were reconstructed using MrBayes 3.2 (Ronquist et al. 2011) run on the Cyber Infrastructure for Phylogenetic Research (CIPRES; Miller et al. 2010). Bayesian MCMC analyses for each gene were run for 50 million generations sampling every 5000 generations. Each analysis was repeated twice to confirm results. ESS values along with -lnL scores were evaluated for convergence using Tracer v.1.6 (Rambaut et al. 2014). A 50% majority rule consensus tree from the resulting posterior distribution was constructed for each gene region. We used genetic clustering in the initial "discovery" phase of species delimitation, a procedure common in the recent literature (reviewed in Carstens et al. 2013). STRUCTURE 2.3.4 (Pritchard et al. 2000) was used to analyze biallelic nuclear genotypic data. Each unique haplotype was treated as a single allele and was called using SNAP-map (Aylor et al. 2006). STRUCTURE runs were conducted with 20 iterations for multiple K values (K = number of genetic clusters). Separate runs were conducted with coastal clade species (see Results) with a K of 1-9, and desert clade species (see Results) with a K of 1-8. Analyses were run for 100,000 generations, with the first 10,000 generations removed as burnin. Analyses were run using both a noadmixture model (assumes each individual comes from one of the K distinct populations) and an admixture model (allows for population admixture). All other priors were left as default. Structure Harvester (Earl 2012) was used to find the best-fit K utilizing the ∆K method of Evanno et al. (2005). Data were summarized using the FullSearch algorithm of CLUMPP (Jakobsson and Rosenberg 2007), and visualized with DISTRUCT (Rosenberg 2004).

Species trees and divergence times
Multilocus species trees were reconstructed using *BEAST 1.8 . A priori species limits were based on genetic clusters identified by STRUC-TURE admixture and no-admixture model results. For each analysis, substitution models, clock models, and trees were unlinked across loci. A best-fit model of molecular evolution was applied to each gene region (Table 1), with mitochondrial COI partitioned by codon and a Yule process applied for the species tree prior. Clock rate parameters were examined using Tracer v1.6 (Rambaut et al. 2014). If the 95% highest posterior density (HPD) of the coefficient of variation for any individual gene included zero (indicating that a strict molecular clock could not be statistically rejected), a strict clock was used for that gene region. Subsequent analyses were run with the appropriate uncorrelated relaxed or strict clock models (Table 1). Analyses were run for 200 million generations logging every 2000 generations. ESS values along with -lnL scores were evaluated for convergence using Tracer. Each species tree analysis was repeated twice to ensure consistent results. Posterior distributions of repeated runs were combined using LogCombiner, and a maximum clade credibility (MCC) tree was constructed from the resulting combined posterior distributions using TreeAnnotator.
Divergence time analyses were performed in *BEAST ). Because of general uncertainty in rates of molecular evolution in Laniatorean harvestmen, we conducted three separate analyses to provide bounds on possible dates. First, a well-accepted arthropod COI clock rate of 3.54% per million years (Ma) (arithmetic mean of branch rates (ucld) setting in *BEAST = 0.0178) was specified for the partitioned COI data (Papadopoulou et al. 2010). Second, the unpartitioned COI rate of Papadopoulou et al. (2010) was used (ucld mean = 0.0169). Finally, a slower COI rate calculated for the laniatorean genus Sclerobunus (ucld mean = 0.01115) was specified for the COI data partition (Derkarabetian et al. 2010). This latter rate was inferred indirectly based on combined biogeographic and fossil calibrations. For each analysis, a Yule process was specified for the species tree prior and a pairwise linear and constant root was applied for the population model, with a uniform distribution on the upper and lower bounds of the age. Substitution models, clock models, and trees were unlinked across loci. The appropriate model of molecular evolution and clock-like rate was applied to each gene fragment (see above). Analyses were run for 50 million generations sampling every 500, and repeated to confirm consistent results. Repeated runs were combined using LogCombiner, and a MCC tree was constructed using TreeAnnotator.

Genetic species delimitation
Bayes factor delimitation (BFD, Grummer et al. 2014) allows testing of alternative species delimitation models by assigning individuals to different lineages, with subsequent comparisons of marginal likelihoods of alternative models. We tested the alternative species delimitation models summarized in Table 2 -these alternative models are based on a combination of STRUCTURE results, geographic considerations, and prior taxonomy. Marginal likelihood estimations were run for 100,000 generations, sampled every 1000 with 100 steps using path sampling (Lartillot and Philippe 2006) and stepping stone (Xie et al. 2011) methods. A comparison of marginal likelihoods was conducted using Bayes factors, with values above 10 considered as decisive support (Kass and Raftery 1995).

Niche modeling
With incomplete knowledge of the full distribution of the S. sura group, possible regions with additional new populations or species were identified via ecological niche modeling (ENM). DIVA-GIS (Hijmans et al. 2004) and Maxent (Phillips and Dukid 2008) were used to construct niche models, using known localities for the species group as input. Niche models have previously successfully predicted the potential habitats of animal taxa with small and fragmented distributions (e.g., Rissler and Apodaca 2007, Bond and Stockman 2008, Bond 2012. Predicted distributions of the S. sura group were reconstructed using altitudinal and climatic layers for nineteen quarterly and annual measurements of temperature and precipitation (Bioclimatic layers 1-19), obtained from the WorldClim dataset (Hijmans et al. 2005). A jackknife analysis was conducted in Maxent to discover the most likely environmental factors impacting species' distributions. Output from Maxent was converted to raster format and reclassified to binary presence/absence in ArcMap 10.3 (ESRI) using the 10 percentile training presence logistic threshold.

Study of morphology
Male penises that were not protruding from the genital operculum were physically extracted using a blunt insect micro pin. Exposed penises were placed in room temperature (or hot) 10% KOH for 1-2 minutes for expansion. Female ovipositors were Putative species grouped with geographic neighbors exposed using the same blunt pin procedure. Specimens were imaged using a Quanta 450 scanning electron microscope (SEM) after being mounted and coated with 20nm platinum. One or two specimens were used for SEM as needed. Whole specimen digital images were captured using a Visionary Digital BK plus system (http://www.visionarydigital.com). Individual images were merged into a composite image using Helicon Focus 6.2.2 software (http://www.heliconsoft.com/heliconfocus.html). Specimen measurements were taken with an SZX12 Olympus dissecting scope equipped with an ocular micrometer, at 50× magnification.

Data characteristics
DNA sequence data were collected for two outgroup taxa (data for S. lobata from transcriptomes), eight described species from the S. sura group, and four novel geo-graphic populations with morphological features placing them in the S. sura group (see diagnostic features below). Not all specimens were successfully amplified for all gene regions, resulting in some missing data (Table 1). All heterozygous nuclear sequences were PHASED with 100% certainty. GenBank accession numbers for unphased sequences are reported in Suppl. material 1: Table S3, and phased data matrices have been submitted to the Dryad Digital Repository (http://dx.doi.org/10.5061/dryad.4gk4f ).

Gene trees and genetic clustering
Mitochondrial COI gene trees were reconstructed using both a codon partitioned and un-partitioned model, and are generally topologically similar (Figure 3). A clade that includes southern California desert populations together with Arizona sky-island populations (hereafter named the "desert" clade) is recovered in both analyses. Geographically proximate populations from canyons in the Anza Borrego desert (Borrego Palm Canyon, S. borregoensis) are surprisingly genetically divergent, with average Kimura 2-parameter distances exceeding 17% (calculated using MEGA6, Tamura et al. 2013), and do not form a clade on COI trees. An important incongruence between codon partitioned versus un-partitioned results concerns the uncertain placement of the S. flava + Palos Verdes clade. These populations are consistently placed with other coastal CA populations (hereafter named the "coastal" clade, see below) in nuclear gene tree analyses, but group with the desert clade with low support (posterior probability (PP) = 0.59) in codon partitioned COI analyses. Despite expected variance in nuclear gene tree topologies, several general trends are apparent. First, a coastal CA clade is recovered (generally with strong support, PP > 0.95) in all six nuclear gene trees (Figure 4). This clade includes S. sura, S. seca, S. chalona, S. flava, Santa Ynez, and Palos Verdes populations. Specimens from Santa Ynez are genetically distinct in all nuclear gene trees. A desert clade is strongly supported in four of six nuclear gene trees; paraphyly of this group in two gene trees may reflect a gene tree rooting issue. Within the desert clade, Borrego Palm Canyon specimens are genetically distinct from neighboring S. borregoensis in four gene trees where both populations were sampled (Figure 4). STRUCTURE analyses for coastal clade specimens favor a K = 6 model (Figure 5A). Both no-admixture and admixture models place S. sura, S. seca, S. chalona, and Santa Ynez individuals into four separate genetic clusters. Specimens from Palos Verdes and S. flava from Piuma Road group together into a fifth genetic cluster, while Topanga Canyon S. flava specimens represent a sixth cluster. The STRUCTURE admixture model favors all discrete desert populations (either found on isolated mountain ranges or in isolated desert canyons) as unique genetic clusters (K = 7, Figure 5B). A no-admixture model favors K = 6, grouping geographically adjacent Glen Oaks and S. rothi specimens ( Figure 5C).

Species trees and divergence times
*BEAST analyses conducted with admixture STRUCTURE clusters as a priori species recover desert and coastal clades with strong support (PP > 0.95, Figure 6A). Relationships among coastal taxa are well supported, with an internal topology that mirrors geography (early-diverging southern lineages, Santa Ynez intermediate, derived northern lineages). Relationships between desert taxa are less strongly supported, with CA desert canyon populations forming a clade sister to montane Arizona populations.
Using the Papadopoulou et al. (2010) partitioned COI rate, the estimated divergence time for the most recent common ancestor (tMRCA) of the S. sura group is ~18 Ma, the tMRCA for the coastal clade is ~10 Ma, and the tMRCA for the desert clade is ~7 Ma ( Figure 6B). The Papadopoulou et al. (2010) unpartitioned COI rate provides very similar time estimates, while the slower Sclerobunus rate results in clearly older time estimates for these same nodes (Table 3). In general, we have no a priori reason to favor one of these rates, given the general lack of knowledge of rates of molecular evolution in laniatorean harvestmen. However, we favor the younger dates for three reasons. First, the Papadopoulou et al. (2010) rates result in date estimates that coincide with important biogeographic events in the region (see Discussion). Second, various temporal analyses conducted by Emata and Hedin (2016) for the California phalangodid genus Calicina also showed that the Papadopoulou et al. (2010) rate provided biologically realistic dates for a California taxon, while slower rates suggested unrealistically old divergence dates. Finally, we note that recent studies of Sclerobunus using whole genome SNP data ) provide younger divergence ages than previously hypothesized in that system, suggesting that the Derkarabetian et al. (2010) COI rate may have been underestimated. All of these arguments are ad hoc, with resolution in the Sitalcina system ultimately requiring additional data.

Niche modeling
Precipitation in the coldest quarter (BIO19) and precipitation in the driest month (BIO14) were the best predictors of S. sura group distributions. ENMs based on these variables provide a visualization of possible sampling gaps within the S. sura group (Figure 7). In California and Baja California Norte, areas of high suitability that represent possible gaps include the Santa Ynez Mountains northwest of Santa Barbara, and the southern Santa Lucia Mountains. Some of these habitats are occupied by S. californica, which has an apparently exclusive distribution with members of the S. sura group (i.e., no sympatry has ever been recorded, Figure 1). Similarly, we expect S. lobata to occupy the predicted suitable habitats of coastal southern CA and northern Baja (Figure 1). In these particular cases, we hypothesize that the ENM has over-predicted the distribution of the S. sura group. California desert canyon habitats, both north and south of our sampled populations, have predicted suitability. In Arizona and northern mainland Mexico (Sonora), unsampled populations are predicted to occur in montane habitats between known S. catalina -S. rothi populations, and south of our current sample. We are not aware of Sitalcina records from northern mainland Mexico, but

Genetic species delimitation
BFD results support a 13 species model, following the K = 6 STRUCTURE model for coastal taxa and the K = 7 model for desert taxa (Tables 2 & 4, Figure 6). This 13 spe-cies model supports three new populations as putative species (Borrego Palm Canyon, Santa Ynez, Glen Oaks), and the division of both S. peacheyi and S. flava into separate species. This model is decisively supported over alternative species delimitation models (Table 4). In translating BFD results into formal taxonomy, we have taken a conservative approach, also considering support from other lines of evidence. Because Glen Oaks is known only from a single adult female specimen, we defer formal species description to a later date, when additional specimens can be collected and studied. Further attention should be also directed at the geographically disjunct Palos Verdes population, which is genetically distinct from S. flava for multiple genes (Figures 3, 4). Finally, because we failed to find obvious morphological differences between disjunct populations of S. peacheyi and S. flava, and because multispecies coalescent models might oversplit genetically structured populations (Niemiller et al. 2012, Hedin et al. 2015, we do not further split S. peacheyi and S. flava at this time. Genetically distinct and morphologically diagnosable populations from Borrego Palm Canyon and Santa Ynez are formally described below. Also, previously unknown males of S. borregoensis are described. Following Ubick and Briggs (2008), we emphasize both somatic and genitalic characters in our descriptions.

Sitalcina sura Group
Diagnosis. Members of the S. sura group are distinguished from related S. californica and S. lobata by the following characters (Ubick and Briggs 2008): Both sexes possess a short row of dorsomesal asetose tubercles on the palpal Fm. Males lack an AS, possess a bilobed PSL, and a curved to straight ectal spur on TrIV. Females with imbricate OVM (absent/ weak in the (S. chalona, (S. sura, S. seca)) clade), and curved OVS with brush-like tips. Description. FEMALE. As in Ubick and Briggs 2008. MALE. Integument color pale orange, appendages lighter. Body finely rugose with a few large tubercles on posterior tergites, one pair anteriorly on EM; 3 pairs of AT. EM low, flattened, eyes present. Palpal Fm with median dorsobasal row of 4 asetose tubercles and one small mesal tubercle. Palpal megaspines: trochanter one ventral and small; Fm 3 ventrobasal, one mesodistal; patella 2 mesal, one ectal; tibia and tarsus 2 mesal, 2 ectal. TC 3-5-5-5.
Distribution and habitat. Known only from the vicinity of Mountain Palm Springs, Anza Borrego Desert State Park. New collections are from a north-facing slope, under the first layer of granite rocks in a small ravine adjacent to a palm grove (Figure 2A). Description. Integument color pale orange with lighter appendages. Body finely rugose with larger tubercles along tergal margins and small tubercles anteriorly on EM; 2-3 pairs of AT. EM low and rounded, eyes present. Palpal Fm with median dorsobasal row of 3 asetose tubercles and one small mesal tubercle. Palpal megaspines: trochanter 2 ventral; Fm 3 ventrobasal, one mesodistal; patella 2 mesal, one ectal; tibia and tarsus 2 mesal, 2 ectal. TC 3-5-5-5.
Other material examined. See Suppl. material 1: Table S1 for locality information for specimens examined.
Distribution and habitat. Known only from the type locality. Specimens were collected from sparse desert chaparral habitat, under small rocks amongst larger granite boulders, NE-facing slope above palm oases.  Etymology. This species is named in honor of Darrell Ubick (CAS) whose foundational taxonomic research with Sitalcina made this project possible.
Distribution and habitat. Known only from the vicinity of Montecito, Santa Ynez Mountains, Santa Barbara County. Specimens were found under stones, in narrow ravine with stream, in a Quercus and Platanus forest.
Note. This population was mentioned by Ubick and Briggs (2008), but specimens were damaged and left undescribed. The new material collected for this study is from the same general vicinity as specimens studied by Ubick and Briggs (2008); northern edge of Montecito, off East Mountain Drive.

Sitalcina remains poorly known
Species limits in the S. sura group were first studied by Ubick and Briggs (2008) based on consideration of somatic and genital morphology. In the current study, these and other species hypotheses (Table 2) were formally tested using multilocus DNA sequence data and a Bayes Factor Delimitation approach. With two newly described species, Sitalcina now includes a dozen species, and we expect future research to uncover additional undescribed taxa. Ultimately, we predict that Sitalcina diversity might approach that observed in the phalangodid genus Calicina, a taxon that includes over 25 short-range endemic species from central and northern California (Ubick andBriggs 1989, Emata and. Our results imply that some described species (e.g., S. peacheyi, S. flava) may actually comprise multiple cryptic species. This pattern is expected in a morphologically conservative group that occupies naturally fragmented habitats. This prediction would also apply to the geographically widespread S. californica (see Ubick and Briggs 2008), which has an interesting disjunct distribution in California (Figure 1). We have also discovered new geographic populations (e.g., Glen Oaks, Palos Verdes) that may represent new species. For both new populations and potential cryptic species, future research implementing genomic-scale datasets would be informative. Finally, our results confirm highly localized geographic speciation in these sedentary animals. Examples include geographically adjacent S. seca, S. sura, and S. chalona, as well as S. borregoensis plus S. oasiensis. If we extrapolate this pattern of localized divergence to regions where Sitalcina is likely, but currently unknown, this suggests many additional species. If Sitalcina is distributed like "beads on a string", we are missing many of the beads.
Sitalcina specimens are sometimes very challenging to collect. Specimens only occur near the surface (i.e., under accessible rocks or woody debris) at certain times of the year, apparently migrating deep into the soil matrix as surface conditions dry (Ubick and Briggs 2008). This means that collections must take place in "good" years (e.g., non-drought years), in appropriate months, by experienced collectors in suitable microhabitats. An example is the Mountain Palm Springs population of S. borregoensis, which was first collected in 1967 (Ubick and Briggs 2008), but not recollected until 2012. Glen Oaks is another example, with only a single specimen collected despite over 3 hours of search time by a team of six experienced collectors. In this context, new populations of Sitalcina in poorly collected regions are likely, and we highlight the montane sky islands of northern Mexico and the desert canyons of southern California plus Baja California Norte as potential undiscovered diversification hotspots.
If multiple known species occur on the handful of sky-islands in southeastern Arizona, we expect additional undescribed species on the dozens of montane sky islands in northern Sonora. For many lineages, montane populations in southern Arizona constitute a "northern tip of a southern iceberg", with centers of distribution found in northern Mexico (e.g., montane jumping spiders - Maddison andMcMahon 2000, montane rattlesnakes -Bryson et al. 2011). The desert canyons of southern California plus Baja California Norte constitute another distributional area with potentially high unknown diversity. This region is tectonically active. For example, the population at Mountain Palm Springs lays almost directly on the Elsinore Fault, while the Borrego Palm Canyon population inhabits the area between the San Felipe Fault and the more easterly San Andreas Fault (Steely et al. 2009). Furthermore, relictual canyon populations are expected to be strongly isolated by unsuitable xeric habitats, leading to active speciation. As an example, Bond (2012) described multiple species of cryptic trapdoor spider species from a small region in Anza-Borrego Desert State park, with each species known only from one or two locations. Similar desert canyon local endemics likely occur in Sitalcina.

Historical biogeography
Niche modeling reveals a conspicuous geographic gap between habitats that are suitable for Sitalcina, separating eastern habitats in Arizona from western habitats in California (Figure 7). This gap is not only consistent with biotic features that separate these regions (i.e., low elevation Sonoran desert), but also coincides with major modern day and historical biogeographic barriers. These biogeographic barriers include marine incursions into the Salton Trough associated with the opening of the Gulf of California (ca. 6.3 Ma, Oskin andStock 2003, Dolby et al. 2015), and the more recent drainage of the Colorado River into the Gulf (ca. 4.1 Ma, Dorsey et al 2007, Dolby et al. 2015. Moreover, populations on either side of this biogeographic gap are currently found on different continental plates -all S. sura group populations from California are found on the Pacific plate, while all populations from Arizona are found on the North American plate (we note here that Anza-Borrego canyon populations are close to plate boundaries, further discussed below). For these reasons, a logical prediction would include a primary phylogenetic division in this region, resulting in California versus Arizona clades. Major phylogenetic or phylogeographic splits coincident with the Colorado River (or Salton Trough) have been found in many animal taxa (e.g., Riddle et al 2000, Riddle and Hafner 2006, Crews and Hedin 2006, although we emphasize that these studies have focused mostly on low elevation desert species. We failed to recover this predicted pattern, but instead recovered a well-supported phylogenetic pattern in both gene and species trees that includes the separation of S. sura group members into coastal versus desert clades (Figures 3, 4, 6). Here, southern Californian populations from the Anza Borrego desert region are phylogenetically allied with Arizonan populations, rather than with geographically proximate Californian populations. Because all lines of evidence indicate that Sitalcina is dispersal-limited, we hypothesize that vicariance has dominated the biogeographic history of this lineage. Below we present a biogeographic model for Sitalcina that emphasizes the role of plate tectonics and habitat connections to elements of a MTG.
The fundamental spatial premise of our model is that the S. sura group has a center of origin in northwestern Mexico, perhaps centered around the region which is today the northern Gulf of California. The fundamental temporal premise is that initial diversification in the group happened before the Pacific plate began to migrate actively northwestwards (along the San Andreas Fault system) against a stationary North American plate (Figure 13). Geological evidence indicates older rifting beginning as early as ca. 30 Ma in northwestern Mexico, with more active northwestern plate movement sometime after ca. 12.3 Ma (Atwater andStock 1998, Dolby et al. 2015). Furthermore, our model postulates initial diversification along a west to east axis, with westerly coastal clade members (plus S. lobata and S. californica), and easterly desert clade members. Importantly, this eastern lineage is hypothesized to have straddled both continental plates prior to active plate movements ( Figure 13). Our estimates for the overall age of the S. sura group, and for the tMRCA of desert versus coastal clades are generally consistent with this predominantly Miocene timeframe. Also, estimated times for the separation of CA canyon taxa from remaining desert clade members are consistent with the evolution of the Gulf of California (tMRCA of 7.09, Figure 6B, versus ca. 6.3 Ma, Oskin andStock 2003, Dolby et al. 2015).
Our model has interesting connections to the botanical literature. As noted in the Introduction, most members of the S. sura group (and Sitalcina more generally) are found in habitats dominated by sclerophyllous woody plant taxa (e.g., Quercus, Pinyon pine, Arctostaphylos, Ceanothus, etc.). These plant lineages are part of a MTG, as originally defined by Raven and Axelrod (1978). Although this concept was invoked prior to modern phylogenetic analyses, recent molecular studies have generally supported core aspects of this idea (Lancaster andKay 2013, Baldwin 2014). Important properties of plant lineages included in this flora include southwestern ancestry, deepest divergences occurring in the Oligocene and/or Miocene (e.g., Lancaster and Kay 2013, figure 2), a proposed role for plate tectonics (e.g., Raven and Axelrod 1978, fig. 4), and biogeographic connections between California and Arizona upland taxa (Raven and Axelrod 1978, table 7). These aspects of the MTG are consistent with our biogeographic model. Another interesting parallel is that many Madro-Tertiary taxa now found in coastal southern California are hypothesized phylogenetic relicts, found only in suitable habitats directly adjacent to a marine influence (Raven and Axelrod 1978). In this regard, we note that coastal clade S. sura group species are all found in moist mountains with direct marine influence (e.g., Palos Verdes peninsula, Santa Monica Mountains, Santa Ynez Mountains, Santa Lucia Range), with the single exception of S. chalona.
Several predictions can be derived from our biogeographic model. First, we expect undescribed species in upland habitats of western Sonora, some of which may be related to California desert canyon taxa. Second, we predict undiscovered populations in desert canyons of northern Baja, and again predict possible mixed phylogenetic affinities (some coastal clade, some desert clade). Third, we predict that northern populations of S. californica on the North American plate (north of San Francisco) should be phylogenetically derived from southern populations on the Pacific plate ( Figure 1). Finally, if our model has generality, we predict that future studies of other regional dispersal-limited animal lineages will show similar biogeographic patterns.
Both California and deserts of the American southwest are active areas for modern biogeographic research in animals. Despite this fact, we have found relatively few animal taxa that show biogeographic patterns similar to those observed in Sitalcina. Salamanders in the Batrachoseps pacificus group have a distribution much like the coastal clade (from Monterey south to northern Baja, almost all populations restricted to Pacific plate, Jockusch et al. 2014), and the group includes isolated southern California desert populations, but these populations have apparent western ancestry (Martínez-Solano et al. 2012,). Aptostichus trapdoor spiders include many California coastal, California desert canyon, and Arizona sky island species (Bond 2012), but phylogenetic details remain uncertain. One compelling parallel is found in night lizards (Xantusia). Leavitt et al. (2007) found that taxa from mostly transmontane southern California and northern Baja (including X. wigginsi, "San Jacinto" Clade, and "Yucca Valley" Clade) are sister to X. bezyi from uplands of central Arizona. Western species are found in pinyon-juniper and desert chaparral, while eastern X. bezyi occur in desert chaparral. These lizard taxa are hypothesized to have diverged 6.4 Ma, coincident with the evolution of the Gulf of California (Leavitt et al. 2007). All of these biogeographic patterns are similar to those observed in the Sitalcina desert clade.

Conclusions
We have uncovered an apparently novel phylogenetic pattern in a biogeographically wellstudied region. Our biogeographic model can be tested with additional research, both in Sitalcina, but also in other plant and animal lineages with similar geographic distributions. Ubick and Briggs (2008) suggested that Sitalcina likely included additional undescribed species -our data support this contention, and suggest even more undiscovered richness. Basically all members of the S. sura group have very small geographic distributions, and many represent apparently old lineages found in often disturbed habitats. Examples include the Palos Verdes population, known only from a single impacted canyon in a matrix of urban development. Similarly, the two known species from Anza-Borrego Desert State Park are each known only from single desert canyons, both located in close proximity to popular hiking trails. These short-range endemic taxa should receive more conservation attention than currently afforded (Harvey 2002, Harvey et al. 2011.

Supplementary Tables 1-3
Authors: Angela DiDomenico, Marshal Hedin Data type: Excel Table  Explanation note:  Table S1. Voucher and locality data Table S2. PCR primer and reaction information Table S3. GenBank accession information Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.