Dispersal Capabilities Do Not Parallel Ecology and Cryptic Speciation in European Cheliferidae Pseudoscorpions (Pseudoscorpiones: Cheliferidae)

Abstract

The ability to disperse has continually shaped both the distribution and diversification of biota, and it affects the survival of the species with respect to wide-ranging habitat loss. As a response, organisms unable to spread by their own means often developed surrogate dispersal strategies. Pseudoscorpions possess small body sizes and cannot actively disperse over large distances and geographic barriers; therefore, they have adopted other ecological strategies. They are either sedentary and remain confined to stable environments or passively disperse via phoresy and are capable of inhabiting a wide variety of habitats, including temporary ones. In this paper, we use barcoding data to investigate the genetic diversity of four widely distributed and relatively morphologically uniform Cheliferidae genera Chelifer, Dactylochelifer, Rhacochelifer and Hysterochelifer. We aim to (i) test whether the genera harbor cryptic diversity and (ii) evaluate whether the genetic structure of the species parallels their dispersal capabilities and habitat preferences (i.e., ecological strategies). In general, we uncovered independent lineages within all analyzed genera, which indicates their need for a thorough and integrative taxonomic revision. More specifically, we detected a varying degree of genetic structuring among the lineages. Known phoretic species, as well as some species and delimited lineages that are not known to use this manner of dispersal, showed a complete lack of geographical structure and shared haplotypes over large distances, while other taxa had restricted distributions. We argue that genetic structure can be used as a proxy to evaluate species’ dispersal manner and efficacy. Our results also suggest that taxa inhabiting stable environments might use phoresy for their dispersal.

1. Introduction

One of the central topics in biodiversity research is to uncover the underlying evolutionary mechanisms that created the current distributions of biota [1,2]. Congruencies among phylogeographic patterns in unrelated taxa that share a similar level of dispersal capability can, for example, inform us of past geological and climatic events that played a key role within a geographic region [3,4,5]. As such, the capability to disperse (or the lack thereof) has continually shaped taxa distributions [6,7,8] and critically affects the survival of species in the case of wide-ranging habitat loss [9]. Organisms unable to disperse by their own means often developed passive methods of dispersal, such as the wind-mediated dispersal of eggs and dispersal via attachment to a more vagile carrier, in order to compensate for their small body size or sedentary lifestyle [10,11,12]. However, relatively little is known about the genetic structure of arthropods that rely on other organisms for their own dispersal [13,14,15,16].

Pseudoscorpions, with about 4000 described species [17], belong among smaller arachnid orders that have historically been neglected due to their morphological uniformity, small body size and rather cryptic lifestyle [18]. Due to their small size, pseudoscorpions typically cannot actively disperse over large distances and geographical barriers. As a result, pseudoscorpion taxa adopted one of the two distinct ecological strategies: humicolous species mainly inhabit soil and leaf litter, while arboricolous species are associated with bark, tree hollows and other temporary habitats [19,20,21]. The strategies also presumably have an opposite impact on the diversification and dispersal capabilities of the taxa. The soil-dwelling species, or species inhabiting isolated microhabitats with stable conditions such as caves, are usually sedentary. Consequently, they show high tendencies toward microendemism [22,23,24,25] and might be highly sensitive to habitat degradation [26]. On the other hand, pseudoscorpion species that inhabit temporary and unstable habitats (e.g., tree hollows, bird nests, decaying plant material, insect galls, etc.) were able to surpass their low vagility and disperse via phoresy, i.e., hitchhiking on the body of a more vagile carrier [27], because of the strong evolutionary pressure to disperse in order to find other suitable habitats. Pseudoscorpions’ phoretic hosts are typically a larger and mobile insect, which naturally exists in the same type of patchy microhabitat [27,28,29]. As a result, phoretic species typically have greater distribution ranges [13,16,30,31,32]. Additionally, large distributions are known from synanthropic species that at least partially owe their dispersal to human activity, e.g., transport of goods and landscape transformation. Such distributions may even span across different continents [33].

The assessment of pseudoscorpion diversity and phylogenetic relationships within the order Pseudoscorpiones has traditionally been studied via morphological approaches [33,34]. However, the evaluation of morphological characters and their interpretation in an evolutionary context are challenging because of the compounded effect of the overall external uniformity of the group and our poor understanding of intra- and interspecific variability in the given characters [35,36]. As a result, species routinely oscillated among genera and sometimes even families [17,37]. Our estimation of pseudoscorpion diversity is also likely to be skewed for the same reasons. The notion that morphological approaches typically underestimate the real number of species in morphologically uniform taxa [38] evidently also applies to pseudoscorpions. Molecular methods revealed the presence of cryptic diversity in many genera [25,39,40,41,42], including those from well-studied regions, such as Central and Western Europe [16,21]. Additionally, the implementation of genetic/genomic data [43,44] provided stability to phylogenetic relationships and created the necessary backbone for the targeted investigation of ecologic and evolutionary patterns in pseudoscorpions [13,16,21,45,46,47,48].

