# Relevance of ddRADseq method for species and population delimitation of closely related and widely distributed wolf spiders (Araneae, Lycosidae)

Vladislav Ivanov, Yuri Marusik, Julien Pétillon, Marko Mutanen

► To cite this version:

Vladislav Ivanov, Yuri Marusik, Julien Pétillon, Marko Mutanen. Relevance of ddRADseq method for species and population delimitation of closely related and widely distributed wolf spiders (Araneae, Lycosidae). Scientific Reports, 2021, 11 (1), pp.2177. 10.1038/s41598-021-81788-2. hal-03128210

HAL Id: hal-03128210

https://hal.science/hal-03128210

Submitted on 25 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN

# Relevance of ddRADseq method for species and population delimitation of closely related and widely distributed wolf spiders (Araneae, Lycosidae)

Vladislav Ivanov<fnref n="1" />✉, Yuri Marusik²,³, Julien Pétillon⁴ & Marko Mutanen<fnref n="1" />

Although species delimitation is often controversial, emerging DNA-based and classical morphology-based methods are rarely compared using large-scale samplings, even less in the case of widely distributed species that have distant, allopatric populations. In the current study, we examined species boundaries within two wolf spider species of the genus *Pardosa* (Araneae, Lycosidae), *P. riparia* and *P. palustris*. Wolf spiders constitute an excellent model for testing the relevance of traditional vs. modern methods in species and population delimitation because several closely related species are distributed over cross-continental geographic ranges. Allopatric populations of the two *Pardosa* species were sampled across Europe to Far East Russia (latitudinal range > 150°) and several dozen individuals were studied using morphological characters (morphometry of three measures for both sexes, plus five in males only and two in females only), DNA barcoding (COI sequencing) and double-digest restriction site associated DNA sequencing (ddRADseq). The results obtained allow for changing the taxonomic status of two Far East Russian populations to subspecies and ddRADseq proved to be a powerful tool for taxonomic research despite scarce sampling and inherent subjectivity of species delimitation in allopatry. Overall, this study pleads for both multi-criteria and more population-based studies in taxonomy.

Species delimitation is often challenging and controversial. Species is a central taxonomic category in various fields of biological research, yet the reality of species is still debated<fnref n="1" />⁻³. One of the reasons is the application of different operational criteria or species concepts results in various numbers of species under identical circumstances⁴. Combining different types of data in an integrative approach has the advantage of providing more informed decisions⁵⁻⁸. However, conflicts between results obtained through different species delimitation methods impose significant challenges for reaching uncontroversial conclusions³,⁷.

Traditionally, comparative morphology has served as the backbone in taxonomic routine⁹. Currently, broad datasets of DNA barcodes (portion of mitochondrial cytochrome c oxidase I in animals, COI) are often available along with morphological information. Studies implementing DNA barcodes and morphology on large geographic scales have shown correlations between COI divergence, morphological differences and geographical distance¹⁰. Operational species such as Barcode Index Numbers (BINs)¹¹,¹² provide a useful proxy for species boundaries. Yet, they do not always reflect them accurately as there are groups with high intraspecific or low interspecific variation¹²⁻¹⁷. The lumping of morphologically distinct spider species is reported from large German and Canadian datasets¹²,¹⁸. However, cases of shared mitochondrial DNA (mtDNA) between spider species have been known for a long time and are often explained by introgression or recent radiation¹⁹⁻²³. As mtDNA can be identical between closely related species regardless of isolation or extent of geographical distance, further genetic data is required to elucidate true boundaries of species.

Reduced-representation genome methods such as restriction site associated DNA sequencing (RADseq), transcriptomes or genotyping-by-sequencing provide an unprecedented opportunity to “peek” into the genome²⁴ and

<footnote n="1">¹Department of Ecology and Genetics, University of Oulu, Oulu, Finland. ²Institute for Biological Problems of the North, RAS, Magadan, Russia. ³Department of Zoology and Entomology, University of the Free State, Bloemfontein 9300, South Africa. ⁴UMR CNRS ECOBIO, Université de Rennes 1, Rennes, France. ⁵email: vladislav.ivanov@oulu.fi</footnote>
estimate divergence between populations, test for gene flow as well as compare closely related species. RADseq and similar methods have been successfully applied to species delimitation in different groups of organisms<sup>25–30</sup> demonstrating a promising step toward the standardization of taxonomic routine. Despite the abundance of molecular data challenges have been reported<sup>31,32</sup>.

Geographical sampling plays a crucial role in species delimitation. Across a species range, morphological and/or genetic variations may be observed in different populations, thus further presenting challenges in determining the taxonomic status of a particular population<sup>31-34</sup>. Incomplete and biased sampling can easily lead to incorrect assignment of sampled populations. However, the reality of taxonomic routine is that there is often no possibility to obtain representatives from key locations although there is a demand to assign collected diversity to new or existing species<sup>35</sup>. DNA-based methods are prone to incorrect species inference in the absence of intensive sampling<sup>36</sup>, but they have the advantage of providing insights into historical relationships between target populations, estimating gene flow and placing examined specimens in a context of closely related taxa<sup>34</sup>. Biological processes, such as horizontal gene transfer, introgression and incomplete lineage sorting, further complicate species inference<sup>37</sup>.

The target taxa of the study are representatives of the family Lycosidae, wolf spiders. These spiders constitute an excellent model for testing the relevance of traditional versus modern methods in species and population delimitation as several closely related species are distributed over cross-continental geographic ranges<sup>38</sup>. Recent advances in spider molecular systematics<sup>39,40</sup> suggest the family is relatively young with high to moderate diversification rates. Specialists in Lycosidae have taken advantage of molecular methods to study species boundaries and population genetics patterns in different genera but most studies were based on single or a few genes<sup>22,41-53</sup>. Species delimitation in Lycosidae is complicated by morphological homogeneity. However, courtship behavior<sup>54-56</sup>, ecological differentiation or geographical distribution<sup>57</sup> may be a more reliable source of evidence for taxonomic decisions. A previous study<sup>44</sup> suggested a clear association between morphospecies and genetic clades in sympatry when using double-digest restriction site associated DNA sequencing (ddRADseq) despite mitonuclear discordance.

In this study, the main interest was to evaluate the taxonomic status of allopatric populations in two separate cases of wolf spiders (Lycosidae) of genus *Pardosa* based on ddRADseq, DNA barcodes and morphology. Genus *Pardosa* is notorious for exhibiting shared mtDNA haplotypes in distinct morphological, ecological and behavioral species<sup>12,18,44,58-62</sup> which complicates the use of DNA barcodes for specimen identification and species delimitation. The focus was on the following two wolf spider species, *Pardosa riparia* (C. L. Koch, 1833) and *P. palustris* (Linnaeus, 1758). Both species are represented in the current study by two distant populations from Finland and the Magadan region in Far East Russia (FER). The first species, *P. riparia*, belongs to the *pullata* species group and inhabits Europe, Turkey, Russia, Central Asia and Japan. The second, *P. palustris*, belongs to the *monticola* species group and has an even wider distribution, occurring in Eurasia and North America but not in Japan<sup>63</sup>. Initially, it was believed that the distribution of both species is uninterrupted as indicated by data in the World Spider Catalog. However, after an intense survey of the literature and consultation with colleagues in Russia, it is certain that populations of both species in the Magadan region are truly allopatric in relation to the rest of Russia and Europe (Fig. S1 and associated references, Supplementary material (SM) 1). Geographical barriers in this case are vast mountain ranges on the border of the Eurasian and North American tectonic plates (e.g., Verkhoyansky, Chersky, Kolymsky ranges). Spiders are capable of dispersal over thousands kilometers by ballooning<sup>64-68</sup> but high mountains can impose a significant challenge for constant gene flow between populations. Thus, divergence between populations is likely to be driven by genetic drift as it is typical for the allopatric speciation process<sup>69</sup>. The initial hypothesis is that in both cases there are isolated populations of the same species which will be tested with molecular and morphological evidence. The current study suggests that ddRADseq could be superior or similar in its utility for species delimitation compared to classical morphological methods, especially in cases of genetically isolated populations where morphological evidence is of limited use.

