Alpine endemic spiders shed light on the origin and evolution of subterranean species

Sampling

Populations were collected in caves, abandoned mines and other hypogean habitats across the known distribution range of Troglohyphantes vignai and Pimoa rupicola in the Western Alps. The distribution range of T. vignai was entirely covered (including type locality and former localities of T. rupicapra, indicated hereinafter as T. vignai sensu rupicapra), while for Pimoa we only sampled Alpine populations, thus excluding French and Apenninic populations. The complete list of localities is reported in Table 1 and Fig. 1. The toponomastics and classification of the different sectors and sub-sectors of the Alps follows the partition of the Alpine chain (SOIUSA: Marazzi, 2005). Specimens were hand-collected, preserved in 95% ethanol and stored in freezer. Given that the sampled environments were highly oligotrophic, in certain localities we were able to detect and collect only few individuals of the two investigated species. The number of individuals collected for each locality ranged from 2 to 8 in P. rupicola and P. n. sp. and from 1 to 7 in T. vignai (Table 2). Overall, 119 Pimoa specimens from 25 localities and 37 Troglohyphantes specimens from 8 localities were used in this study.

DNA extraction, amplification and sequencing

One leg was removed from each specimen for DNA extraction. Whole genomic DNA was extracted from the samples using the SpeedTools Tissue Extraction Kit (Biotools) following the manufacturer’s protocol. A 676 bp region of the mitochondrial cytochrome oxidase subunit I (cox1) gene and a 400 bp region of the nuclear second internal transcribed spacer region (ITS-2) gene were amplified using polymerase chain reaction (PCR). We utilized the primers C1-J-1490 (5′-GGTCAACAAATCATAAAGATATTGG-3′; Folmer et al., 1994) and C1-N-2191 (5′-CCCGGTAAAATTAAAATATAAACTTC-3′; Simon et al., 1994) for the cox1 and the ITS-5.8s (5′-GGGACGATGAAGAACGGAGC-3′) and the ITS-28s (5′-TCCTCCGCTTATTGATATGC-3′) for the ITS-2 (White et al., 1990).