In recent years, DNA barcoding has gained popularity as an accessible tool for both specimen identification and species discovery [49,50]. As a result, public databases, such as BOLD [51] and GenBank [52], have amassed large amounts of barcoding data for a broad range of taxa. Such availability facilitates specimen identification via comparison with known reference sequences [53], useful in wide array of applications e.g., environmental DNA analyses [54,55], prey/food analysis [56,57] and linking different developmental stages, sexes or castes within one species [58,59,60]. Barcoding is commonly used for regional diversity assessment [61,62], producing genetic data for widely distributed species, thus allowing the large-scale assessment of population structure and distribution patterns and identification of potential cryptic diversity [63,64,65,66].

In this paper, we investigate the genetic structures of four widely distributed European genera of pseudoscorpions (Chelifer, Dactylochelifer, Rhacochelifer and Hysterochelifer) from the family Cheliferidae. Cheliferidae is cosmopolitan and one of the most diversified pseudoscorpion families (312 described species; [17]); however, its diversity and intrafamiliar relationships remain poorly understood. Cheliferidae have never been sufficiently sampled in any phylogenetic analyses [43,44], and their taxonomy frequently relies on male-specific morphological characters [67]. Interestingly, the species of our target genera show significant differences in terms of the extent of their distributional ranges and habitat preferences [17]. While some species are known to be phoretic [27], as well as synathropic in case of Chelifer cancroides [35], the dispersal capabilities of most taxa are not known. We combine barcoding data [cytochrome c oxidase subunit 1 (Cox1) [68,69]], proceeding both from our targeted sampling and public databases, in order to (i) test whether the genetic structure reflects habitat preferences and dispersal capabilities (i.e., the opposite ecological strategies) and (ii) whether the genera in question harbor cryptic diversity.

2. Materials and Methods

2.1. Taxonomic Sampling

We analyzed 217 samples of pseudoscorpions from the family Cheliferidae collected mostly across Europe and adjacent areas (Table S1). We predominantly targeted common genera, namely Chelifer, Dactylochelifer, Hysterochelifer and Rhacochelifer, in order to assess distribution patterns and genetic structures within their widely distributed species and test for the potential presence of cryptic diversity. The specimens were collected via leaf litter sifting or hand-collected under tree bark, under stones, in tree hollows, in galls and in synanthropic environments. The intensity of sampling and the used collection methods were similar at each site; all specimens have been included in the analyses. However, none of the individuals was collected during phoresy. All samples were stored in 96% ethanol at −20 °C for future DNA analyses. Additional Cheliferidae Cox1 sequences (62 in total) were downloaded from the GenBank and BOLD databases (Table S1) to both broaden geographic sampling of our target groups and increase family-level taxon sampling. Specimens of Cheiridium museorum (Cheirididae), Titanatemnus natalensis (Atemnidae), Chernes hahnii and C. nigrimanus (both Chernetidae) were used as outgroups in our analyses, representing pseudoscorpion families recovered as the closest relatives of Cheliferidae [44]. All distribution maps were created via SimpleMappr [70], and geographic distances among sample localities were calculated via Geographic Distance Matrix Generator [71].

2.2. DNA Extraction and PCR Amplification

Genomic DNA was extracted from the whole specimens using a Geneaid Genomic DNA Extraction Kit (Tissue) (Taipei, Taiwan) following manufacturer’s protocol. To maximize the extraction yield, specimens were carefully opened either in the pleural region or between the carapace and first tergite. The vouchers were subsequently returned to 96% ethanol for future morphological determination and stored in the collections of the Department of Zoology, Faculty of Science, Charles University, Prague. A fragment of mitochondrial protein coding cytochrome c oxidase I gene (Cox1, the animal barcode) was amplified via PCR (95 °C 3 min; 35 cycles of 95 °C 30 s, 45 °C 1 min, 72 °C 1.5 min; 72 °C 10 min) using primer pair C1-J-1490/C1-N-2198 [72]. The PCR products were Sanger-sequenced using the Macrogen (Amsterdam, Netherlands) and Biocev (Vestec, Czech Republic) sequencing facilities.

2.3. Phylogenetic Analyses