## Results

**Morphological analysis.** No considerable qualitative differences in the structure of the copulatory organs were observed between the populations of *P. riparia* nor between the populations of *P. palustris*. Figures S2 and S3 in SM1 show the copulatory organs and habitus of the studied species. Some variations were detected, the first being the angle of the entrance duct openings relative to the septum in female *P. riparia* suggesting prezygotic isolation due to copulatory organs structural incompatibility (Fig. S2, SM1). The second variation was the terminal apophysis of *P. palustris* where in Finnish males it had long and sharp teeth on the distal margin, while in the Russian population it was smoother (Fig. S3, SM1). The PCA did not suggest differences between populations based on the measurements of morphological characters, except in the case of *P. riparia* females (Fig. 1). Boxplots (Fig. S5 and S6, SM1) indicated considerable overlap in copulatory organ measurements, while somatic characters showed greater divergence. Leg length (PatTib) showed the largest differences between populations in both species but failed the Bartlett test for homoscedasticity and could not be utilized for further analysis. The ANCOVA analysis suggested significant differences in the variances of male copulatory organ lengths in both species when taking into account body size. In females the size differences in copulatory organs were mostly insignificant. Additional t-test with unequal variances mostly corresponded with the ANCOVA (Table S6).

**ddRADseq assembly and datasets.** Many datasets with various parsing were assembled to examine the effect of the amount of missing data and to investigate how clustering at different values of *m* affected the *p*-distance results (Table 1). The amount of missing data ranged from 15.6 to 82.3%. Percentage of SNPs was between 4.3 and 9.6 for different datasets. Rip\_pal datasets included all sequenced specimens with outgroups, Pal\_mon datasets (except _m10c85) included available specimens from the *monticola* group, Rip datasets included only
<figure><img src="image_3.png" /><figcaption>**Figure 1.** PCA of morphological measurements data for *P. riparia* and *P. palustris*. The proportion of the overall variance explained by each PC is shown in parentheses. Areas bordered by the line represent 95% confidence intervals.</figcaption></figure>

<table><caption>**Table 1.** Summary of the ddRADseq datasets used in the downstream analysis. n: number of individuals; VAR: percentage of variable sites; PIS: percentage of parsimony informative sites. Dataset names are organized as following: Rip_pal—*pullata* and *monticola* groups combined plus outgroup species; Rip_pal_IN—same as previous but no outgroup; Rip—*pullata* group only; Pal_mon—*monticola* group only; m—minimum taxon coverage; c—clustering threshold.</caption><thead><tr><th>Matrix</th><th>n</th><th>Number of loci</th><th>Number of unlinked SNPs</th><th>Consensus sequences (bp)</th><th>VAR (%)</th><th>PIS (%)</th><th>Missing data (%)</th></tr></thead><tbody><tr><td>Rip_pal_m3c85</td><td>36</td><td>64,629</td><td>57,072</td><td>11,844,992</td><td>5.4</td><td>1.5</td><td>82.3</td></tr><tr><td>Rip_pal_m6c85</td><td>36</td><td>23,616</td><td>22,974</td><td>4,350,581</td><td>7.4</td><td>2.9</td><td>75</td></tr><tr><td>Rip_pal_m9c85</td><td>36</td><td>8950</td><td>8906</td><td>1,660,667</td><td>8.8</td><td>3.7</td><td>68.3</td></tr><tr><td>Rip_pal_m12c85</td><td>36</td><td>3146</td><td>3129</td><td>588,319</td><td>9.4</td><td>4.1</td><td>59.5</td></tr><tr><td>Rip_pal_m15c85</td><td>36</td><td>1262</td><td>1256</td><td>236,744</td><td>9.6</td><td>4.3</td><td>51.2</td></tr><tr><td>Rip_pal_m18c85</td><td>36</td><td>568</td><td>567</td><td>106,350</td><td>9.6</td><td>4.1</td><td>44.7</td></tr><tr><td>Rip_pal_m22c85</td><td>36</td><td>124</td><td>124</td><td>23,168</td><td>9.0</td><td>3.9</td><td>35.8</td></tr><tr><td>Rip_Pal_IN_m6c85</td><td>32</td><td>23,372</td><td>22,675</td><td>4,306,371</td><td>6.8</td><td>2.9</td><td>72.8</td></tr><tr><td>Rip_m10c85</td><td>18</td><td>684</td><td>670</td><td>127,926</td><td>6.6</td><td>2.6</td><td>40.6</td></tr><tr><td>Pal_mon_m10c85</td><td>12</td><td>2123</td><td>2071</td><td>395,237</td><td>5.4</td><td>2.7</td><td>15.6</td></tr><tr><td>Rip_m6c85</td><td>18</td><td>5431</td><td>5340</td><td>1,016,811</td><td>6.4</td><td>2.4</td><td>58.4</td></tr><tr><td>Rip_m4c85</td><td>18</td><td>11,914</td><td>11,609</td><td>2,232,854</td><td>5.7</td><td>1.8</td><td>67.9</td></tr><tr><td>Pal_mon_m4c85</td><td>14</td><td>39,831</td><td>36,284</td><td>7,266,287</td><td>4.8</td><td>1.8</td><td>52.9</td></tr><tr><td>Pal_mon_m6c85</td><td>14</td><td>22,115</td><td>21,213</td><td>4,058,676</td><td>5.5</td><td>2.3</td><td>43.9</td></tr><tr><td>Pal_mon_m3c85</td><td>14</td><td>54,495</td><td>46,703</td><td>9,929,008</td><td>4.3</td><td>1.3</td><td>57.7</td></tr></tbody></table>
the *pullata* group while the Pal_mon_m10c85 dataset contained only *P. palustris* and *P. ovtchinnikovi*. The datasets Rip_m10c85 and Pal_mon_m10c85 were used for STRUCTURE analysis to measure the extent of gene flow between populations and between closely related species. High *m* values were used in these analyses to diminish the proportion of missing data. A value of 85 was found to be optimal for the clustering threshold and was used for all datasets. Summary statistics for the assemblies are provided in SM2.

**Phylogenetic inference, haplotype networks and species trees.** Maximum likelihood trees for COI that included only specimens attempted for ddRADseq confirmed the previously reported DNA barcode sharing within the *pullata* group (Fig. S7, SM1). The same was observed in the ML COI tree that included all public records in BOLD (SM3). In the case of *P. palustris*, the pattern was similar, populations from Russia were clustered separately to the Finnish and Faroe Island specimens for the ddRADseq mimicking dataset (Fig. S7, SM1). The same was observed in the ML tree built in IQ-Tree and species ID tree built in BOLD for the *monticola* group (SM3 and SM4). However, when building a species ID tree in BOLD for the larger dataset of the *monticola* group, FER and American specimens fell within one of the European clusters (SM5). When the same was done for the *pullata* group, *P. riparia* becomes a sister group to all other species in the *pullata* group (SM4).

The haplotype network built for the *pullata* group (Fig. 2a) indicated that all sequenced *P. riparia* from FER belong to the same haplotype, while Finnish specimens had several haplotypes mixed with other representatives of the species group. In the case of *P. palustris* (Fig. 2b), its position was similar to *P. pullata* in its own species group, i.e. it did not share haplotypes with other representatives of the *monticola* group. Specimens from FER and America nevertheless could be separated from the *P. palustris* haplotypes from other countries, thus indicating divergence in COI.