PCR amplifications were carried out in 25 µL reaction volume in a final concentration of 0.1 µL Taq polymerase (Promega), 5 µL buffer (Promega), 2.25 µL MgCl2 (Promega), 0.2 mm of each dNTP, 0.5 µL of each primer and 1.5 µL of DNA sample. PCR conditions for amplification were as follows: initial denaturing step at 95 °C for 5 min, 35 amplification cycles (94 °C for 30 s, 45 °C for 35 s, 72 °C for 45 s cox1 fragment and 94 °C for 45 s, 48 °C for 1 min, 72 °C for 60 s for ITS-2 fragment) fallowed by a final extension at 72 °C for 5 min. For certain populations of Pimoa rupicola (localities #23 and #26), a slightly different annealing protocol for the cox1 was utilized (94 °C for 30 s, 42 °C for 35 s, 72 °C for 45 s). PCR products were visualized on agarose gels.

PCR product were cycle-sequenced at Macrogen, Inc. (Seoul, Korea; http://www.macrogen.com). The DNA sequences obtained were preliminary assembled and edited using Geneious 7.1 (Kearse et al., 2012; http://www.geneious.com).

The alignment of the cox1 sequences was trivial, as they showed no evidence of indel mutations. The ITS-2 fragments were aligned with the online version of MAFFT (Katoh & Toh, 2008; http://mafft.cbrc.jp), using the Q-INS-I strategy with default options. We explored best partitioning schemes and substitution models simultaneously using PartitionFinder v.1.0.1 (Lanfear et al., 2012) under a Bayesian information criterion (BIC).

Genetic analyses

Ecological niche modeling

We relied on ecological niche modeling (ENM; see, e.g., Elith et al., 2006) to model the ancestral distribution of the target species. Detailed methodological protocol is provided in Supplemental Information 1. In a first step, we collected all records of the target species available in the literature. We managed to track down 22 localities for Troglohyphantes vignai, most of which clustered together. On the other hand, for Pimoa we recovered 110 localities (61 for Pimoa n. sp. and 49 for P. rupicola), including new unpublished records discovered during the present study. Given the low number of localities for T. vignai, we only inferred the ENM model for the Pimoa lineages. The dataset was corrected for potential spatial autocorrelation and haphazard sampling (Oliveira et al., 2014 and references therein). We obtained present day climatic data (19 ‘Bioclim variables’, see Table 1 of Supplemental Information 1) and altitude a.s.l. from the WorldClim website (www.worldclim.org). We obtained downscaled and calibrated Paleoclimatic data for the Last Glacial Maximum (∼22,000 years ago; hereinafter LGM) from three different simulations available from Global Climate Models (GCMs; Coupled Model Intercomparison Project phase 5; http://cmip-pcmdi.llnl.gov/cmip5). The climatic preferences of the two species were investigated via Principal Component Analysis (PCA) in the Vegan R package (Oksanen et al., 2013). We investigated collinearity among covariates and obtained a final set of uncorrelated variables (Annual mean temperature (Bio1), Temperature annual range (Bio7) and Mean temperature of the driest quarter (Bio9)).

We generated presence-only models with the Maximum Entropy Distribution Models available in MaxEnt (Phillips, Anderson & Schapire, 2006), as implemented in the dismo R package (Hijmans et al., 2014). Firstly, we computed the models on the present climate and on the occurrence points within the M area (sensu Barve et al., 2011; details in Supplemental Information 1). To generate the prediction, we ran each niche model 20 times using a loop script in R, keeping in all cases a random partition of 20% of the occurrence points, which was used to evaluated model performance via the Area Under the Curve (AUC) of the Receiver Operating Characteristic (ROC) plot (Fielding & Bell, 1997). We projected the MaxEnt models in the LGM climate, under each of the three GCMs climatic scenarios. We used a conservative approach to identify potential Pleistocene refugia: we first applied a threshold of 0.6 to the continuous probability surface of presence estimated after the projection. We then combined the three projections, by sub-sampling only those pixel classified as potentially occupied (p > 0.6) in each LGM forecast.

Population structure

The new sequences obtained in the present study are available in GenBank (Supplemental Information 2). A fragment of the mitochondrial cox1 gene of 676 bp was obtained for 37 specimens of Troglohyphantes vignai in 8 localities, corresponding to 14 unique haplotypes. The cox1 data set had 79 segregating sites and 9 parsimony informative sites. The overall mean distance (p-distance among haplotypes) was 0.0495 ± 0.0059. Sequences of the nuclear intron ITS-2 were obtained from the same individuals. The alignment was 400 positions long, 10 additional absence/presence gap characters were scored, corresponding to 10 ITS-2 sequence types. The ITS-2 had 16 segregating sites and 9 parsimony informative sites. The overall p-distance among caves was 0.0389 ± 0.0020. We obtained 676 bp cox1 sequence fragments of 119 Pimoa individuals from 25 localities. The 93 individuals from 19 localities of Pimoa n. sp. yielded 43 haplotypes (35 segregating sites and 7 parsimony informative sites) and the 27 individuals from 6 localities of P. rupicola yielded 7 haplotypes. The average p-distance within populations in Pimoa n. sp. and P. rupicola was 0.0076 ± 0.0017 and 0.0052 ± 0.0017, respectively, and the maximum p-distance between the two lineages was 0.1164 ± 0.0111. The nuclear ITS-2 sequences were obtained from 118 Pimoa specimens. The final alignment included 411 positions and 4 additional gap characters. Individuals of Pimoa n. sp. (90 individuals) and P. rupicola (28 individuals) yielded 34 and 10 sequence types, respectively. The average p-distance within populations in Pimoa n. sp. and P. rupicola was 0.0102 ± 0.0027 and 0.0035 ± 0.0017, respectively, and the maximum p-distance between the two lineages was 0.0701 ± 0.0121.

The standard genetic diversity indices calculated for the cox1 and the ITS-2 for each locality are summarized in Table 2. Pimoa n. sp. showed high levels of nucleotide diversity and low levels of haplotype diversity in most of the populations, while both P. rupicola and T. vignai showed low levels of haplotype and nucleotide diversity. This pattern was especially obvious in P. rupicola, since most individuals within populations were identical.

Pairwise F_ST values calculated for the localities of the three species are reported in Table 3. Pairwise F_ST values between localities revealed contrasting patterns of gene flow. In T. vignai, pairwise F_ST values between localities were always higher than 0.8, and significant (p < 0.05), except for localities #14 and #16. A relatively strong population structure was also found in P. rupicola. Pairwise F_ST values between localities were always higher than 0.6, although significant comparisons involved exclusively southernmost localities (#28 and #29). Pimoa n. sp. showed instead a more shallow population structure, with several F_ST values below 0.5, generally corresponding to nearby localities.

Pimoa haplotypes were resolved as two independent networks, corresponding to Pimoa n. sp. (1) and P. rupicola (2), respectively, separated by 42 steps. P. rupicola haplotypes were limited to single populations, except H26 and H27, which were shared across populations. Generally, populations had low haplotype diversity, and in one case (#26) there was one single haplotype. Two divergent haplotypes (11 steps), however, were found in locality #29.

In Pimoa n. sp. haplotypes from nearby localities clustered together or were separated by only few steps (1–4). Several haplotypes (H7, H8, H12, H14) were shared among closely located localities (e.g., occurring in the same Alpine valley or in adjacent valleys). In some instances, however, haplotypes from distant localities were found to be very similar (few steps), e.g., locality #1 and #2 with #12. Moreover, two haplotypes (H18 and H22) were found in individuals occurring in distant populations.

In T. vignai, haplotypes were not shared between localities. Haplotypes from closely located populations were generally separated by few mutations. Haplotype diversity was low within populations, and localities #7, #14 and #16 showed single haplotypes. Two localities (#9 and #25), on the other hand, had more than two haplotypes. The haplotypes of T. vignai sensu rupicapra were separated by 17–20 steps from the nearest locality (#14), which is actually the locus typicus of T. vignai (Buco di Valenza cave). A higher number of steps (38–56) separated this cave from the remaining localities.

Phylogenetic tree and estimation of the divergent time

Partition Finder selected the full codon as the best partition scheme for the alignments of both species. The models for each partition are reported in Table 2 of Supplemental Information 1. Both the Bayesian and ML analyses of the concatenated data matrix of T. vignai yielded in similar tree topologies, and most branches were highly supported (i.e., posterior probabilities (PP) > 0.95, bootstrap support (BS) < 75%; Fig. 2). T. vignai was split in two main clades: one including the southern populations (#14, #16, #25 and #27), and the second one including the remaining northern populations (#5, #6, #7 and #9). The Bayesian and the ML analyses also recovered similar tree topologies in Pimoa (Fig. 3). Two well-supported clades were detected, corresponding to Pimoa n. sp. and P. rupicola, respectively. Individuals from geographical adjacent localities were closely related, although basal branches within Pimoa n. sp. were poorly supported.

The GMYC algorithm identified 2 coalescent clusters within Pimoa cox1 sequences (ML = 100.8932; LR = 28.59767; p < 0.000), one including all sequences/localities of P. rupicola and the other including all sequences/localities of Pimoa n. sp. Troglohyphantes vignai cox1 sequences were resolved as 7 coalescent clusters (ML = 424.2824; LR = 207.4615; p < 0.000; Fig. 4). Except for cluster A2, which included individuals from two caves, each cluster corresponded to individuals sampled from single caves.

The species trees and the embedded cox1 gene tree recovered for each spider genus are shown in Fig. 4. The substitution rate estimated for the Troglohyphantes cox1 was 0.0218 substitutions per lineage/million years (95% HPD = 0.010–0.033), and for the ITS-2 was 0.0024 (95% HPD = 0.0005–0.0031). The split between T. vignai and T. nigraerosae was traced back to approximately 7.2 million years ago (Ma, 95% Highest Posterior Density, HPD = 13.7–3.5 Ma). Diversification of the extant T. vignai lineages occurred 2.9 million years ago (95% HPD = 5.4–1.5 Ma), while the diversification of the extant northern populations (D2, E2, F2, G2 clusters) occurred approximately 0.5 Ma (95% HPD = 0.9–0.2 Ma). The substitution rate estimated for Pimoa cox1 was 0.0217 substitutions per lineage/million years (95% HPD = 0.011–0.033), and for the ITS-2 was 0.006 (95% HPD = 0.002–0.013). The basal split between Pimoa n. sp. and P. rupicola was estimated to have occurred around 5.7 Ma (HPD = 12–2 Ma). The origin of the extant diversity of each lineage was estimated approximately at 0.4 Ma for both lineages (HPD = 1–0.15 and HPD = 0.85–0.15 Ma, for P. rupicola and Pimoa n. sp. respectively).

Climatic segregation

We studied the climatic preference for the two Pimoa linages with PCA, using bioclimatic variables (see Table 1 of Supplemental Information 1). The bi-plot of scores for the first two axes of the PCA is shown in Fig. 5. The first two axes explained 86.2% of the variation in the data. The variance explained by other axis was negligible (<1% each). In respect to the first two axes, the two species segregate into two distinct clusters. The first axis (eigenvalues = 12.28; variance explained = 62.0%) mostly reflect a gradient of temperature. The second axis (eigenvalues = 4.79; variance explained = 24.2%) was positively correlated with variables reflecting diurnal and annual thermic excursion (Bio2, Bio4, Bio7) and anticorrelated with variables reflecting seasonal precipitation (Bio13, Bio16, Bio19). Although the first axis (PC1) explains most of the variance in the data, the localities of the two Pimoa species cluster in two groups according to the second axis (PC2), which combines bioclimatic variables referable to continentality. According to the original definition (see Currey, 1974 for more details), continentality is intended as a measure of how the climate is affected by its remoteness from the sea. Specifically, the distance from water masses influences the climate in terms of higher seasonality (Bio4), increasing diurnal (Bio2) and annual (Bio7) temperature ranges, as well as decreasing precipitation in the coldest (Bio19) and wettest (Bio13 and 16) periods. In light of our results, Pimoa n. sp. occurs in areas characterized by higher continentality in respect to P. rupicola.

Current and past distribution of Alpine Pimoa

The predictive performance of our bioclimatic models was fairly high both in Pimoa n. sp. (mean ± SD AUC of the 20 runs = 0.845 ± 0.053) and in P. rupicola (0.908 ± 0.089). Overall, the suitable areas predicted by the model were congruent with the known distribution of the two species, at least in the Western Alps and Apennines. Current predictions identified suitable areas for P. n. sp. around the medium mountain belt (500–1500 m a.s.l.) of the Central and Western Alps, from the Camonica valley (province of Bergamo) down to the margin between Cottian and Maritime Alps (Stura valley, province of Cuneo). In respect to the known distribution, the predicted range of P. n. sp. extended northwards over the known limit of the species (see dotted line in Fig. 6A). More suitable areas were also detected in the northern edge of the Tuscan Apennines (Fig. 6A). The suitable areas of P. rupicola corresponded to the southern border of the Western Alps, in the coastal belt that spreads from Côte d’Azur (SW France) to the Ligurian eastern coast and Tuscan Apennines (Italy). Additional suitable areas were also found in Tuscany, Lazio and Corsica. The projection of the distribution model to the environmental condition of the Last Glacial Maximum (LGM) revealed that most of the current suitable areas were unsuitable in the LGM (Fig. 6B). Our threshold approach identified one main refugia for P. n. sp. (RF1) and two for P. rupicola (RF2 and RF3). All refugia corresponded to areas that were devoid from glaciers (in accordance with Ehlers, Gibbard & Hughes, 2011). The RF1 extended outside the southern edge of the actual distribution of Pimoa n. sp., in the hills and plains surrounding roughly the Northern border of Maritime and Ligurian Alps. RF2 corresponded to small areas along the French and Italian Riviera. RF3 extended over a wider geographic area in the Apennine, and in the northern part of the Corsica.

The history of two cave-dwelling spiders

The confounding effects of adaptation, biogeography and dispersal ability on the origin and the distribution patterns of cave organisms pose a stimulating challenge to biogeographers (Porter, 2007; Juan et al., 2010). According to Culver & Pipan (2010), long term climatic changes can be claimed as the main factors that prompted invertebrate species to colonize the subterranean habitats. In this regard, the Miocene climatic transitions and the Pleistocene glaciations are considered among the most important events. Here, by reconciling phylogeographic patterns and predictive ENMs, we provide support to this view, pointing out the Cenozoic climatic transitions has the most important factors shaping the present day genetic diversity and the distribution range of our model species.

Although special caution should be exercised when considering time estimates based on molecular data, especially in the absence of fossil record (Hipsley, 2014), our results fits well with some of the major climatic event undergone in the Western Alps during the Cenozoic. Accordingly, the two Pimoa lineages and Troglohyphantes vignai originated from the Middle (Serravallian) to the Late Miocene (Messinian), namely from 13 to 3.5 mya. More precisely, the isolation of T. vignai and P. rupicola from their northern sister groups (T. nigraerosae and P. n. sp.) dates back 7.2 and 5.7 mya, respectively. This time period approximately corresponds to the closure of the Gibraltar Strait and the onset of the so-called Messinian Salinity Crisis (MSC; after Ruggieri, Adams & Ager, 1967). However, given the large confidence intervals around our time estimates (∼10 million years), it is difficult to draw precise conclusions about the event that exactly determined the split of the different lineages.

It is worth noticing that the onset of a climatic transition in the Middle Miocene, marked the decline of the last global climate optimum conditions, leading to a progressive deterioration of the dominant subtropical climatic conditions (Suc, 1984; Shevenell, Kennett & Lea, 2004; Jiménez-Moreno, Fauquette & Suc, 2010). It is arguable that in parallel to the slow climate deterioration and the increase of seasonality, isolation of Pimoa and Troglohyphantes occurred.

Being possibly pre-adapted to shallow moist humid habitats (Deeleman-Reinhold, 1978; Hormiga, 1994; Zhang & Li, 2013; Wang et al., 2008), both species progressively colonized the subterranean habitat, most likely during the Pleistocene. Given their contrasting level of troglomorphism, it is likely that the process began earlier in Troglohyphantes.

The known distribution range of T. vignai stretches discontinuously along the Cottian (Chisone, Po and Varaita valleys) down to Ligurian (Valle Pesio) Alps. Conversely, the distribution of the northern sister group T. nigraerosae is adjacent southern Graian Alps (Isaia et al., 2011). The lack of shared haplotypes and the F_ST values close to 1 between the sampled localities of T. vignai indicates a strong isolation of the populations (Holsinger & Weir, 2009). The same idea is further corroborated by the identification of each population as an independent coalescent lineage (i.e., GYMC cluster). These results are in agreement with other studies on subterranean arachnids (Hedin, 1997; Hedin & Thomas, 2010; Dixon & Zigler, 2011) that support the “caves as islands” scenario (sensu Snowman, Zigler & Hedin, 2010), in which dispersal is virtually absent and the different populations diverge in allopatry. Under such conditions, it seems likely the subterranean habitat acted as an evolutionary cul-de-sac (see Fišer, Blejec & Trontelj, 2012) for T. vignai. According to our time estimates the diversification of extant T. vignai lineages (especially northern population) occurred approximately during the Pleistocene glaciations. In this respect, it is worth noting that subterranean localities inhabited by T. vignai lie at the periphery of the Pleistocene glaciers (Ehlers, Gibbard & Hughes, 2011; see also local glacial limits reconstructed in Motta, 2014). Because subterranean populations most likely cannot survive under the ice cover (Culver & Pipan, 2010), we suggest that the present day distribution range of T. vignai is the shadow of a wider ancestral distribution. Populations inhabiting the northern valleys in the Cottian Alps (Germanasca and Chisone valleys; see Fig. 4), where the ice shield was more compact, provide further evidences of the effect of Quaternary ice sheets.

In this area the distribution of T. vignai overlaps with the range of several hypogean species of Doderotrechus beetles (Carabidae, Trechini) (Giachino, 1993; Giachino & Vailati, 1997; Casale & Giachino, 2008). The similarity of distribution of both groups is most likely the result of a similar response to glacial dynamics. Gaining further knowledge on the biogeography of other subterranean species may provide further confirmation for the patterns here recovered.

In contrast with T. vignai, for which we did not detect any evidence of current population expansion, the topology of the chronogram obtained for Pimoa lineages hints at a recent expansion following a bottleneck (see Fig. 4). The putative population expansion in both lineages of Pimoa would fall within the Quaternary Glacial Cycles, between 2.8 and 2.5 Ma (Gibbard, Head & Walker, 2010). The movement of glaciers as well as the continuous formation and melting of new ice sheets may have deeply affected the different populations, altering profoundly the local habitat suitability. Such transformations prompted either the migration of populations to more suitable areas at lower altitudes or latitudes or the local extinction of resident populations. Such scenario is congruent with the genetic fingerprint found in both Pimoa lineages. Accordingly, we hypothesize a glacial cycle-driven extinction of ancestral populations during cooler periods, followed by the expansion of populations which survived in climatic refugia during warmer periods. A similar pattern was observed in the Asian species P. clavata (Wang et al., 2008). The ENMs projected into the paleoclimatic reconstruction pointed out three putative areas devoid from glaciers that may have acted as glacial refugia for the surviving populations of each lineage (RF1-3; Fig. 6B). In the case of P. n. sp., we detected one main macrorefugia (RF1) associated to the southernmost offshoots of the Alpine glacial masses. Notably, similar areas have also been classified as peripheral refugia for several Alpine plants, such as Phyteuma globulariifolium (Schönswetter et al., 2002) and Ranunculus glacialis (Schönswetter et al., 2003). For P. rupicola, potential macrorefugia where located along the SW French coast (RF2) and the Tuscan coast (RF3), including the northern part of Corsica. Even though the presence of P. rupicola in Tuscany is confirmed by literature records (Brignoli, 1971; Hormiga, 1994), recent investigation conducted by the authors did not confirm the present occurrence of the species in Corsica. Because of the larger spatial resolution of the LGM stacked rasters compared to the present day data (2.5 min versus 30 arc-seconds (∼1 km)), it should be borne in mind that we may have not detected small point-like microrefugia (sensu Rull, 2009) within the interior of the Pleistocene ice shield covering the Alps.

As suggested by the lower F_ST values compared to Troglohyphantes, the relatively fast recolonization of ice-free areas is probably the result of the more effective dispersal ability of Pimoa. In particular, the population expansion followed a south-north direction, leading to the present distribution ranges of both lineages. Concerning P. n. sp., we additionally predict suitable areas up to the Central Alps. Given the continuity of the suitable habitats predicted by the ENM and the supposed high dispersal ability of Pimoa, we hypothesize an ongoing expansion of the populations northwards. Indeed, the occurrence of Pimoa in outer shaded and humid habitats such as beech forests and other broadleaved forests (Bertkau, 1890; Jackson, 1926; Thaler, 1976; Isaia et al., 2015; Isaia, Paschetta & Chiarle, 2015) provides empirical evidence of the existence of epigean dispersal. Because of the sex bias among the specimens collected in superficial habitats (Isaia et al., 2015; Isaia, Paschetta & Chiarle, 2015) and the general trend observed for spiders (Foelix, 1996), gene flow appears mostly mediated by males.

It should be borne in mind, however, that in light of the relatively low sample-size and the potential bias linked to the possible male-mediated dispersal, our cox1 data (and the associated F_ST results) may not fully reflect patterns of gene flow (see, e.g., Willing, Dreyer & Van Oosterhout, 2012; Davalos & Russell, 2014). Since low sample size is known to impact this kind of estimations, caution should be exercised when interpreting the values of nucleotide and haplotype diversity calculated for each locality (e.g., Goodall-Copestake, Tarling & Murphy, 2012). To minimize these potential bias, we left out from the calculation of the F_ST populations consisting of less than three individuals. Nevertheless, it is worth noticing that during data exploration we obtained very similar patterns of gene flow when comparing cox1 and ITS-2 results.

At present, the two lineages of Pimoa identified in this study show allopatric distributions (Fig. 6). P. n. sp. populations occur preferentially in areas characterized by higher continentality, and seem to tolerate cooler temperatures at higher altitudes and latitudes, as suggested by the results of the PCA (Fig. 5). On the other hand, Pimoa rupicola occurs in less continental areas, characterized by relatively small seasonal variations and high mean annual temperatures (i.e., Mediterranean climate). Similar complete niche partitioning between congeneric subterranean spiders has been reported elsewhere (Ribera, 1978; Gasparo & Thaler, 2000; Mammola & Isaia, 2014).

The application of ENM techniques has become a widespread practice to answer biogeographical and evolutionary questions (Franklin, 2009). In particular, ENM have been extensively used to identify Pleistocene refugia (e.g., Waltari et al., 2007; Rodriguez-Sanchez & Arroyo, 2008; Peterson, 2009; Planas et al., 2014). In constructing our scenarios, we have adopted a conservative approach, as our goal was to generate predictions under different levels of uncertainty. Although we relied on this approach, we are aware that ENMs have been rarely—and only recently—applied to study the hypogean ecosystems (see, e.g., Bryson et al., 2014; Camp et al., 2014; Naranjo, Moreno & Martiín, 2014). This is probably because, in first approximation, the link between the climatic variables (i.e., the external climate) and the subterranean habitat is not so straightforward. However, temperature of the underground compartment generally reflects the climatic regimen on the surface (Smithson, 1991; Badino, 2010). Although less intuitive, the regimen of rainfall plays an equally crucial role—if not more important—in determining such conditions (see details in Badino (2004) and Badino (2010)).

Overlooked diversity

In light of our results, some consideration regarding the overlooked diversity of our model species can be drawn. Concerning T. vignai, in this study we have included specimens from eight different localities, including topotypical material and material of T. vignai sensu rupicapra. The low levels of genetic variability observed between the latter and the topotypic material of T. vignai (Fig. 2, p-distance = 0.0022), provide further support for the synonymy between the two species proposed by Pesarini (2001). Therefore specimens of T. vignai sensu rupicapra have to be regarded as a population of T. vignai isolated in an area characterized by different climatic conditions (major Mediterranean influence). Such isolation could tentatively be related to the higher development of troglomorphism in T. vignai sensu rupicapra, as already observed by Brignoli (1971) and Isaia & Pantini (2010). Moreover, climatic factors provide a further line of interpretation for the presence of two main lineages within T. vignai, corresponding to the northern and the southern clade (Figs. 1 and 2). Slight differences in the shape of the lamella significativa of the male palp could lead to consider the two lineage as candidate species (Vieites et al., 2009), however the genetic distance between the two lineages, is about half of the value observed in the two nominal species T. vignai and T. nigraerosae.

The genus Pimoa is represented worldwide by 28 species (World Spider Catalog, 2015), many of which have only been described recently (Xu & Li, 2007; Xu & Li, 2009; Trotta, 2009; Hormiga & Lew, 2014). The application of molecular tools to investigate fine scale phylogeographic patterns in this group may uncover additional hidden diversity (Wang et al., 2008). It is generally accepted that species delimitation and eventually species description should be based on the integration of multiple lines of evidence (Padial et al., 2010). Here, we uncovered two deeply divergent genetic lineages (Fig. 3, GYMC clusters with p-distance above 9%) within Pimoa rupicola, which are further delimited both by genitalic morphology (S Mammola, G Hormiga, MA Arnedo, M Isaia, 2015, unpublished data) and different ecological requirements (Fig. 5).