The chromatograms were assembled, edited and checked for stop codons in Geneious v. R9 [73]. The sequences were subsequently aligned using the MUSCLE [74] alignment algorithm included among the Geneious built-in options. The resulting matrix (“Cheliferidae_out_total”) comprising all sequences was reduced to unique haplotypes in TCS [75]. The reduced matrix comprising 197 unique haplotypes (“Cheliferidae_out_haplo”) was used for all downstream phylogenetic analyses. The best partitioning scheme (Table S2) was determined via PartitionFinder v. 2.1 [76] using the greedy algorithm and BIC selection criterion. Maximum likelihood (ML) analyses were conducted via RAxML v. 8.2 [77]. The GTRGAMMA evolutionary model was assigned to each partition, and the best tree was selected from 1000 iterations. Nodal support was recovered from 1000 replicates of bootstrap resampling. Bayesian inference (BI) was performed using BEAST 2.5.2 [78]. Two independent runs, each of 10⁸ MCMC (Markow chain Monte Carlo) generations, were remotely run via CIPRES Science Gateway [79]. The input file, created in BEAUTi, was set with a single GTR+G partition and a strict clock in order to facilitate convergence and speed up the analyses. The convergence and effective sample size (ESS) was checked via Tracer v. 1.7.2 [80]; 25% of the samples was discarded as a burn-in for the analyses. The maximum clade credibility tree was inferred via TreeAnnotator v. 2.6.6. All of the resulting trees were visualized and manipulated via FigTree v. 1.4.4 [81]

2.4. Molecular Species Delimitation

Four genera with wide distributions (Chelifer, Dactylochelifer, Hysterochelifer and Rhacochelifer) were subjected to species delimitation analyses in order to assess (i) whether species boundaries were congruent with the morphological characterization of given taxa and (ii) whether widely spread species potentially harbored cryptic diversity. Individual matrices (“Chelifer_out”, “Dactylo_out”, “Hystero_out” and “Rhaco_out”) comprising the unique haplotypes of each taxon and selected outgroups (C. hahnii + D. latreillii for “Chelifer_out” and “Hystero_out” matrixes, and C. hahnii + C. cancroides for “Dactylo_out” and “Rhaco_out” matrixes) were created for each taxon.