The ML tree based on ddRADseq data distinguished between Russian and Finnish samples for both species but also recognized morphospecies of *pullata* and *monticola* groups as separate lineages (Fig. 3). Species trees calculated using SVDquartets implemented in PAUP* showed high support for splitting distant populations into two different species. However, the bootstrap support in inner nodes remained quite low (50), thus SVDquartets did not support *P. pullata* as a distinctly separate species. The same was observed within the *monticola* group, where the split between *P. palustris* and the rest of the species had a bootstrap value of 50 (Fig. S8, SM1).

**Species delimitation.** BIN analysis conducted in BOLD in both cases assigned Russian specimens to already known European species (*P. riparia* BIN: BOLD:AAF7515; *P. palustris* BIN: BOLD:ACH2627). Species discovery methods GMYC and bPTP returned contradictory results based on COI analysis and were largely dependent on the dataset analyzed. GMYC oversplit both species when analyzed alone or in combination with species from the same species groups. When the two species groups were combined, the FER populations were lumped into one species in both cases. Results of bPTP were similar in that the number of suggested species decreases when more distant lineages were added but they were largely oversplit. Summary of the GMYC and bPTP analyses are in Tables S7 and S8 (SM1), figures of the trees are available upon request.

The Bayes factor calculated on the SNAPP results favors the two species hypothesis in both cases (381.01 for *P. riparia* and 295.22 for *P. palustris*).

**Testing for admixture.** The STRUCTURE results suggested that *P. riparia* from Finnish and Russian populations belongs to the same gene pool in relation to other species of the *pullata* group. The highest probability was observed at K5 (Fig. S9, SM1). Population genetic analysis suggested that representatives of *P. riparia* from FER and Finland belong to the same species (Fig. 3). The situation with *P. palustris* is different. STRUCTURE results suggested that the Russian and Finnish populations had limited shared ancestry (Fig. 3). The mean estimated probability was higher for K3 than for K2 (Fig. S9, SM1), thus *P. palustris* populations appeared more diverged than *P. riparia* populations.

**Pairwise distance patterns in ddRADseq datasets.** Intra- and interspecific *p*-distances varied dramatically depending on *m* value. The clear difference between intra- and interspecific distances was observed from *m* = 3 to *m* = 9 in the *pullata* group and from *m* = 3 to *m* = 15 in the *monticola* group (Fig. 4). Concurrently, *p*-distances between the FER and Finnish populations for both species remained within the 95% confidence interval of intraspecific distance irrespective of *m*. In addition, means of *p*-distances between *P. palustris* populations remained non-overlapping while *p*-distances between *P. riparia* populations were closer.

The *p*-distances were visualized similarly to the barcode gap, i.e., as plotted frequencies of *p*-distances for interspecific, target populations and intraspecific values. There was evident gap between intra- and interspecific distances and target populations are likely to belong to intraspecific values distribution at *m* = 6 (Fig. 5).

**Discussion**

While species delimitation in spiders using molecular methods has a long tradition (e.g.⁷⁰), there are only a few studies that incorporate both morphology and genetics in Lycosidae and even less so using multiple genes (e.g.⁴¹,⁴³,⁵⁶,⁵⁹,⁷¹). Therefore, this study is one of the first combining ddRADseq, DNA barcodes and morphology to study the problem of allopatric species in wolf spiders and, to our best knowledge, the first to compare FER populations of target species to European populations. While results of different analyses are not entirely conclusive, similar trends are observed for both *P. riparia* and *P. palustris*. It is acknowledged that limited sampling imposes difficulties in drawing completely unbiased conclusions. However, we suggest that the abundance of data collected should elucidate the questions asked. The congruence of the results are discussed below, their biological meaning, taxonomic decisions that can be made based on the presented evidence and implications of this research for further studies.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** COI haplotype networks. (a) *pullata* species group, (b) *monticola* species group. Each circle represents a haplotype and circle size is proportional to strain frequency. Different colors represent different species. Lines between haplotypes and short solid lines are single mutational steps. Dotted lines outline haplotypes by populations. Black dots are predicted or missing haplotypes.</figcaption></figure>

All types of data (morphological, COI and ddRADseq) clearly show divergence between distant populations. It is to be expected for the case of true allopatric populations where interaction between populations is prohibited
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Maximum likelihood tree (RipPal_m6c85 dataset) and STRUCTURE barplots (Rip_m10c85 and Pal_mon_m10c85 datasets) based on ddRADseq data. Numbers indicate bootstrap support values for nodes. STRUCTURE barplots represent the partition with highest K supported where K = 5 for *pullata* group and K = 3 for *monticola* group.</figcaption></figure>

by a natural barrier. Therefore, one line of evidence that supports the two-species hypothesis for both cases is geographical isolation that lasted long enough to promote detectable divergence in morphological characters and DNA. Distribution maps provided (Fig. S1, SM1) are the most up to date for both species.

The traditional morphological approach that was followed in this study was aiming to find characters that would be indicative of divergence to the extent of prezygotic reproductive isolation and/or considerable divergence between populations. In entomology and arachnology, genital differences have played a particular role in often showing rapid divergent evolution⁷². Differences in genital structures have been thought to form a mechanical reproductive barrier between species through character displacement (the lock-and-key hypothesis), but genital evolution is more likely driven by sexual selection⁷³. Qualitative assessment of morphological characters indeed shows several indicative features that could be used for reliable identification including the copulatory
Further intensive sampling is required to confirm if this character could be used as a diagnostic feature. PCA based on the same measurements suggests greater divergence in *P. riparia* than in *P. palustris*. Variation reported here coupled with the geographical distribution could be enough to erect FER populations to at least subspecies status. However, we believe that morphology alone is insufficient in this case as subjectivity of character choice is an inherent problem for traditional taxonomic studies. DNA sequences were thus examined for congruence with the morphological results, as the DNA sequence analyses rely on more general evolutionary theories and follow more standardized procedures for acquiring and analyzing data.

DNA barcoding again proved to be limited source of evidence for species delimitation in *Pardosa*. Nevertheless, both FER populations have diverged haplotypes from European ones, especially in *P. riparia*. The eastern population of *P. palustris* is clearly closer to Alaskan populations than to European (one mutation difference against five) while *P. riparia* populations form separate clusters in both ML trees and haplotype networks. Species delimitation methods such as GMYC and bPTP are heavily criticized for oversplitting species in spiders<sup>33,75</sup> and despite the results being reported here we refrained including them in the final taxonomic conclusion. Briefly, in the context of closely related species both populations in both species are split into three to 41 species by GMYC or bPTP. If two species groups are combined, then the *P. palustris* populations appear as one homogeneous cluster irrespective of geography. While *P. riparia* populations are lumped with other species of the *pullata* group much like the species ID trees from BOLD and the COI ML trees. Overall, COI remains inconclusive for *Pardosa* species delimitation, but slight differences observed corresponded to the distribution data, i.e., isolated FER populations have different haplotypes from European populations, and it can be used as proxy for quick identification using DNA barcodes.

Genomic data provided evidence that FER populations in both species belong to a separate cluster based on ML trees and the species delimitation results from SNAPP for both populations in both species. It should be noted here that multispecies coalescent methods (MSC) have been criticized for delimiting structure rather than species<sup>76</sup>, while it is suggested that they are better justified from an evolutionary point of view<sup>77</sup>. Therefore, given the limited number of scenarios that could be tested with SNAPP and the lack of agreement about the utility of MSC for species delimitation, the supported two species hypotheses should be treated as a reflection of divergence between distant populations rather than decisive evidence for species boundaries. The species tree computed in PAUP* has high support for the split but low support for separation of the target species from their respective species groups. Thus confirming the presence of genetic differences between populations but telling little about actual species statuses. Simultaneously, p-distances and population genetic analysis suggest that distant populations are a single entity rather than two diverged species when compared to closely related species. Only in the case of STRUCTURE results for *P. palustris*, distant populations are forming two clusters with limited shared ancestry. However, this could be an artefact of insufficient sampling. Overall, ddRADseq data suggest detectable divergence between European and FER populations but do not unequivocally support the presence of separate species in the Magadan region.