We implemented four independent species-discovery approaches. Firstly, we used Assemble Species via Automatic Partitioning (ASAP) [82], which is a method that operates using hierarchical clustering principles. The analyses were performed via the ASAP delimitation server (https://bioinfo.mnhn.fr/abi/public/asap/ [accessed on 2 September 2022]) using both simple p-distances and the Kimura (K80) substitution model [83]. Secondly, we used statistical parsimony (SP) haplotype network analysis in TCS v 1.21 [75]; putative species were delimited under the 95% parsimony criterion. Finally, we used two independent tree-based methods: the Poisson Tree Process (PTP) [84] and the Generalized Mixed Yule coalescent analysis [85]. The PTP analyses were remotely run via the bPTP webserver (https://species.h-its.org/ptp/ [accessed on 15 October 2022]) for 5 × 10⁵ generations. The topologies for the PTP analyses were obtained via RAxML analyses carried out on taxa-specific matrices (“Chelifer_out”, “Dactylo_out”, “Hystero_out” and “Rhaco_out”) with the same parameter settings as the phylogenetic analyses performed on “Cheliferidae_out_haplo” (see above). The best partitioning schemes for each codon position were determined via PartitionFinder (Table S2). The outgroups were removed from each topology prior to the delimitation analyses. The GMYC analyses with a single threshold were carried out for taxa-specific matrices in the R environment (http://www.r-project.org [accessed on 15 February 2023]) using the SPLITS package [86]. The topologies were obtained via BEAST analyses performed on taxa-specific matrices with removed outgroups. Each analysis comprised two independent runs of 5 × 10⁷ generations, remotely performed via CIPRES Science Gateway and by implementing the best partitioning schemes for each codon position determined via PartitionFinder (Table S2). Delimitation analyses yielding lower numbers of putative species were preferred in the results interpretation stage (see Section 4.2).

Inter- and intraspecific genetic distances among the delimited Chelifer, Dactylochelifer, Hysterochelifer and Rhacochelifer species were calculated using Mega 11 [87] via both uncorrelated p-distances and the Tamura–Nei distance model [88].

3. Results

3.1. Taxonomic Sampling and Phylogenetic Analyses

A total of 282 specimens (278 Cheliferidae samples, four outgroups) representing Cheliferidae species worldwide was used in this study (Table S1). The Cox1 fragment (654 bp) was newly acquired for 218 individuals, while the remaining 64 sequences were obtained from public databases (Table S1). The “Cheliferidae_out_haplo” matrix comprising 198 unique haplotypes (193 ingroup and 5 outgroup taxa) was used as an input for ML and BI analyses. The matrix was further used as a source of reduced individual taxon-specific matrices used in species delimitation. The partitions and evolutionary substitution models selected via PartitionFinder are reported in Table S2.

The phylogenetic analyses performed via ML (−lnL= −15,350.541) and BI yielded topologies that lacked support for deeper nodes (Figure 1 and Figure S1); however, the supported nodes (mostly at crown level) were topologically congruent between the results of both approaches. Overall, the topology recovered in BI was better supported. In both analyses, Cheliferoidea outgroups (T. natalensis, C. hahnii and C. nigrimanus) formed a clade sister to the monophyletic Cheliferidae. Nannochelifer and Protochelifer were recovered in a clade (supported in both analyses) sister to all of the remaining Cheliferidae taxa. Within the remaining Cheliferids, the relationships between genera were largely unresolved; all genera, with the exception of Beierius and Philomaoria, formed a clade supported in BI. The sister relationships of Beierochelifer and Hansenius, Parachelifer and Mesochelifer and Socio-chelifer and Chelifer were supported in BI, which also supported a clade formed by Parachelifer, Mesochelifer, Sociochelifer, Chelifer and Hysterochelifer. All genera included in the analyses, with the exception of Mesochelifer, were supported at least in the BI (Figure 1 and Figure S1).

3.2. Molecular Species Delimitation and Geographic Distribution of the Lineages

3.2.1. Chelifer

Our “Chelifer_out” matrix comprised 21 unique haplotypes of 32 individuals, all identified as C. cancroides (Figure 2) and two outgroups (C. hahnii and D. latreillii). The results of our species delimitation analyses were not congruent between the delimitation approaches, with the number of recovered independent lineages ranging from three to six (Figure 2). The lowest number of lineages was recovered via ASAP1. The method delimited three independent lineages at the first ranked parsimony scheme (threshold distance 0.067). The second ranked ASAP2 scheme (threshold distance 0.020) recovered six evolutionary lineages congruent with the results of statistical parsimony (SP) performed in TCS. The tree-based methods (PTP and GMYC) recovered between five and six independent lineages, respectively. However, the distribution between taxa slightly differed from the SP results. The independent lineages delimited from our Chelifer dataset were spatially overlapping and did not show any intra-lineage geographical structure (Figure 2). The most sample-rich evolutionary lineage (lineage C in ASAP1 top ranked scheme) contained 27 individuals (16 unique haplotypes) distributed across Europe and North America. The largest distance between individuals sharing the same haplotype was 6837 km (Table S3).

3.2.2. Dactylochelifer

Our Dactylochelifer dataset (“Dactylo_out” matrix) comprised 41 unique haplotypes proceeding from 81 individuals. The species delimitation results varied (8–18 independent lineages) depending on the implemented method (Figure 3). The lowest number of lineages was recovered via ASAP2. The second ranked scheme (threshold distance 0.067) identified eight evolutionary lineages, while the top ranked scheme (threshold distance 0.039) yielded ten lineages. Twelve evolutionary lineages were recovered via SP (95% parsimony criterion) performed in TCS. The tree-based PTP approach yielded congruent results via ASAP1, i.e., ten independent lineages. GMYC delimited the highest number of lineages (18).

The majority of our delimited lineages was detected in relatively large geographic areas; the three most densely sampled lineages belonging to D. latreillii morphotype (F, G, H) showed large and geographically overlapping distributions across Europe and Western Asia with no apparent geographic structure. The largest geographic distance (1443 km) between individuals sharing the same haplotype was detected in lineage F (Table S3).

3.2.3. Hysterochelifer

The “Hystero_out” matrix comprised 61 unique haplotypes corresponding to 77 Hysterochelifer individuals. Nine to seventeen independent lineages were identified via the species delimitation process (Figure 4). The results from the top ranked ASAP1 score (threshold distance 0.062) indicated the existence of nine independent lineages (Figure 4), while second ranked ASAP2 score (threshold distance 0.032) detected 13 lineages. Statistical parsimony (95% parsimony criterion) performed in TCS delimited 17 lineages. The PTP approach delimited 15 lineages, while the GMYC approach yielded 13 lineages congruent with the results of the second ranked ASAP score.

Two distributional patterns were detected within the Hysterochelifer diversity. Some independent lineages exhibited local endemism to some extent (e.g., lineages B, C and D), while the lineages A (H. meridianus morphotype) and F (H. tuberculatus morphotype) had wider distributions that often overlapped with the distributions of other, more geographically isolated lineages. The distribution area of lineage A (H. meridianus morphotype), for example, spanned from eastern Spain to eastern Turkey, and the largest distance between individuals sharing the same haplotype was 851 km (Table S3).

Figure 2. Chelifer cancroides species delimitation. (A) Phylogenetic tree of Chelifer cancroides with the results of different species delimitation approaches (right). Topology was obtained via the BI conducted in BEAST, the numbers on the nodes denote PP support values, and the nodes with PP values > 0.95 are marked with black circles. White star marks a node that had a different topology in the ML analyses used for species delimitation in PTP. Terminal tree taxa are color-coded according to the country of their collection. Bottom left corner insert shows an adult male of C. cancroides; scale bar: 1 mm. (B) Map depicting the sampling locations of C. cancroides delimited lineages according to our preferred species delimitation outcome obtained via ASAP1 (upper right). Top left corner of the map shows the sampling locations in North America.

Figure 2. Chelifer cancroides species delimitation. (A) Phylogenetic tree of Chelifer cancroides with the results of different species delimitation approaches (right). Topology was obtained via the BI conducted in BEAST, the numbers on the nodes denote PP support values, and the nodes with PP values > 0.95 are marked with black circles. White star marks a node that had a different topology in the ML analyses used for species delimitation in PTP. Terminal tree taxa are color-coded according to the country of their collection. Bottom left corner insert shows an adult male of C. cancroides; scale bar: 1 mm. (B) Map depicting the sampling locations of C. cancroides delimited lineages according to our preferred species delimitation outcome obtained via ASAP1 (upper right). Top left corner of the map shows the sampling locations in North America.

Figure 3. Dactylochelifer species delimitation. (A) Phylogenetic tree of the genus Dactylochelifer with the results of different species delimitation approaches (right). Topology was obtained via the BI conducted in BEAST, the numbers on the nodes denote PP values, and the nodes with PP values > 0.95 are marked with black circles. Terminal tree taxa are color-coded according to the country of their collection. Specimen DL2042 (marked with asterisk) lacks locality information. Bottom left corner insert shows an adult of Dactylochelifer latreillii; scale bar: 1 mm. (B) Map depicting sampling locations of the Dactylochelifer delimited lineages according to our preferred species delimitation outcome obtained via ASAP2 (upper right). Top left corner of the map shows the sampling locations in the whole sampling area and the upper right corner insert shows locations from Central Europe in detail.

Figure 3. Dactylochelifer species delimitation. (A) Phylogenetic tree of the genus Dactylochelifer with the results of different species delimitation approaches (right). Topology was obtained via the BI conducted in BEAST, the numbers on the nodes denote PP values, and the nodes with PP values > 0.95 are marked with black circles. Terminal tree taxa are color-coded according to the country of their collection. Specimen DL2042 (marked with asterisk) lacks locality information. Bottom left corner insert shows an adult of Dactylochelifer latreillii; scale bar: 1 mm. (B) Map depicting sampling locations of the Dactylochelifer delimited lineages according to our preferred species delimitation outcome obtained via ASAP2 (upper right). Top left corner of the map shows the sampling locations in the whole sampling area and the upper right corner insert shows locations from Central Europe in detail.

Figure 4. Hysterochelifer species delimitation. (A) Phylogenetic tree of the genus Hysterochelifer with the results of different species delimitation approaches (right). Topology was obtained via the BI conducted in BEAST, the numbers on the nodes denote PP values, and the nodes with PP values > 0.95 are marked with black circles. Terminal tree taxa are color-coded according to the country of their collection. Bottom left corner insert shows an adult female of H. tuberculatus; scale bar: 1 mm. (B) Map depicting sampling the locations of Hysterochelifer delimited lineages according to our preferred species delimitation outcome obtained via ASAP1 (upper right).

Figure 4. Hysterochelifer species delimitation. (A) Phylogenetic tree of the genus Hysterochelifer with the results of different species delimitation approaches (right). Topology was obtained via the BI conducted in BEAST, the numbers on the nodes denote PP values, and the nodes with PP values > 0.95 are marked with black circles. Terminal tree taxa are color-coded according to the country of their collection. Bottom left corner insert shows an adult female of H. tuberculatus; scale bar: 1 mm. (B) Map depicting sampling the locations of Hysterochelifer delimited lineages according to our preferred species delimitation outcome obtained via ASAP1 (upper right).

3.2.4. Rhacochelifer

The “Rhaco_out” dataset comprised 26 unique haplotypes derived from 36 Rhacochelifer individuals. The results of our species delimitation analyses ranged from three to six delimited lineages, depending on the method used (Figure 5). The most conservative results (i.e., three lineages) were obtained from the second ranked ASAP2 scheme (threshold distance 0.128), and the top ranked scheme (threshold distance 0.073) delimited four lineages (ASAP1), which was congruent with the results of the PTP approach. Statistical parsimony (95% parsimony criterion) performed in TCS delimited six lineages, while the GMYC method detected five lineages within our dataset.

Two distinct distribution patterns were detected within the Rhacochelifer diversity. Three unidentified independent Rhacochelifer lineages (A, B, C) showed tendencies toward local endemism, while the most sample-rich lineage D (R. maculatus morphotype) had a wider distribution that overlapped with the remaining geographically isolated lineages. Within R. maculatus, the longest geographical distance recorded among individuals with a shared haplotype was 1512 km (Table S3).

4. Discussion

4.1. Phylogenetic Relationships among Cheliferidae

Although Cheliferidae have been included in phylogenetic analyses of the order Pseudoscorpiones, their intrafamiliar relationships remain poorly understood because of the relatively sparse taxon sampling [43,44]. In our study, we analyzed a dataset comprising a significant portion of known European Cheliferidae diversity (66% of genera in Europe and 24% worldwide). Our Cox1 data did not possess sufficient informativeness to resolve deeper nodes of Cheliferidae relationships. However, the relationships that were supported contradicted the traditional views of Cheliferidae classification [34], i.e., the basal division into Philomaoriinae (comprising only Philomaria) and Cheliferinae (comprising the remaining genera) [67] subfamilies. Given the ancient origins of pseudoscorpions and the deep divergences between their lineages [44], it is clear that genomic or transcriptomic data will be needed in order to resolve the pseudoscorpion relationships at the intrafamilial level, as in other arachnid groups [89,90]. On the other hand, all sampled genera, with the exception of Mesochelifer, were recovered as monophyletic in our analyses (at least by BI; Figure 1), implying that traditional morphological characters used in higher level taxonomy [34,36] are likely only suitable for correct delimitation to the genus level in this family.

4.2. Cheliferidae Diversity and Delimitation

Pseudoscorpions belong to morphologically uniform arthropod taxa [18], and the application of molecular methods has highlighted that morphology-based delimitation approaches largely underestimate species diversity. The presence of cryptic diversity has been found in many families [21,25], and our results suggest that Cheliferidae are no exception to these reports. We identified several independent lineages in each of the studied genera, but the delimitation outcomes differed among the used methods (Figure 2, Figure 3, Figure 4 and Figure 5). Different delimitation methods are known to give different results, particularly in relation to how prone they are to oversplitting [38,91,92]. The delimitation should, thus, be based on an integration of multiple independent methods and the congruence between their outcomes [92,93]. In our case, it was impossible to integrate molecular approaches with morphologic or morphometric analyses [16,94]. We had limited sampling of males, which is the sex bearing species-specific characters; therefore, the assessment of character variability was not possible. For this reason, we cannot discard the notion that some of the independent lineages delimited in our analyses may correspond to previously described or synonymized species [17].

As pseudoscorpions have an accelerated substitution rate [95] and single-locus species delimitation methods tend to over-split taxa [38], we adopted a conservative approach while interpreting the species delimitation results. We, thus, favored the delimitation outcomes proposing a smaller number of putative species (ASAP threshold 0.062 to 0.072). Under the conservative approach, the divergences between our delimited lineages reached up to 20% (Table S4), which, on average, corresponds to interspecific distances recovered in other arachnid groups, e.g., spiders (17.4%) and harvestmen (19.4%) [96], as well as pseudoscorpions (19.7–25.6%) [16,21]. Similarly, intraspecific variability (1–4% in Chelifer, 1–3% in Dactylochelifer, 1–3% in Rhacochelifer and 0–4% in Hysterochelifer; Table S5) was congruent with the values reported from other arachnids e.g., spiders (0.7%), harvestmen (1.3%) and other pseudoscorpions (0.9–8%) [16,21,40,41,97].

4.2.1. Chelifer cancroides

The monotypic genus Chelifer, comprising cosmopolitan species C. cancroides, belongs to the most recognizable and frequently studied pseudoscorpion species [98,99,100,101]. It is ecologically plastic. Chelifer cancroides is most commonly found inside buildings, less frequently found in bee hives, under bark and inside tree hollows and nests and occasionally found in leaf litter (Table S1). According to our results, C. cancroides harbors three cryptic lineages that are not geographically segregated. All lineages co-occur in Central Europe; however, the most geographically widespread lineage (lineage C) is present across both Europe and North America (Figure 2). We even detected haplotype sharing among individuals from Canada and Czech Republic (up to 6837 km apart; Table S3). The lack of geographic structure, thus, supports our previous knowledge of C. cancroides ecology [27,100]. Its cosmopolitan distribution was traditionally attributed to a combined effect of phoretic dispersal and synanthropy [18,27,100,102]. Although phoresy in C. cancroides may be efficient, the haplotype sharing between continents may be more plausibly explained by human-mediated activities than long-distance trans-oceanic dispersal [103].

There are currently ten species synonymized to C. cancroides [17] as a consequence of previous morphological revisions of the genus (e.g., [100,104,105]). We detected three independent lineages in our analyses, thus contradicting the taxonomic status of C. cancroides. Unfortunately, most of the original descriptions of the formerly synonymized species are rather uninformative and lack the necessary details (e.g., [106]) to match with our delineated lineages. Furthermore, the challenges affiliated with analyzing historical pseudoscorpion materials [16] compound the complexity of the taxonomic revision of this species. Future studies utilizing C. cancroides as a model organism [101] should, thus, take into account the species most likely comprises cryptic diversity. Working with unrecognized species complexes may bias the conclusions drawn from the data, which unfortunately was the case in previously used animal models, such as ascidians [107], amphipods [108] and annelids [109].

4.2.2. Dactylochelifer

The genus Dactylochelifer with 52 described species represents the most specious genus of the family Cheliferidae [17]. Most of the described species have, with a few exceptions, geographically restricted occurrences [17]. Dactylochelifer latreillii, with its distribution spanning across Europe, North Africa and Central Asia, represents such an exception. The results of our study, based on samples covering a large portion of the species range, uncovered three independent lineages with overlapping (in some cases syntopic) distributions (lineages F, G, H) (Figure 3). These lineages neither can be distinguished based on the analysis of the traditional morphological characters nor exhibit strong habitat preferences. The representatives of a single lineage can be found in leaf litter, under tree bark or in cigar-like galls induced on the common reed (Phragmites australis) by the fruit fly genus Lipara (Table S1). The galls provide habitat for a wide assemblage of arthropods, including hymenopterans, hemipterans and beetles [110,111,112], that, in turn, could serve as phoretic carriers [27]. Compared to C. cancroides, D. latreillii does not have tendencies for synanthropy and likely owes its larger distributional patterns to effective phoretic dispersal. Shared haplotypes were detected in all of its lineages across substantial geographic distances (F: 1443 km, G: 90 km, H: 880 km; Figure 3, Table S3). The pattern is consistent with haplotype sharing distances in Chernes hahnii and Lamprochernes [13,16], which has largely been attributed to phoresy. However, some portions of D. latreillii distribution were not sampled in our analyses (e.g., North Africa and Mediterranean islands). Given the challenging taxonomy of Dactylochelifer, it is, thus, possible that the wide distribution [17] of the D. latreillii morphotype could also stem from the incorrect identification and underestimation of local species diversity.

Besides D. latreillii, we detected five additional lineages within the genus Dactylochelifer (lineages A–E, Figure 3) that likely correspond to previously described species with more geographically restricted ranges. However, the diagnostic characters are mainly male specific, and their variability is poorly understood. The combination of these factors prevented us from performing specific species identification.

4.2.3. Hysterochelifer

The genus Hysterochelifer comprises 13 nominal species distributed across North America, Europe and Central Asia; however, its center of diversification lies in the Mediterranean region [17]. Our sampling covers most of its known Mediterranean distribution (Figure 4). Our results recovered two morphologically distinct groups, one (lineage A) corresponding to H. meridianus morphotype and a second to several lineages identified as a H. tuberculatus morphotype (lineages B–I, Figure 4). The results of our analyses support the traditional view concerning the wide distribution of H. meridianus [17]. Although some delimitation approaches yielded between four and six independent lineages, the most conservative ASAP results support the existence of a single species (Figure 4). Its wide distribution, including the islands, and the fact that the specimens were exclusively collected under tree bark indicate that the species most likely disperses via phoresy. We detected shared haplotypes among localities 509 km apart (Table S3), but overall, 3% p-distance within the lineage A (H. meridianus morphotype) corresponds to intraspecific variability in other pseudoscorpion groups [25,40] and other species in this study (Table S5).

Interestingly, the lineages forming the H. tuberculatus morphotype likely have different ecological strategies. According to the literature, H. tuberculatus inhabits both leaf litter and tree bark [35]. However, it is clear from our results that this morphotype comprises multiple lineages, which may differ in terms of their ecology. All of our samples stem from leaf litter. It is, thus, interesting that at least some lineages (e.g., E, F, I) show wider distributions. For example, the lineage F is spanning between Italy and Turkey, and both lineages F and E are also present on the Mediterranean islands of Malta and Sicily. Such distributions would be consistent with phoretic dispersal [13,16], which would suggest that this manner of dispersal may be common in some soil-dwelling lineages. The remaining H. tuberculatus morphotype lineages (B, C, D, G, H) show restricted distributions that would most likely correspond to a more sedentary lifestyle. In this case, the lack of dispersal likely promoted the diversification of the H. tuberculatus morphotype. The inability to cross geographical barriers and maintain isolated populations is known to result in genetic diversification in many sedentary organisms. Similar patterns were detected among a wide variety of arachnids, such as harvestmen [113], scorpions [114] and pseudoscorpion families that do not possess phoretic dispersal capabilities [23,48].

4.2.4. Rhacochelifer

The distribution of Rhacochelifer is restricted to the Old World (Europe, Africa, and Central Asia), with the largest portion of the diversity (70% species) originating from the Mediterranean [17]. Our sampling is sparser relative to previously discussed genera (see above); however, our results suggest similar ecological and distribution patterns as those in Hysterochelifer (Figure 4 and Figure 5). Similar to H. meridianus, the R. maculatus morphotype (lineage D) spans most of the Mediterranean, with individuals sharing the same haplotype across 1575 km. Such a pattern also implies phoretic dispersal [13,16,45], which was historically documented in the literature [27]. Our samples of the R. maculatus morphotype were collected under tree bark, where pseudoscorpions typically encounter a wide variety of potential carriers. At the same time, we detected few soil-dwelling lineages (A–C) with restricted distributions (Figure 5), which implies their dispersal capabilities are limited. Despite our intensive and standardized sampling efforts across localities and different types of habitats inhabited by pseudoscorpions (Figure 6, Table S1), we were not able to collect a large number of locally distributed species, which may suggest that R. maculatus are significantly more abundant than other species.

4.3. Effects of Dispersal Strategies on Distribution and Diversity of Cheliferids

Our study targeted several related pseudoscorpion taxa with large areas of occurrence [17]; however, the results do not show geographically congruent distribution patterns. Interestingly, we detected varying degrees of haplotype sharing and population structuring among the target genera, which gave us insight into their ecology [16,21] and dispersal capabilities. For example, R. maculatus and C. cancroides, known phoretic species [27], show a similar lack of geographic structure compared to other species that are known to phoretically disperse, e.g., Dinocheirus arizonensis [45] and Lamprochernes spp. [16] (both Chernetidae). The same lack of geographic structure and haplotype sharing was detected in the H. meridianus morphotype (Hysterochelifer lineage A), H. tuberculatus (lineages E, F, I) and the D. latreillii morphotype (Dactylochelifer lineages F-H) (Figure 3 and Figure 4), which are not known for this type of dispersal. Documenting phoresy in pseudoscorpions is challenging due to their small size and cryptic lifestyle. Insights into their genetic background could, thus, be used as a proxy for studying their dispersal capabilities, supplementing our lack of direct observations [13]. The lack of geographic structure is even more pronounced in C. cancroides, which likely owes its distribution to the compound effect of synanthropy and phoretic dispersal [100]. At the same time, it seems that varying degrees of genetic structuring may exist between closely related taxa, e.g., the Hysterochelifer lineages A, E, F and I with wide distributions vs. the remaining Hysterochelifer lineages (B, C, D, G and H), that are geographical constrained (Figure 4). Such discrepancies may result from an association with rare and geographically constrained habitats [45,115] rather than different dispersal capabilities. Our findings suggest that the lifestyles of Cheliferidae are more flexible than the two previously followed strategies suggest [19,20,21]. We found support for a strong propensity of phoretic dispersal in lineages inhabiting temporary habitats; however, lineages from presumably stable environments (leaf litter and soil) might not be as limited dispersal-wise as would be expected.

Sedentary pseudoscorpions are well known for their abundant cryptic diversity (e.g., [40,41,42]); however, there is an increasing amount of evidence that the same phenomenon is realized between species with phoretic dispersal that were traditionally perceived as being widely distributed [13,16,39,115]. Due to the challenging taxonomy of cheliferids (see above) and their often cryptic lifestyles, we may still lack relatively basic knowledge concerning their diversity, distribution and ecology. Barcoding data can, thus, be used as a first step in research that will subsequently guide targeted taxonomical efforts and help to identify an appropriate methodological framework for detecting better diagnostic characteristics to define stable taxonomic units [16].

5. Conclusions

We investigated the genetic diversity of four widely distributed and relatively morphologically uniform Cheliferidae genera (Chelifer, Dactylochelifer, Rhacochelifer and Hysterochelifer) in order to (i) uncover the potential presence of cryptic diversity and (ii) evaluate whether the habitat preferences (i.e., ecological strategies) of species go hand in hand with their dispersal capabilities. We uncovered unidentifiable lineages within all genera, some of which might represent either cryptic diversity or previously synonymized species. Some lineages showed a complete lack of geographical structure and shared haplotypes over large distances, and some did not. A lack of geographical structure was detected among all known phoretic species, as well as among lineages that are not known to use this manner of dispersal. We, therefore, argue that genetic structure can be used as a proxy to evaluate the species’ dispersal manner and efficacy. Lineages from temporary habitats are most likely phoretic; however, lineages inhabiting presumably stable environments might be more capable in terms of dispersal than previously thought. Our results highlight the challenges involved in Cheliferidae research, particularly the poor definition of characters used for identification and the presence of male specimens needed for correct species determination. It is clear that all of the analyzed genera included in this paper are in need of thorough and integrative taxonomic revision.