The taxonomic decision about FER populations is straightforward given a certain level of congruence between DNA and morphological data. All lines of evidence despite some discrepancies suggest a level of divergence of FER populations of *P. riparia* and *P. palustris* detectable by morphological examination and DNA sequencing. However, not to the extent to confidently claim species level divergence if using closely related species as a yardstick. Therefore, subspecies status for eastern populations of *P. riparia* and *P. palustris* is the most appropriate and compliant to the data analyzed. The main point of distinguishing differences between populations in scientific name is that it shows detected diversity in a quick and straightforward way, easily accessible for various specialists thus enabling proper faunistic, taxonomic and evolutionary studies planning and provoking further research. There seems to be a tradition of avoiding subspecies status in arachnology when no clear morphospecies can be described<sup>78</sup>, thus we do not suggest names at this stage to avoid inflation. Nevertheless, this current research implies that described cases of *P. riparia* and *P. palustris* could serve as model populations for studying recent divergence and speciation as such. However, a more favorable outcome would be a clear taxonomic decision that we are planning to present when additional data is collected. A somewhat similar view is stated in a considerably more data rich study of the *Habronattus tarsalis* species complex (Salticidae). There, multiple isolated populations exhibited deep genetic and morphological differences but the authors were reluctant to discuss their taxonomic status as the process of divergence is believed to be more important for future studies and naming species is "somewhat missing the point"<sup>79</sup> While for evolutionary and population genetics studies such a view can be justified, it does not add to the documentation of biodiversity which is one of the major unfulfilled tasks for humanity. As genome scale sequencing is becoming more and more accessible it can strongly facilitate species discovery and validation, we would encourage researchers to be more bold in using DNA evidence alone for taxonomic decisions especially if genome level sequencing is available.

To conclude, *ddRADseq* is a powerful tool for research in taxonomy as it can supplement traditional morphological methods with historical perspective and provide data on current genetic processes. Therefore, genomic tools can give insight into the old problem of species delimitation in allopatry by adding additional layers of objectivity despite inherently arbitrary final decisions. However, at the current stage of the research we would like to restrain from giving formal description of subspecies before the major populations across the whole area of distribution are studied even though the amount of evidence is more than sufficient for such a conclusion.

## Materials and methods

**Sampling.** Specimens for morphological (155 specimens), COI (25 specimens) and *ddRADseq* (17 specimens) analyses were sampled from target populations of *P. riparia* and *P. palustris* in Finland and the Magadan Region, FER. These were coupled with representatives of *pullata* and *monticola* species groups from Europe and Asia, several other *Pardosa* species and outgroup species from Lycosinae subfamily [*Alopecosa aculeata*]
(Clerck, 1757) and *Trochosa spinipalpis* (F. O. P.-Cambridge, 1895)]. Wider sampling of species from target species groups was done by mining public records from the Barcode of Life Database (BOLD) for COI analysis (212 specimens, 14 species). Detailed information about all the specimens used in the study is in Tables S1, S2 and S3 (SM1). Specimens were collected by hand or using pit-fall traps and preserved in absolute or 70% ethanol. No experiments were performed on live animals hence ethical approval was not required. Protected species were not utilized in this study and no permits were needed to collect the specimens.

**Morphological examination.** Specimens were carefully examined to identify species, and the whole body and copulatory organs of the target species were photographed. Differences in structures of copulatory organs are traditionally used as the most reliable source of qualitative information for spider species delimitation¹². Characters that are reported to be informative for *Pardosa*⁶¹ were photographed and measured, with the following modifications: CarLen and CarWid—dorsal length and width of cephalothorax in males and females, PatTib—length of patella + tibia in males and females, CymbLen—length of male palpal tarsus, CymbTip—length of tip of male palpal tarsus, BulbLen and BulbWid—length and width of bulbus in males, ApophLen—length of tegular apophysis for *P. riparia* males, SeptLen and SeptWid—length and width of septum in females. Males (*P. riparia* 15 specimens from each population; *P. palustris* 26 from Finland, 22 from FER) and females (*P. riparia* 15 specimens from each population; *P. palustris* 21 from Finland, 22 from FER) were analyzed separately. The measurement scheme can be found in Fig. S4 (SM1) and the measurements in Tables S4 and S5 (SM1). Principle Component Analysis (PCA) was used to analyze variation in the morphological data using RStudio with *pcomp* command⁸⁰ and the results were visualized with *ggplot2* package⁸¹. Each character between two populations was checked for homoscedasticity in R (bartlett.test) and compared with ANCOVA using CarLen as a covariate to control for body size variation following the same rationale as in⁵⁷. In addition, a t-test with unequal variances was utilized as an additional comparison between measurements. To visualize overlap between measured characters boxplots were produced with the *ggplot2* package.

**Mitochondrial DNA sequencing.** The COI sequencing was performed in the Canadian Centre for DNA barcoding (CCDB). One leg from each specimen was placed into individual wells in a 96-well plate prefilled with ethanol and the plate was sent to CCDB. There, DNA was extracted and the COI gene (654 bp) was sequenced following standard protocol⁸². All collection data and specimens' images were uploaded into BOLD together with successful sequences by corresponding Sample ID (SM1). In addition, publicly available sequences were included. For the full list of samples see Table S3 (SM1).

**ddRADseq library preparation.** Genomic DNA (gDNA) was extracted with DNeasy Blood & Tissue Kit (Qiagen) according to the manufacturer's protocol. One to four legs provided the required amount of gDNA for subsequent whole genome amplification (WGA). The rationale for using WGA was that the specimens were quite small and for morphological analysis it was important to keep them as intact as possible. Furthermore, the amount of gDNA extracted from the legs would not suffice for ddRADseq library preparation. The REPLI-g Mini Kit (Qiagen) for WGA was used according to the manufacturer's instructions. It was assumed that WGA did not bias the results significantly as there was evidence of a negligible effect on the study outcome⁸³⁻⁸⁵ and examples of other studies successfully implementing WGA⁸⁶⁻⁸⁹.

Previously published protocols for the ddRADseq library preparation were followed⁹⁰,⁹¹. Briefly, a combination of PstI and MseI restriction enzymes were used for shearing the DNA. Ligated adapters designed for this pair of enzymes were ligated and the samples were purified with AMPure XP magnetic beads (Agencourt). Samples were pooled into six sub-pools based on concentration measurements using the PicoGreen Kit (Molecular Probes). Automated size-selection was performed using Blue Pippin (Sage Science, 1.5% agarose cartridge) to obtain a library with a mean of 300 bp long fragments. Selected fragments were amplified with the Phusion High-Fidelity PCR Master Mix (Finnzymes) and the products were purified with AMPure XP magnetic beads. Prior to sequencing, the size distribution and concentration of sub-pools and the final pool were checked using Bioanalyzer (Agilent). The final library was sequenced on an Illumina HiSeq 2500 machine, 100 PE in FIMM, (Institute for Molecular Medicine, Finland). Raw DNA reads from ddRADseq are available at the NCBI Sequence Read Archive (BioProject ID: PRJNA345307 and PRJNA595572).

**ddRADseq data bioinformatics.** Quality control was performed with FastQC⁹². Paired-end reads were assembled de novo using the ipyrad workflow (https://ipyrad.readthedocs.io/). The software provides strict filtering steps to ensure the absence of paralogs and low-quality base calls in the assembly. The most important parameter altered from the default was the minimum number of samples per locus (*m*). Test datasets were assembled with different combination of the parameter to find optimum values as well as to investigate the influence of missing data on consecutive analyses. Values of *m* from a minimum of three to a maximum of 22 were tested. Two types of datasets were assembled, the first included all species with an outgroup and the second contained only target species with their closest relatives. The first dataset was used for phylogenetic inference and *p*-distance computation and the second set for species delimitation and population genetic analyses. All strict filtering steps allowed for additional control of possible WGA bias. The proportion of missing data was calculated with Mesquite version 3.51⁹³.

**Phylogenetic analyses and haplotype networks.** Phylogenetic analyses were required to reveal historical relationships among taxa and to test the validity of the prevailing species hypotheses. Maximum likelihood (ML) trees for ddRADseq datasets were inferred with the IQ-Tree program⁹⁴ with 5000 ultrafast bootstrap support for branches⁹⁵. Haplotype networks for COI were computed with a modified script in RStudio (origi-
nally from https://johnbhorne.wordpress.com/2016/09/15/still-making-haplotype-networks-the-old-way-how-to-do-it-in-r/) using packages *ape*/5.3⁹⁶ and *pegas*/0.11⁹⁷ and PopArt/1.7 using the TCS method⁹⁸,⁹⁹. Maximum likelihood trees for COI datasets were constructed using MEGA7¹⁰⁰ with bootstrap support estimated from 500 replicates and the GTR+G substitution model.

**Population genetics.** STRUCTURE version 2.3.1 was used to detect admixture between populations¹⁰¹. To optimize the runs, StrAuto program¹⁰² was used with the following parameters: 500,000 replicates with 50,000 burn-in, 10 replicates for each K (K = 1–6 for the *pullata* group and K = 1–4 for the *monticola* group). An optimal K was estimated with STRUCTURE HARVESTER¹⁰³ based on the ad hoc ΔK statistics¹⁰⁴. Replicates were permuted in the program CLUMPP¹⁰⁵ and bar plots were visualized with the program Distruct¹⁰⁶.

**Species delimitation.** For COI datasets, species discovery methods were utilized⁷ to assess the overall efficiency of DNA barcodes for species delimitation purposes in the focal species. Firstly, the BIN system was used for the initial assignment of specimens to operational taxonomic units (OTUs) as implemented in BOLD¹¹. The species ID trees for all available *P. palustris* and *P. riparia* are presented in the SM3.

Secondly, a browser version of bPTP¹⁰⁷ was utilized with the following parameters: unrooted tree, 100,000 Markov chain Monte Carlo (MCMC) generation, 0.1 burn-in. Input trees used were COI ML trees built with IQ-Tree for five datasets: *P. riparia* only (n = 21), *P. palustris* only (n = 45), *pullata* group (n = 117), *monticola* group (95 specimens) and combined *pullata* and *monticola* groups (n = 212).

The third species discovery method used was a browser version of GMYC¹⁰⁸. The same datasets employed for bPTP were utilized. Ultrametric trees required as GMYC input were computed in BEAST 2.5¹⁰⁹ with Strict Clock Model, 10 million MCMC chain length, 100,000 pre-burn. Priors: coalescent constant population model, uniform frequency parameter, exponential gamma rate, population size 1/X, uniform rates from 1 to 100. Single threshold was used for the GMYC run¹¹⁰.

For species delimitation with the ddRADseq data, Bayes-factor species delimitation (BFD*) was implemented in SNAPP, BEAST 2.60³⁴. Specimens were subsampled from the initial ddRADseq assemblies to form two datasets. The first dataset included 18 specimens from the *pullata* group with 543 unlinked SNPs (40.6% missing data) while the second included ten specimens from the *monticola* group with 1153 unlinked SNPs, (1% missing data). Dataset selection was based on a minimum missing data parameter. Two alternative scenarios were tested for both datasets, where representatives of *P. riparia* or *P. palustris* from different populations belong to the same or different species. Bayes factor was calculated to evaluate the alternative hypotheses. Run parameters and specimens included in each dataset can be found in SM1.

Species trees for ddRADseq datasets were estimated in a coalescent framework using all data with the program SVDquartets¹¹¹ implemented in PAUP*¹¹². Sub-sampling was done at 10,000 quartets. Sampled quartets were assembled into a species tree using a variant of Quartet FM¹¹³, which was the implementation recommended by the developers of SVDquartets.

**Pairwise distances.** The estimation of *p*-distances was performed using PAUP* for the ddRADseq datasets. The goal of the analysis was to determine if the levels of variation in ddRADseq corresponded to the different species hypotheses derived from the COI data, morphology and ddRADseq data analyses. The *p*-distances for ddRADseq data were calculated for the dataset that included all species and an outgroup. Measured *p*-distances were grouped into intraspecific, interspecific and *p*-distances between compared geographically distant populations. Mean, quartiles and 95% confidence intervals were calculated and absolute values of *p*-distances were visualized as boxplots using *ggplot2* package in RStudio for ddRADseq datasets with different *m* values. In addition, frequencies of *p*-distance values within species, between target populations and between distinct species were computed and plotted using *ggplot2* package in RStudio.

**Data availability**

Raw ddRADseq reads can be found in NCBI Sequence Read Archive (BioProject ID: PRJNA345307 and PRJNA595572).

Received: 14 December 2019; Accepted: 11 January 2021

Published online: 26 January 2021

**References**

1. Lee, M. S. Y. Species concepts and species reality: salvaging a linnaean rank. *J. Evol. Biol.* **16**, 179–188 (2003).

2. Seifert, B. A pragmatic species concept applicable to all eukaryotic organisms independent from their mode of reproduction or evolutionary history. *Soil Org.* **86**, 85–93 (2014).

3. Reydon, T. A. C. & Kunz, W. Species as natural entities, instrumental units and ranked taxa: new perspectives on the grouping and ranking problems. *Biol. J. Linn. Soc.* **126**, 623–636 (2019).

4. De Queiroz, K. Species concepts and species delimitation. *Syst. Biol.* **56**, 879–886 (2007).

5. Dejaco, T., Gassner, M., Arthofer, W., Schlick-Steiner, B. C. & Steiner, F. M. Taxonomist's nightmare ... evolutionist's delight: an integrative approach resolves species limits in jumping bristletails despite widespread hybridization and parthenogenesis. *Syst. Biol.* **65**, 947–974 (2016).

6. Padial, J. M., Miralles, A., De la Riva, I. & Vences, M. The integrative future of taxonomy. *Front. Zool.* **7**, 16 (2010).

7. Carstens, B. C., Pelletier, T. A., Reid, N. M. & Satler, J. D. How to fail at species delimitation. *Mol. Ecol.* **22**, 4369–4383 (2013).

8. Gratton, P. et al. Testing classical species properties with contemporary data: how 'bad species' in the brassy ringlets (*Erebia tyndarus* complex, Lepidoptera) turned good. *Syst. Biol.* **65**, 292–303 (2016).

9. Fujita, M. K., Leaché, A. D., Burbrink, F. T., McGuire, J. A. & Moritz, C. Coalescent-based species delimitation in an integrative taxonomy. *Trends Ecol. Evol.* **27**, 480–488 (2012).
10. Mutanen, M. et al. Allopatry as a Gordian knot for taxonomists: patterns of DNA barcode divergence in Arctic-Alpine Lepidoptera. *PLoS ONE* 7, e47214 (2012).

11. Ratnasingham, S. & Hebert, P. D. N. A DNA-based registry for all animal species: the Barcode Index Number (BIN) system. *PLoS ONE* 8, e66213 (2013).

12. Astrin, J. J. et al. Towards a DNA barcode reference database for spiders and harvestmen of Germany. *PLoS ONE* 11, e1–24 (2016).

13. Morinière, J. et al. Barcoding Fauna Bavaria: 78% of the Neuropterida Fauna Barcoded! *PLoS ONE* 9, e109719 (2014).

14. Pentinsaari, M., Hebert, P. D. N. & Mutanen, M. Barcoding beetles: a regional survey of 1872 species reveals high identification success and unusually deep interspecific divergences. *PLoS ONE* 9, e108651 (2014).

15. Hausmann, A. et al. Genetic patterns in European geometrid moths revealed by the Barcode Index Number (BIN) system. *PLoS ONE* 8, e84518 (2013).

16. Hausmann, A., Haszprunar, G. & Hebert, P. D. N. DNA barcoding of the geometrid fauna of Bavaria (Lepidoptera): successes, surprises, and questions. *PLoS ONE* 6, 1–9 (2011).

17. Pohjoismäki, J. L. O., Kahanpää, J. & Mutanen, M. DNA barcodes for the northern European tachinid flies (Diptera: Tachinidae). *PLoS ONE* 11, e1–23 (2016).

18. Blagoev, G. A. et al. Untangling taxonomy: a DNA barcode reference library for Canadian spiders. *Mol. Ecol. Resour.* **16**, 325–341 (2016).

19. Johannesen, J. & Veith, M. Population history of *Eresus cinnaberinus* (Araneae: Eresidae) colour variants at a putative species transition. *Heredity* (Edinb). **87**, 114–124 (2001).

20. Croucher, P. J. P., Oxford, G. S. & Searle, J. B. Mitochondrial differentiation, introgression and phylogeny of species in the *Tege-naria atrica* group (Araneae: Agelenidae). *Biol. J. Linn. Soc.* **81**, 79–89 (2004).

21. Lattimore, V. L., Vink, C. J., Paterson, A. M. & Cruickshank, R. H. Unidirectional introgression within the genus *Dolomedes* (Araneae:Pisauridae) in southern New Zealand. *Invertebr. Syst.* **25**, 70 (2011).

22. De Busschere, C., Van Belleghem, S. M. & Hendrickx, F. Inter and intra island introgression in a wolf spider radiation from the Galápagos, and its implications for parallel evolution. *Mol. Phylogenet. Evol.* **84**, 73–84 (2015).

23. Spasojevic, T., Kropf, C., Nentwig, W. & Lasut, L. Combining morphology, DNA sequences, and morphometrics: revising closely related species in the orb-weaving spider genus *Araniella* (Araneae, Araneidae). *Zootaxa* **4111**, 448 (2016).

24. Andrews, K. R., Good, J. M., Miller, M. R., Luikart, G. & Hohenlohe, P. A. Harnessing the power of RADseq for ecological and evolutionary genomics. *Nat. Rev. Genet.* **17**, 81–92 (2016).

25. Herrera, S. & Shank, T. M. RAD sequencing enables unprecedented phylogenetic resolution and objective species delimitation in recalcitrant divergent taxa. *Mol. Phylogenet. Evol.* **100**, 70–79 (2016).

26. Dierickx, E. G., Shultz, A. J., Sato, F., Hiraoka, T. & Edwards, S. V. Morphological and genomic comparisons of Hawaiian and Japanese Black-footed Albatrosses (*Phoebastria nigripes*) using double digest RADseq: Implications for conservation. *Evol. Appl.* **8**, 662–678 (2015).

27. Hou, Y. et al. Thousands of RAD-seq loci fully resolve the phylogeny of the highly disjunct arctic-alpine genus *Diapensia* (Diapensiaceae). *PLoS ONE* **10**, e0140175 (2015).

28. Eaton, D. A. R. & Ree, R. H. Inferring phylogeny and introgression using RADseq data: an example from flowering plants (Pedicularis: Orobanchaceae). *Syst. Biol.* **62**, 689–706 (2013).

29. Razkin, O. et al. Species limits, interspecific hybridization and phylogeny in the cryptic land snail complex *Pyramidula*: the power of RADseq data. *Mol. Phylogenet. Evol.* **101**, 267–278 (2016).

30. Kozlov, M. V., Mutanen, M., Lee, K. M. & Huemer, P. Cryptic diversity in the long-horn moth *Nemophora degeerella* (Lepidoptera: Adelidae) revealed by morphology, DNA barcodes and genome-wide ddRAD-seq data. *Syst. Entomol.* **42**, 329–346 (2017).

31. Dincă, V., Lee, K. M., Vila, R. & Mutanen, M. The conundrum of species delimitation: a genomic perspective on a mitogenetically super-variable butterfly. *Proc. R. Soc. B Biol. Sci.* **286**, 20191311 (2019).

32. Linck, E. et al. Dense geographic and genomic sampling reveals paraphyly and a cryptic lineage in a classic sibling species complex. *Syst. Biol.* **68**, 956–966 (2019).

33. Hamilton, C. A., Hendrixson, B. E., Brewer, M. S. & Bond, J. E. An evaluation of sampling effects on multiple DNA barcoding methods leads to an integrative approach for delimiting species: a case study of the North American tarantula genus *Aphonopelma* (Araneae, Mygalomorphae, Theraphosidae). *Mol. Phylogenet. Evol.* **71**, 79–93 (2014).

34. Leaché, A. D., Fujita, M. K., Minin, V. N. & Bouckaert, R. R. Species delimitation using genome-wide SNP data. *Syst. Biol.* **63**, 534–542 (2014).

35. Avendaño, J. E., Arbeláez-Cortés, E. & Cadena, C. D. On the importance of geographic and taxonomic sampling in phylogeography: a reevaluation of diversification and species limits in a Neotropical thrush (Aves, Turdidae). *Mol. Phylogenet. Evol.* **111**, 87–97 (2017).

36. Puechmaille, S. J. The program STRUCUTRE does not reliably recover the correct population structure when sampling is uneven: subsampling and new estimators alleviate the problem. *Mol. Ecol. Resour.* **16**, 608–627 (2016).

37. Mallet, J., Besansky, N. & Hahn, M. W. How reticulated are species? *BioEssays* **38**, 140–149 (2016).

38. Marusik, Y. M. & Koponen, S. A survey of spiders (Araneae) with Holarctic distribution. *J. Arachnol.* **33**, 300–305 (2005).

39. Garrison, N. L. et al. Spider phylogenomics: untangling the spider tree of Life. *PeerJ* **4**, e1719 (2016).

40. Fernández, R. et al. Phylogenomics, diversification dynamics, and comparative transcriptomics across the spider tree of life. *Curr. Biol.* **28**, 1489–1497.e5 (2018).

41. Vink, C. J. & Mitchell, A. D. 12S DNA sequence data confirms the separation of *Alopecosa barbipes* and *Alopecosa accentuata* (Araneae, Lycosidae). *Bull. Br. Arachnol. Soc.* **12**, 242–244 (2002).

42. Vink, C. J., Mitchell, D. A., Paterson, A. M., Mitchell, A. D. & Paterson, A. M. A preliminary molecular analysis of phylogenetic relationships of Australasian wolf spider genre (Araneae, Lycosidae). *J. Arachnol.* **30**, 227–237 (2002).

43. Murphy, N. P. et al. Phylogenetic reconstruction of the wolf spiders (Araneae: Lycosidae) using sequences from the 12S rRNA, 28S rRNA, and NADH1 genes: Implications for classification, biogeography, and the evolution of web building behavior. *Mol. Phylogenet. Evol.* **38**, 583–602 (2006).

44. Chang, J., Song, D. & Zhou, K. Incongruous nuclear and mitochondrial phylogeographic patterns in two sympatric lineages of the wolf spider *Pardosa astrigera* (Araneae: Lycosidae) from China. *Mol. Phylogenet. Evol.* **42**, 104–121 (2007).

45. Hebets, E. A. & Vink, C. J. Experience leads to preference: experienced females prefer brush-legged males in a population of syntopic wolf spiders. *Behav. Ecol.* **18**, 1010–1020 (2007).

46. Park, Y. C., Yoo, J.-S., Schwarz, M. P., Murphy, N. & Kim, J.-P. Molecular phylogeny of East Asian Wolf Spiders (Araneae: Lycosidae) inferred from mitochondrial 12S ribosomal DNA. *Ann. Entomol. Soc. Am.* **100**, 1–8 (2007).

47. Gotch, T. B., Adams, M., Murphy, N. P. & Austin, A. D. A molecular systematic overview of wolf spiders associated with Great Artesian Basin springs in South Australia: evolutionary affinities and an assessment of metapopulation structure in two species. *Invertebr. Syst.* **22**, 151 (2008).

48. Correa-Ramírez, M. M., Jiménez, M. L. & García-de León, F. J. Testing species boundaries in *Pardosa sierra* (Araneae: Lycosidae) using female morphology and COI mtDNA. *J. Arachnol.* **38**, 538–554 (2010).

49. Planas, E., Fernández-Montraveta, C. & Ribera, C. Molecular systematics of the wolf spider genus *Lycosa* (Araneae: Lycosidae) in the Western Mediterranean Basin. *Mol. Phylogenet. Evol.* **67**, 414–428 (2013).
50. González-Trujillo, R. et al. Pleistocene refugia and their effects on the phylogeography and genetic structure of the wolf spider *Pardosa sierra* (Araneae: Lycosidae) on the Baja California Peninsula. *J. Arachnol.* **44**, 367–379 (2016).

51. Piacentini, L. N. et al. A revision of the wolf spider genus *Diapontia* Keyserling, and the relationships of the subfamily Sosippinae (Araneae: Lycosidae). *Arthropod Syst. Phylogeny* **75**, 387–415 (2017).

52. Ivanov, V., Lee, K. M. & Mutanen, M. Mitonuclear discordance in wolf spiders: genomic evidence for species integrity and introgression. *Mol. Ecol.* **27**, 1681–1695 (2018).

53. Piacentini, L. N. & Ramírez, M. J. Hunting the wolf: a molecular phylogeny of the wolf spiders (Araneae, Lycosidae). *Mol. Phylogenet. Evol.* **136**, 227–240 (2019).

54. Töpfer-Hofmann, G., Cordes, D. & Helversen, O. V. Cryptic species and behavioural isolation in the *Pardosa lugubris* group (Araneae, Lycosidae), with description of two new species. *Bull. Br. Arachnol. Soc.* **11**, 257–274 (2000).

55. Chiarle, A., Kronestedt, T. & Isaia, M. Courtship behavior in European species of the genus *Pardosa* (Araneae, Lycosidae). *J. Arachnol.* **41**, 108–125 (2013).

56. Just, P., Opatova, V. & Dolejš, P. Does reproductive behaviour reflect phylogenetic relationships? An example from Central European *Alopecosa* Wolf spiders (Araneae: Lycosidae). *Zool. J. Linn. Soc.* **185**, 1039–1056 (2019).

57. Puzin, C., Leroy, B. & Pétillon, J. Intra- and inter-specific variation in size and habitus of two sibling spider species (Araneae: Lycosidae): taxonomic and biogeographic insights from sampling across Europe. *Biol. J. Linn. Soc.* **113**, 85–96 (2014).

58. Zehethofer, K. & Sturmbauer, C. Phylogenetic relationships of Central European wolf spiders (Araneae: lycosidae) inferred from 12S ribosomal DNA sequences. *Mol. Phylogenet. Evol* **10**, 391–398 (1998).

59. Muster, C. & Berendonk, T. U. Divergence and diversity: Lessons from an arctic-alpine distribution (*Pardosa saltuaria* group, Lycosidae). *Mol. Ecol.* **15**, 2921–2933 (2006).

60. Muster, C., Maddison, W. P., Uhlmann, S., Berendonk, T. U. & Vogler, A. P. Arctic-alpine distributions—metapopulations on a continental scale? *Am. Nat.* **173**, 313–326 (2009).

61. Sim, K. A., Buddle, C. M. & Wheeler, T. A. Species boundaries of *Pardosa concinna* and *P. lapponica* (Araneae: Lycosidae) in the northern Nearctic: morphology and DNA barcodes. *Zootaxa* **3884**, 169–178 (2014).

62. Slowik, J. & Sikes, D. S. Molecular systematics of the *Pardosa groenlandica* species complex (Araneae: Lycosidae): large sample sizes fail to find monophyletic species. *Can. Entomol.* **147**, 643–664 (2015).

63. World Spider Catalog. World Spider Catalog. Natural History Museum Bern. Version 21.0 http://wsc.nmbe.ch (2020) doi:10.24436/2.

64. Duffey, E. Aerial dispersal in spiders. In *Proc. 17th Eur. Colloq. Arachnol.* 187–191 (1998).

65. Bell, J. R., Bohan, D. A., Shaw, E. M. & Weyman, G. S. Ballooning dispersal using silk: world fauna, phylogenies, genetics and models. *Bull. Entomol. Res.* **95**, 69–114 (2005).

66. Krehenwinkel, H. & Tautz, D. Northern range expansion of European populations of the wasp spider *Argiope bruennichi* is associated with global warming-correlated genetic admixture and population-specific temperature adaptations. *Mol. Ecol.* **22**, 2232–2248 (2013).

67. Postiglioni, R., Eisenberg, A., Carlozzi, A. & Bidegaray-Batista, L. The dark side of ballooning: nocturnal aerial dispersal in wolf spiders from the South American Coastline. *Arachnology* **17**, 312–316 (2017).

68. Cho, M., Neubauer, P., Fahrenson, C. & Rechenberg, I. An observational study of ballooning in large spiders: nanoscale multifibers enable large spiders' soaring flight. *PLOS Biol.* **16**, e2004405 (2018).

69. Gavrilets, S. Perspective: models of speciation: what have we learned in 40 years?. *Evolution* (N. Y.) **57**, 2197–2215 (2003).

70. Ramirez, M. G. & Beckwitt, R. D. Phylogeny and historical biogeography of the spider genus *Lutica* (Araneae, Zodariidae). *J. Arachnol.* **23**, 177–193 (1995).

71. Planas, E. & Ribera, C. Description of six new species of *Loxosceles* (Araneae: Sicariidae) endemic to the Canary Islands and the utility of DNA barcoding for their fast and accurate identification. *Zool. J. Linn. Soc.* **174**, 47–73 (2015).

72. Eberhard, W. G. *Sexual Selection and Animal Genitalia* (Harvard University Press, Cambridge, 1985). https://doi.org/10.4159/ harvard.9780674330702.

73. Arnqvist, G. Comparative evidence for the evolution of genitalia by sexual selection. *Nature* **393**, 784–786 (1998).

74. Ballarin, F., Marusik, Y. M., Omelko, M. M. & Koponen, S. On the *Pardosa monticola* species-group (Araneae: Lycosidae) from Middle Asia. *Arthropoda Sel.* **21**, 161–182 (2012).

75. Lopardo, L. & Uhl, G. Testing mitochondrial marker efficacy for DNA barcoding in spiders: a test case using the dwarf spider genus *Oedothorax* (Araneae: Linyphiidae: Erigoninae). *Invertebr. Syst.* **28**, 501 (2014).

76. Sukumaran, J. & Knowles, L. L. Multispecies coalescent delimits structure, not species. *Proc. Natl. Acad. Sci.* **114**, 1607–1612 (2017).

77. Leaché, A. D., Zhu, T., Rannala, B. & Yang, Z. The spectre of too many species. *Syst. Biol.* **68**, 168–181 (2019).

78. Kraus, O. Why no subspecies in spiders? In *Proceedings of the 19th European Colloquium of Arachnology* (eds. Toft, S. & Scharff, N.) 303–314 (Aarhus University Press, 2002).

79. Hedin, M., Foldi, S. & Rajah-Boyer, B. Evolutionary divergences mirror Pleistocene paleodrainages in a rapidly-evolving complex of oasis-dwelling jumping spiders (Salticidae, *Habronattus tarsalis*). *Mol. Phylogenet. Evol.* **144**, 106696 (2020).

80. RStudio Team. RStudio: Integrated Development for R. RStudio (2015).

81. Wickham, H. *ggplot2: Elegant Graphics for Data Analysis* (Springer, Berlin, 2016).

82. deWaard, J. R., Ivanova, N. V., Hajibabaei, M. & Hebert, P. D. N. Assembling DNA barcodes. In *Environmental Genomics* (eds Martin, C. C. & Martin, C. C.) 275–294 (Humana Press, Totowa, 2008). https://doi.org/10.1007/978-1-59745-548-0_15.

83. Barker, D. L. et al. Two methods of whole-genome amplification enable accurate genotyping across a 2320-SNP linkage panel. *Genome Res.* **14**, 901–907 (2004).

84. Pinard, R. et al. Assessment of whole genome amplification-induced bias through high-throughput, massively parallel whole genome sequencing. *BMC Genom.* **7**, 216 (2006).

85. Han, T. et al. Characterization of whole genome amplified (WGA) DNA for use in genotyping assay development. *BMC Genom.* **13**, 217 (2012).

86. Rheindt, F. E., Fujita, M. K., Wilton, P. R. & Edwards, S. V. Introgression and phenotypic assimilation in zimmerius flycatchers (Tyrranidae): population genetic and phylogenetic inferences from genome-wide SNPs. *Syst. Biol.* **63**, 134–152 (2014).

87. Blair, C., Campbell, C. R. & Yoder, A. D. Assessing the utility of whole genome amplified DNA for next-generation molecular ecology. *Mol. Ecol. Resour.* **15**, 1079–1090 (2015).

88. Burford Reiskind, M. O. et al. Development of a universal double-digest RAD sequencing approach for a group of nonmodel, ecologically and economically important insect and fish taxa. *Mol. Ecol. Resour.* **16**, 1303–1314 (2016).

89. Lee, K. M. et al. Information dropout patterns in restriction site associated DNA phylogenomics and a comparison with multilocus sanger data in a species-rich moth genus. *Syst. Biol.* **67**, 925–939 (2018).

90. Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S. & Hoekstra, H. E. Double digest RADseq: an inexpensive method for de novo SNP discovery and genotyping in model and non-model species. *PLoS ONE* 7, e37135 (2012).

91. DaCosta, J. M. & Sorenson, M. D. ddRAD-seq phylogenetics based on nucleotide, indel, and presence—absence polymorphisms: analyses of two avian genera with contrasting histories with contrasting histories. *Mol. Phylogenet. Evol.* **94**, 122–135 (2015).

92. Andrews, S. FastQC: a quality control tool for high throughput sequence data (2010).

93. Maddison, W. P. & Maddison, D. R. Mesquite: A Modular System for Evolutionary Analysis. Version 3.51 (2018).
94. Nguyen, L.-T., Schmidt, H. A., von Haeseler, A. & Minh, B. Q. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. *Mol. Biol. Evol.* **32**, 268–274 (2015).

95. Hoang, D. T., Chernomor, O., von Haeseler, A., Minh, B. Q. & Vinh, L. S. UFBoot2: improving the ultrafast bootstrap approximation. *Mol. Biol. Evol.* **35**, 518–522 (2018).

96. Paradis, E., Claude, J. & Strimmer, K. APE: analyses of phylogenetics and evolution in R language. *Bioinformatics* **20**, 289–290 (2004).

97. Paradis, E. pegas: an R package for population genetics with an integrated-modular approach. *Bioinformatics* **26**, 419–420 (2010).

98. Leigh, J. W. & Bryant, D. Popart: full-feature software for haplotype network construction. *Methods Ecol. Evol.* **6**, 1110–1116 (2015).

99. Clement, M., Snell, Q., Walke, P., Posada, D. & Crandall, K. TCS: estimating gene genealogies. In *Proceedings 16th International Parallel and Distributed Processing Symposium* (IEEE, 2002). https://doi.org/10.1109/IPDPS.2002.1016585.

100. Kumar, S., Stecher, G. & Tamura, K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. *Mol. Biol. Evol.* **33**, 1870–1874 (2016).

101. Pritchard, J. K., Stephens, M. & Donnelly, P. Inference of population structure using multilocus genotype data. *Genetics* **155**, 945–959 (2000).

102. Chhatre, V. E. & Emerson, K. J. StrAuto: automation and parallelization of STRUCTURE analysis. *BMC Bioinform.* **18**, 192 (2017).

103. Earl, D. A. & vonHoldt, B. M. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. *Conserv. Genet. Resour.* **4**, 359–361 (2012).

104. Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of individuals using the software structure: a simulation study. *Mol. Ecol.* **14**, 2611–2620 (2005).

105. Jakobsson, M. & Rosenberg, N. A. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. *Bioinformatics* **23**, 1801–1806 (2007).

106. Rosenberg, N. DISTRUCT: a program for the graphical display of population structure. *Mol. Ecol. Notes* **4**, 137–138 (2004).

107. Zhang, J., Kapli, P., Pavlidis, P. & Stamatakis, A. A general species delimitation method with applications to phylogenetic place-ments. *Bioinformatics* **29**, 2869–2876 (2013).

108. Fujisawa, T. & Barraclough, T. G. Delimiting species using single-locus data and the generalized mixed yule coalescent approach: a revised method and evaluation on simulated data sets. *Syst. Biol.* **62**, 707–724 (2013).

109. Bouckaert, R. et al. BEAST 25: an advanced software platform for Bayesian evolutionary analysis. *PLOS Comput. Biol.* **15**, e1006650 (2019).

110. Talavera, G., Dincă, V. & Vila, R. Factors affecting species delimitations with the GMYC model: insights from a butterfly survey. *Methods Ecol. Evol.* **4**, 1101–1110 (2013).

111. Chifman, J. & Kubatko, L. Quartet inference from SNP data under the coalescent model. *Bioinformatics* **30**, 3317–3324 (2014).

112. Swofford, D. L. PAUP*: *Phylogenetic Analysis Using Parsimony, version 4.0 b10* (Sinauer Associates, 2003).

113. Reaz, R., Bayzid, M. S. & Rahman, M. S. Accurate phylogenetic tree reconstruction from quartets: a heuristic approach. *PLoS ONE* **9**, e104008 (2014).

## Acknowledgements

We thank Kyung Min Lee for help with ddRADseq analysis and comments on the manuscript; Kirill Mikhailov, Sergei Esyunin and Galina Azarkina for unpublished data on distribution of *P. riparia* and *P. palustris*; Sepideh Shafaie for providing *P. pontica* specimens; Mikko Pentinsaari for help with haplotype networks; Helsinki University Biological Station (Kilpisjärvi) and Tvärminne Zoological Station for help with spider sampling in Finland; CSC – IT Center for Science for providing facilities for computation; staff at the Biodiversity Institute of Ontario for their continuous help in generating sequences, entering data into BOLD and aiding the curation of this information; FIMM staff for their high-quality work. Funding was provided by “Suomen Akatemia”, Grant No. 277984 accorded to MM, “Kvantum Institute” (accorded to MM) and “Suomen Hyönteistieteellinen Seura” (personal grant to VI). Finally, we thank the four anonymous reviewers for multiple comments and excellent suggestions for the earlier draft of this paper.

## Author contributions

V.I. and M.M. designed the research; V.I., Y.M. and J.P. did the field collections; V.I. and Y.M. did qualitative morphology check; V.I. did morphometry, prepared samples for DNA barcoding, did ddRADseq library preparation and all statistical and bioinformatics analysis; V.I., M.M., Y.M., J.P. drafted the manuscript. All authors discussed the results and contributed to the final manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-021-81788-2.

**Correspondence** and requests for materials should be addressed to V.I.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_6.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or
format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the
Creative Commons licence, and indicate if changes were made. The images or other third party material in this
article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the
material. If material is not included in the article's Creative Commons licence and your intended use is not
permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021