# Across mountains and ocean: species delimitation and historical connectivity in Holarctic and Arctic-Alpine wolf spiders (Lycosidae, Pardosa)

Vladislav Ivanov, Gergin Blagoev, Samuel Danflous, Peter Gajdos, Toke T Høye, Kyung Min Lee, Yuri Marusik, Cecilie Lohse Mielec, Christoph Muster, Julien Pétillon, et al.

► To cite this version:

Vladislav Ivanov, Gergin Blagoev, Samuel Danflous, Peter Gajdos, Toke T Høye, et al. Across mountains and ocean: species delimitation and historical connectivity in Holarctic and Arctic-Alpine wolf spiders (Lycosidae, Pardosa). Insect Systematics and Diversity, 2023, 7 (5), 10.1093/isd/ixad018. hal-04264643

HAL Id: hal-04264643

https://hal.science/hal-04264643v1

Submitted on 26 Mar 2025

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
Molecular Phylogenetics, Phylogenomics, and Phylogeography

# Across mountains and ocean: species delimitation and historical connectivity in Holarctic and Arctic-Alpine wolf spiders (Lycosidae, *Pardosa*)

Vladislav Ivanov¹,*,º, Gergin Blagoev²,º, Samuel Danflous³, Peter Gajdosš⁴, Toke Thomas Høye⁵, Kyung Min Lee¹,⁶,º, Yuri Marusik⁷,⁸,⁹,º, Cecilie Lohse Mielec⁵, Christoph Muster¹⁰, Julien Pétillon¹¹,¹², Jörg Spelda¹³,º, Marko Mutanen¹,º

¹Ecology and Genetics Research Unit, University of Oulu, Oulu, Finland,²Centre for Biodiversity Genomics, University of Guelph, Guelph, ON, Canada,³Muséum d'Histoire naturelle de Toulouse, Toulouse, France,⁴Slovak Academy of Sciences, Bratislava, Slovakia,⁵Department of Ecoscience, Aarhus University, Aarhus, Denmark,⁶Zoology Unit, Finnish Museum of Natural History, University of Helsinki, Helsinki, Finland,⁷Institute for Biological Problems of the North, RAS, Magadan, Russia,⁸Department of Zoology and Entomology, University of the Free State, Bloemfontein, South Africa,⁹Altai State University, Barnaul, Russia,¹⁰University of Greifswald, Greifswald, Germany,¹¹UMR CNRS ECOBIO, Université de Rennes 1, Rennes, France,¹²Institute for Coastal and Marine Research, Nelson Mandela University, Gqeberha, South Africa,¹³Bavarian Natural History Collections – Bavarian State Collection of Zoology (SNSB-ZSM), Munich, Germany. Corresponding author, mail: uladzislausw@gmail.com, vladislav.ivanov@oulu.fi

Subject Editor: Lauren Esposito

Received on 15 December 2021; revised on 20 June 2023; accepted on 30 August 2023

Holarctic species offer great opportunities to study biogeography, phylogenetics, taxonomy, and local adaptation. Species that are considered conspecific between the Palearctic and the Nearctic realms are often split into 2 or more species when denser sampling and molecular markers are applied. Similar in complexity but at a finer geographical scale are species groups having Arctic-Alpine distributions where lineages have complicated demographic histories due to glacial dynamics. In both cases, allopatric speciation might not result in fast differentiation of morphological characters if environmental conditions in isolated areas are similar and the main driver of variability is genetic drift. Here, we study the Holarctic *Pardosa hyperborea* (Thorell, 1872) and its closest European relatives to assess their taxonomic status and patterns of genetic variability. Based on DNA barcodes and genomic data from double-digest restriction site associated sequencing, we propose that the North American populations should be regarded as a distinct species (*P. luteola* Emerton, 1894, **stat. resurr.), possibly consisting of several independent lineages. With the help of D-statistics, population genetic simulations and phylogenetic networks analysis, we demonstrate historical introgression among European species of the group and a likely explanation for shared DNA barcodes among allopatric and fully differentiated species. Our study exposes a promising model for studying speciation processes and demographic history in parallel on both sides of the Atlantic Ocean and demonstrates the usefulness of genomic tools in elucidating the taxonomy and biogeography of taxa across broad geographic scales.

Key words: ddRAD sequencing, DNA barcoding, Lycosidae, introgression, glaciation

## Introduction

Widely distributed taxa are of great interest to systematists, phylogeographers, and population geneticists because they enable the study of historical and current processes of dispersal and adaptation of organisms as well as the history of speciation. While higher-level taxa (e.g., families or genera) usually have continent-wide or

even cosmopolitan distributions, at the species level such cases are rarer and usually remain poorly studied, especially in terrestrial fauna and flora (Procheş and Ramdhani 2013). The Holarctic realm poses a significant challenge for studying such species because of multiple disjunctions and connections between Europe, Asia, and the eastern and western parts of North and Central America during the
past 200 million yr (Enghoff 1995, Sanmartín et al. 2001). There are cases where a single species indeed has a Holarctic distribution either naturally or mediated by human activities (Landry et al. 2013). In other studies, a distinct split between Palearctic and Nearctic lineages is shown and new species are described (Kraus 2002, Johnsen et al. 2010, Schär et al. 2018, Divakar et al. 2019, Abel et al. 2020, Brunke et al. 2020, Gilligan et al. 2020).

Another example where species have a disrupted and patchy distribution caused by major planetary scale processes are species with Arctic-Alpine (boreo-montane) ranges. Arctic-Alpine distributions can be largely explained by recurring glaciation oscillations during the Quaternary that shaped the genetic diversity of current populations through repeated cycles of species range contraction and expansion, causing bottlenecks, secondary contacts, and allopatric divergence in discontinuous areas (Hewitt 2000, Schmitt 2007, 2009). Numerous studies in different groups of organisms—vascular plants (Hou et al. 2015, Sanz et al. 2017, György et al. 2018), molluscs (Fiorentino et al. 2016), insects (Macher et al. 2015, Schmitt et al. 2016), spiders (Bidegaray-Batista et al. 2014), mammals (Sawyer and Cook 2016)—suggest various scenarios of diversification, including retreats to refugia and consecutive expansions during interglacials. Currently, DNA sequencing data is often used to study such complicated processes. It allows researchers to fine tune the analysis of divergence time and populations/species interaction under complex demographic scenarios as well as to study peculiarities of local adaptation (Wagner et al. 2018, Dufresnes et al. 2020, 2021, Pan et al. 2020, Pittet et al. 2020). As in the case of Holarctic species, Arctic-Alpine species are a challenge for taxonomists due to the allopatric nature of speciation (Funk and Omland 2003, Mutanen et al. 2012, Dincă et al. 2019, Linck et al. 2019, Ivanov et al. 2021).

The wolf spiders of the *Pardosa saltuaria* species group have a remarkable distribution including both Holarctic and Arctic-Alpine species in Europe and Asia. Currently, there are 7 described morphospecies in the group (Wunderlich 1984, Marusik et al. 1996, Vogel 2004). Four of them are known only from Europe—*P. oreophila* Simon, 1937; *P. saltuaria* (L. Koch, 1870); *P. evelinae* Wunderlich, 1984; *P. drenskii* Buchar, 1968 (Nentwig et al. 2020)—1 is from Asia, *P. oksalai* Marusik, Hippa & Koponen, 1996 (Marusik et al. 1996, Khadzieva et al. 2019) and 1 is known only from United States and Mexico—*P. californica* Keyserling, 1887 (Vogel 2004). The most widespread species—*P. hyperborea* (Thorell, 1872) — inhabits Europe (including Iceland), Greenland, and North America but not the eastern Palearctic (Fig. 1, the references are in the Supplementary Material). The 5 species that occur in Europe (including *P. hyperborea*) seem to be differentiated morphologically (Wunderlich 1984), but 4 of them are nonmonophyletic at the mitochondrial DNA (mtDNA) level—*P. hyperborea*, *P. oreophila* (from the Alps), *P. evelinae* and *P. saltuaria*—while *P. drenskii* is clearly differentiated at both genetic and morphological levels (Muster and Berendonk 2006). Previous studies on this group in Europe were performed using 1 mitochondrial and 3 nuclear markers and suggested: (i) strong historical (or ongoing) gene flow between Scandinavian and central European species; (ii) the Pyrenees and Balkans to be genetically highly diverged lineages and the Pyrenean populations of *P. oreophila* are suggested to be a separate undescribed species; (iii) presence of admixture between morphospecies in the Alps; (iv) uncertain taxonomic status for the currently recognized morphospecies (Muster and Berendonk 2006, Muster et al. 2009). In the context of *Pardosa*, nonmonophyly of mtDNA in closely related species was described among other representatives of the genus (Astrin et al. 2016, Blagoev et al. 2016, Ivanov et al. 2018). Similar patterns are observed in other groups of organisms, e.g., birds (Hibbets et al.

2020), bats (Mao and Rossiter 2020), snails (Miura et al. 2020), crustaceans (Wang et al. 2020), millipedes (Spelda et al. 2011). Suggested causes behind mitonuclear discordance include incomplete lineage sorting (Obertegger et al. 2018, Hundsdoerfer and Kitching 2020), introgression (Bonnet et al. 2017, Leduc-Robert and Maddison 2018, Oxford and Bolzern 2018) or bacterial endosymbionts influencing mtDNA distribution patterns (Sahoo et al. 2018, Elias-Costa et al. 2019, Tseng et al. 2019). Based on our previous research, introgression seems to be the most plausible explanation for shared mtDNA between closely related sympatric *Pardosa* species (Ivanov et al. 2018).

In this study, we focus on the taxonomy of spiders from the *saltuaria* species group from Europe and North America as well as patterns of divergence in nuclear and mitochondrial DNA (nDNA and mtDNA respectively) under allopatric conditions. We use double-digest restriction site associated sequencing (ddRADseq) and DNA barcodes (mitochondrial COI) as the main source of evidence. We attempt to answer 3 main questions: (i) Do genomic data support the pattern of connectivity and admixture among species and populations known from previous studies? (ii) How much do genomic data and DNA barcodes contribute to species delineation in this species complex? (iii) Is *P. hyperborea* really a species of Holarctic distribution?

## Materials and Methods

### Samples and COI Sequences

We included 120 specimens for ddRADseq analysis. It included representatives of the *saltuaria* species group and an outgroup from the *monticola* species group. Specimens currently recognized as *P. hyperborea* were collected in Canada (n = 23), Greenland (n = 19) and Finland (n = 26). Other species of the group included *P. californica* from California, USA (n = 1), *P. oreophila* from the Pyrenees, France (n = 21) and *P. saltuaria* from the Tatra Mountains, Slovakia (n = 19). Specimens of *P. palustris* from Faroe Islands, Russia (Magadan region) and Finland (n = 11) were used as outgroups (Supplementary Table S1). All specimens were collected with pit-fall traps or by hand between 2007 and 2018, preserved in absolute or 70% ethanol and put in the freezer at -20 °C. No experiments were done on living individuals.

All specimens used for ddRADseq were sent for COI sequencing to the Canadian Centre for DNA Barcoding (CCDB). Metadata for samples and specimen images were uploaded to BOLD database. We used 419 specimens for DNA barcoding analysis, including newly generated barcodes and publicly available data. Out of 109 ingroup specimens included in ddRADseq library 99 were successfully sequenced for COI. A list of all specimens used for DNA barcode analysis (with BOLD IDs) is in Supplementary Table S2.

### ddRADseq Library

Genomic DNA (gDNA) was extracted with the DNeasy Blood & Tissue Kit (Qiagen) following standard the manufacturer's protocol. One to four legs were taken from each specimen to obtain gDNA, and the REPLI-g Mini Kit (Qiagen) was used for subsequent whole genome amplification (WGA).

We used the same library preparation protocol as in Ivanov et al. (2018) which is largely based on the initial ddRAD protocol (Peterson et al. 2012, DaCosta and Sorenson 2014). Briefly, we used *Pst*I and *Mse*I restriction enzymes that allowed for combining runs from previously sequenced libraries. The library was multiplexed by ligating adapters specific to the pair of restriction
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Distribution map of *P. hyperborea* based on current literature review. Filled squares—approximate locations of *P. hyperborea* records in Palearctic; empty square—dubious record of *P. hyperborea* in Siberia; filled circles—approximate locations of *P. hyperborea* records in Nearctic. The references on which the map is based are in **Supplementary Material**.</figcaption></figure>

enzymes (48 unique individual-specific barcodes and 2 indices that are compatible with a paired-end Illumina sequencing run). Samples were pooled into sub-pools based on DNA concentration estimated with the PicoGreen Kit (Molecular Probes). After ligation, size selection was performed with a BluePippin (Sage Science, 2% agarose cartridge, mean length 300 bp). The selected DNA fragments were then amplified with a Phusion High-Fidelity PCR Master Mix (Finnzymes). PCR products were purified with AMPure XP magnetic beads (Agencourt). The concentration and size distribution of the sub-pools were measured with Bioanalyzer (Agilent). Libraries were finally pooled into a single library in equimolar amounts and sequenced on an Illumina HiSeq 2500 machine, 100 PE in the Institute for Molecular Medicine, Finland (FIMM). Raw DNA reads from ddRADseq are available at the NCBI Sequence Read Archive (BioProject IDs PRJNA345307 and PRJNA595572).

## Ipyrad Assembly Parameters

The raw reads were demultiplexed to samples and assessed for quality using FastQC (Andrews 2010) at the FIMM. Before assembly forward and reverse reads were merged with PEAR 0.9.8 (Zhang et al. 2014). Then, ipyrad 0.9.30 (Eaton and Overcast 2020) was used for reference and *de novo* assemblies. The reference sequence was a set of *P. pseudoannulata* genome scaffolds (BioProject ID PRJNA512163). A clustering threshold of 0.85 was used for all *de novo* and reference assemblies. We started assemblies with all 120 specimens and consequently removed those that had a low number of assembled loci (less than 2% of the total number of recovered loci, see Tables 1–2 and Supplementary Tables S3–S4 for details). Different values of minimum taxon coverage (*m*) were used to explore the balance between the recovered loci and the proportion of missing data. Reference and *de novo* datasets with the same number of samples included the same specimens.
## Genetic Structure and Admixture

To explore DNA barcode variation, a haplotype network was computed using 419 specimens of species of the saltuaria group available in BOLD. It was done in RStudio with packages ape/5.3, pegas/0.11, and PopArt/1.7 using the TCS method (Clement et al. 2002, Paradis et al. 2004, Paradis 2010, Leigh and Bryant 2015).

To understand the geographical patterns of the genetic structure, we performed STRUCTURE v.2.3.1 (Pritchard et al. 2000) analysis using SNP frequency data. In this analysis, 2,538 putative SNPs were identified by randomly selecting a single SNP from each locus from the ipyrad output (De_novo_91 dataset, Tables 1 and 2). Ten replicates were run at each value of *K* between 1 and 7, each consisting of 50k burn-in generation and 500k postburn-in generations. We used StrAuto to automate the STRUCTURE processing of samples (Chhatre and Emerson 2017). Replicates were permuted using CLAMPP (Jakobsson and Rosenberg 2007) and the optimal *K* was inferred using Structure Harvester (Earl and vonHoldt 2012) according to the *ad hoc* Δ*K* statistics (Evanno et al. 2005), which is the second-order rate of change of the likelihood function. The STRUCTURE results were visualized using DISTRUCT (Rosenberg 2004).

To visualize the genetic structure and admixture patterns spatially, we mapped the membership coefficient matrix obtained from STRUCTURE with the highest Δ*K* against the matrix of geographical coordinates. The distribution map was generated through SimpleMappr online web service (www.simplemappr.net) and modified in Adobe Illustrator CS6 (Adobe Inc.).

## Phylogenetic Inference

IQ-Tree program version 1.6.12 (Nguyen et al. 2015) was used to build maximum likelihood (ML) trees for all ddRADseq datasets using the model supported by ModelFinder based on BIC (variations of TVM for all datasets). The program was run with 5,000 ultrafast bootstrap (Hoang et al. 2018), -bnni correction to reduce the risk of overestimating branch support due to severe model violations and -alrt 5000 to perform additional assessment of branch support with single branch test (Guindon et al. 2010). One tree was also built

for the COI dataset using the GTR model and otherwise identical parameters described above.

## Species Delimitations

DNA barcodes were assigned to BINs using the RESL algorithm implemented in BOLD. It has been shown that the RESL algorithm outperforms or is similar to other algorithms in terms of accuracy and computation time (Ratnasingham and Hebert 2013).

Species delimitation for SNP data was done with BFD* implemented in SNAPP, BEAST 2.6 (Leaché et al. 2014). The 2,434 unlinked SNPs from the *de novo* assembly of 25 specimens (10% missing data) were used as input (Table 1). The *n* = 84 datasets could not be used for the analysis due to computational resources limitations. Altogether 6 models were run and Bayes Factor (BF) was calculated for each to compare with the current taxonomy, model 1 (Table 3). The parameters for each model run included: mutation rate *U* = mutation rate *V* = 1, initial coalescent rate 10, *α*, *β*, and *λ*—default, burn-in 10%, 48 steps with 400k MCMC chain length each. The convergence of all runs was inspected with program Tracer 1.7.1 (Rambaut et al. 2018).

## Four-taxon D-statistics

Patterson's D-statistics is a test for archaic admixture (Green et al. 2010, Durand et al. 2011). In essence, it computes the excess of alleles in 1 species (usually denoted as P3) in relation to 2 sister species (P1 and P2) and the fourth outgroup species is a reference (O). Dsuite (Malinsky et al. 2021) was used to compute D-statistics for De_novo_84 and Reference_84 datasets with guide species trees and default parameters. The guide species trees were computed with SVDquartets (Chifman and Kubatko 2014) for both datasets using unlinked SNPs matrices (6,599 and 7,694 SNPs for De_novo_84 and Reference_84, respectively) and 100,000 quartets were estimated for each matrix. Next, the *f*-branch test (Malinsky et al. 2018) was performed for each dataset to assess the branches involved in hybridization using the Fbranch function implemented in the Dsuite program.

<table><caption>**Table 1.** Assemblies summary. Reference and *de novo* refer to the type of assembly. Header, in assembly: *n*—number of samples used; *m*—minimum taxon coverage parameter; *loci*—number of loci; SNPs—total number of SNPs; uSNPs—number of unlinked SNPs; Missing SNPs—percentage of missing data in SNPs matrix; missing sequence—total percentage of missing data over all loci</caption><thead><tr><th>Assembly</th><th>n</th><th>m</th><th>Loci</th><th>SNPs</th><th>uSNPs</th><th>Missing SNPs, %</th><th>Missing sequence, %</th></tr></thead><tbody><tr><td>Reference_84</td><td>84</td><td>43</td><td>7,784</td><td>109,247</td><td>7,694</td><td>35.87</td><td>41.34</td></tr><tr><td>De_novo_91</td><td>91</td><td>60</td><td>2,540</td><td>43,873</td><td>2,538</td><td>26.28</td><td>25.74</td></tr><tr><td>De_novo_84</td><td>84</td><td>43</td><td>6,608</td><td>113,479</td><td>6,599</td><td>33.38</td><td>32.83</td></tr><tr><td>De_novo_21</td><td>25</td><td>21</td><td>2,454</td><td>34,289</td><td>2,434</td><td>10.18</td><td>9.40</td></tr></tbody></table>

<table><caption>**Table 2.** The description of the datasets used in the study and their use in the downstream analyses. *P. hyp.*—*P. hyperborea*; *P. salt.*—*P. saltuaria*; *P. oreo.*—*P. oreophila*; *P. pal.*—*P. palustris*. The numbers indicate the number of specimens retained in the dataset; Under Analysis header there is a list of analyses applied to the corresponding dataset. Detailed information about all the specimens is in Supplementary Table S4</caption><thead><tr><th>Assembly</th><th><i>P. hyp.</i>, Finland</th><th><i>P. hyp.</i>, Canada</th><th><i>P. hyp.</i>, Greenland</th><th><i>P. salt.</i></th><th><i>P. oreo.</i></th><th><i>P. pal.</i></th><th>Analysis</th></tr></thead><tbody><tr><td>Reference_84</td><td>24</td><td>0</td><td>16</td><td>17</td><td>21</td><td>6</td><td>ML trees, D-statistics, fastsimcoal2, PhyloNetworks</td></tr><tr><td>De_novo_91</td><td>24</td><td>5</td><td>16</td><td>17</td><td>21</td><td>8</td><td>ML trees and STRUCTURE</td></tr><tr><td>De_novo_84</td><td>24</td><td>0</td><td>16</td><td>17</td><td>21</td><td>6</td><td>ML trees, D-statistics, fastsimcoal2, PhyloNetworks</td></tr><tr><td>De_novo_21</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>ML trees, SNAPP</td></tr></tbody></table>
<table><caption>**Table 3.** Summary of species delimitation with BFD*, De_novo_21 dataset, 2,434 unlinked SNPs. N—model number; N_Spec—number of species in the model excluding outgroup; sum(ESS)—sum of effective sample size for marginal likelihoods estimates; MLE—marginal likelihood estimate; Rank—species delimitation model rank in relation to current taxonomy; BF—Bayes factor. The model with the highest rank is in bold. *P. palustris* is used as an outgroup in all models</caption><thead><tr><th>N</th><th>Model</th><th>N_Spec</th><th>Sum(ESS)</th><th>MLE</th><th>Rank</th><th>BF</th></tr></thead><tbody><tr><td>1</td><td>Current taxonomy, <i>P. hyperborea</i> from Greenland and Finland treated as 1 species.</td><td>3</td><td>31,156.33</td><td>-25,155.8</td><td>2</td><td>no</td></tr><tr><td>2</td><td>All populations (species) as separate species.</td><td>4</td><td>29,073.61</td><td>-24,302.7</td><td>1</td><td>1,706.2</td></tr><tr><td>3</td><td><i>P. hyperborea</i> and <i>P. saltuaria</i> as 1 species, <i>P. oreophila</i> is a separate species.</td><td>2</td><td>33,008.87</td><td>-26,279.5</td><td>5</td><td>-2,247.4</td></tr><tr><td>4</td><td>All representatives of saltuaria group treated as 1 species.</td><td>1</td><td>33,828.27</td><td>-27,723.5</td><td>6</td><td>-5,135.4</td></tr><tr><td>5</td><td>Greenland population as separate species, Finnish <i>P. hyperborea</i> combined with <i>P. saltuaria</i></td><td>3</td><td>31,284.13</td><td>-25,171.1</td><td>3</td><td>-30.6</td></tr><tr><td>6</td><td>Greenland population as separate species, all European lineages treated as 1 species.</td><td>2</td><td>33,269.97</td><td>-26,240.5</td><td>4</td><td>-2,169.4</td></tr></tbody></table>

## Admixture Estimates Through Demographic Inference

Demographic inference for scenarios with or without admixture was performed with fastsimcoal2 version 2.6 (Excoffier et al. 2013). The goal was first, to assess support for introgression events in comparison to no introgression scenario and second, to compare different introgression scenarios. Input files were created from the ipyrad output vcf files by converting them to a minor allele site frequency spectrum (MSFS) with easySFS program (https://github.com/isaacovercast/easySFS). Reference and *de novo* assemblies with 84 specimens and *m* = 43 were tested for each designed scenario.

Altogether 14 scenarios were tested. Scenarios followed 4 scheme types: (i) no admixture, 1 scenario; (ii) one admixture event after species/populations split, 5 scenarios; (iii) 2 admixture events after species/populations split, 5 scenarios; (iv) 1 admixture event that preceded species/populations splits, 3 scenarios (fastsimcoal2 section, Supplementary Material). Each scenario was run 100 times with 1,000,000 simulations each. The only fixed parameter was divergence time between *P. palustris* and *P. saltuaria* group representatives (3,000,000 generations). The number of generations was approximated using 2 sources. First, by using the substitution rate for mtDNA in Lycosidae (Piacentini and Ramírez 2019) and applying it to computed *p*-distances (Supplementary Table S5) between *P. palustris* and ingroup species (approximately 1.5 m.y.a. which equals roughly to number of generations). In the same study, a time calibrated phylogenetic tree suggests divergence time between closely related *Pardosa* species groups to be 2–6 Myr. Thus, the value of 3 million generations, however arbitrary, was the estimate used for fastsimcoal2 runs. The rest of the parameters including effective population size ($N_e$), divergence (merge) time between species/populations ($t_m$), admixture time ($t_{adm}$), admixture percentage (adm) and mutation rate ($\mu$) were estimated by the program.

Maximum likelihoods for each scenario were compared using the Akaike Information Criterion (AIC) and the relative likelihood of the parameters was computed with exp((minimum AIC obser-ed— AIC for each model)/2), ΔAIC. Next, estimated parameters that had the highest maximum likelihood for each model were used as input for another run with fastsimcoal2 to obtain maximum likelihood estimates for each model (100 independent runs with 1,000,000 simulations each). Then, ΔAIC was computed to estimate the relative likelihood of each model.

## PhyloNetworks

To further explore possible admixture between the species, PhyloNetworks v.0.12.0 (Solís-Lemus et al. 2017) was run for the *De_novo* _84 and Reference_84 datasets with all recovered loci

(6,608 and 7,784 loci respectively). The PhyloNetworks estimates a possible hybrid origin of the edge taking into account phylogenetic relationships between edges. To prepare the input species tree, individual trees for each locus were computed with RAxML 8.2.12 (Stamatakis 2014) with HKY + Gamma model and 100 bootstrap (default parameters for the script; the script is available at https://github.com/crs14/InconsistencySpeciesTreeGeneFlow/blob/master/scripts/estGeneTrees/raxml.pl). Then the species tree was built with Astral 5.7.8 (Mirarab and Warnow 2015) with *P. hyperborea* from Greenland as a distinct species. Finally, networks were estimated by PhyloNetworks with 0–4 admixture events following the developers’ manual (https://crsl4.github.io/PhyloNetworks.jl/latest/). Each scenario was run 50 times. The best scenario was chosen based on the highest loglikelihood. Because the datasets contained missing data it was impossible to run bootstrap support for recovered hybrid edges.

## Nomenclature

This paper and the nomenclatural act(s) it contains have been registered in Zoobank (www.zoobank.org), the official register of the International Commission on Zoological Nomenclature. The LSID (Life Science Identifier) number of the publication is: urn:lsid:zoobank.org:pub:39DBAD94-6DC6-4F0E-A5CF-4B9B469A24A4.

## Results

### Assembly Statistics

Overall assembly statistics for the 4 datasets used in downstream analyses are in Table 1 and the number of specimens of each species included is in Table 2. The overall statistic for all assemblies are in Supplementary Table S3 and individual number of loci for each specimen and the percentage of missing loci are in Supplementary Table S4. Due to the low number of reads, 29 samples were completely excluded from ddRADseq matrices. They included practically all *P. hyperborea* specimens from the mainland Canada (Supplementary Tables S1 and S4). Thus, datasets with *n* = 84 were used for the majority of downstream analyses. Only the *De_novo* _91 dataset used for STRUCTURE and ML tree inference included 5 Canadian samples of *P. hyperborea* with lowest fraction of loci of 2.9% out of total loci count. As expected, *m* value in ipyrad pipeline was the main reason for the decrease in the number of recovered loci (Supplementary Table S3).

### ML Trees

ML trees were inferred for all 4 assemblies. All trees show a similar pattern where each ingroup species and populations of *P. hyperborea* from Greenland (and Canada for the 91 samples datasets) formed
a monophyletic group (Fig. 2b–d). The only discrepancy in comparison to all other trees was observed for the Reference_84 dataset. There, *P. saltuaria* and *P. oreophila* were clustered together while in all the rest of the trees *P. oreophila* was a sister lineage to both *P. hyperborea* populations and *P. saltuaria* (Fig. 2d). Except for the Reference_84 dataset tree bootstrap supports for each large clade were 100.

The ML tree for the COI dataset showed different results with 3 major lineages: the first was mixed consisting of *P. saltuaria*, *P. oreophila* from the Alps and *P. hyperborea* from Finland and Norway, while *P. oreophila* from the Pyrenees and *P. hyperborea* from North America were 2 other monophyletic clades (Fig. 2a). Specimens of *P. oreophila* from the Pyrenees were a sister lineage to the other 2.

## Geographic Structure of Genetic Variation

The COI haplotype network showed similar patterns as the ML tree for this locus. Namely, Alpine, central, and northern European lineages showed close affinity, while the Pyrenees formed a separate lineage and North American barcodes were completely independent from the European ones (Fig. 3).

STRUCTURE analysis based on SNP data unequivocally suggested a 5 cluster partition. The *K* = 5 received the highest value of Δ*K* with a logarithm probability of the data much lower for partitions with fewer clusters and it steadily diminished for 6 and 7

clusters (Supplementary Fig. S1). Admixture was observed between western and central Canadian and Finnish *P. hyperborea* samples, while samples from Greenland and northern Canada did not show any signs of shared ancestry with European lineages (Fig. 4).

## Species Delimitation

All analyzed DNA barcodes were assigned to 5 BINs. There were 2 BINs reported for specimens identified as *P. hyperborea* in North America (BOLD:AAA6658, common and BOLD:ACE5253, rare) (Blagoev et al. 2016). The common type was found throughout the country and the rare type occurred in the center and west. Greenland samples belonged to the common BIN. The split between North American haplotypes was also apparent on the haplotype network with a minimum number of mutations between species from different BINs equal to 11 (Fig. 3). In Europe, BINs followed previously reported patterns (Muster and Berendonk 2006, Muster et al. 2009) where Scandinavian, Alps and Tatra mountain species—*P. oreophila*, *P. saltuaria*, *P. hyperborea*—belonged to BOLD:ACK1850, while *P. oreophila* from the Pyrenees formed its own BIN (BOLD:AED0758).

The BFD* analysis supported a model with a 5 species partition (Model 2, Table 3) where Finnish *P. hyperborea*, Greenland *P. hyperborea*, *P. saltuaria*, and French *P. oreophila* were assessed as separate species. The BF of 1,706.2 was overwhelmingly higher for the maximum split model than for any other model merging the species. It should be noted that 1–8 steps in each model run could

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Maximum likelihood trees: (a) COI, 426 specimens; (b) ddRADseq, 91 specimens, 43,873 SNPs, *de novo* assembly; (c) ddRADseq, 84 specimens, 113,479 SNPs, *de novo* assembly; (d) ddRADseq, 84 specimens, 109,247 SNPs. CA—Canada, FI—Finland, FR—France, GL—Greenland, NO—Norway, SL—Slovakia, US—United States of America, *n*—number of specimens, *m*—minimum taxon coverage in ipyrad assembly.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Haplotype network for 419 specimens of *P. saltuaria* species group based on mitochondrial COI data. Haplotypes are colored based on morphological identification. Each line is 1 mutation step.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Pie chart map of SNPs, showing genetic admixture of *Pardosa* species based on geographic location. STRUCTURE results of *K* = 5 based on 42,215 SNPs plotted in the sampling sites.</figcaption></figure>
not converge according to the individual Effective Sample Size (ESS). Nevertheless, sum(ESS) for marginal likelihood ranged from 29,073.61 to 33,828.26 suggesting that overall estimates were reliable. It was not feasible to run separately the un-converged steps as it would have taken months of heavy load on high performance computing cluster with a likely low contribution to the result. The species tree of the model with the highest BF can be found in Supplementary Fig. S2.

## D-statistics

The D-statistics for the De_novo_84 and the Reference_84 datasets suggested introgression between all European species but had different support values (D-statistics, Z-score, and P-value) for particular triplets (Table 4). In both datasets the excess of ABBA loci between *P. saltuaria* and *P. oreophila* had the highest support and similar statistics while introgression between *P. saltuaria* and Finnish *P. hyperborea* and *P. oreophila* and Finnish *P. oreophila* was significant only for the *de novo* dataset (Z-score > 3 and P-value << 0.01). For the reference dataset only introgression between *P. oreophila* and *P. saltuaria* was significant. The *f*<sub>4</sub>-statistics assessed admixture levels between 0.02 and 0.06 (2–6%). The *f*-branch analysis suggests introgression occurred between European lineages after the split between Finnish and Greenland *P. hyperborea* as no internal nodes showed any signs of admixture for both datasets (Fig. 5a and b).

## Demographic Scenarios

All scenarios were successfully simulated and ranked based on AIC except for scenario D1 for the *de novo* dataset (Table 5). Scenario D1, which assessed admixture between *P. saltuaria* and a common ancestor of both *P. hyperborea* populations, was constantly giving bad parameter estimates and lineages could not coalesce in designated order despite the high number of initiation attempts. Moreover, it showed low support for the reference dataset. Therefore, it was discarded from further discussion. The best estimated parameters for each scenario are presented in the fastsimcoal2 section of the Supplementary Material.

For the *de novo* dataset, the scenario describing 2 admixture events from *P. oreophila* to Finnish *P. hyperborea* and *P. saltuaria* (C3) was ranked the highest (Table 5). The probability of the second-best scenario (1 admixture event between *P. oreophila* and *P. saltuaria*, B5) was 7 orders lower in magnitude. The best estimated parameters from the C3 scenario suggested recent admixture in comparison to species merges and a lower level of admixture to *P. hyperborea* than to *P. saltuaria* (10 and 1% correspondingly, Table 5 and Fig. 5c).

For the reference dataset, 3 scenarios scored almost equally: B1, C2, and C4 with C2 having the lowest AIC (Fig. 5d and fastsimcoal2 section, Supplementary Material). All of them suggested admixture

between *P. saltuaria* and Finnish *P. hyperborea* and, in addition, between *P. saltuaria* and *P. oreophila* from the Pyrenees (C2) and between *P. saltuaria* and *P. hyperborea* from Greenland (C4) (Table 5). However, the order of coalescent events was violated as the merge between *P. saltuaria* and *P. oreophila* was estimated to happen earlier than expected admixture events which contradicts the phylogenetic order of events (Fig. 2). This happened possibly due to high estimated admixture proportions (up to 24%) that forced lineages to coalesce faster as much of the variation was already shared. Overall, for both *de novo* and reference datasets, admixture within European lineages had the highest support.

## PhyloNetworks

The species trees computed in ASTRAL supported delimitation inferred by BFD* with a local posterior probability of 1.0 for all branches (tree scores for the De_novo_84 and Reference_84 datasets are 0.6 and 0.53 respectively, Supplementary Figs. S3 and S4). The PhyloNetworks estimate of the highest number of admixture events was 1 for both datasets (*de novo* and reference -loglikelihood 0.013 and 0.004 respectively, Supplementary Figs. S5 and S6) suggesting a hybrid origin of *P. saltuaria* in both cases. In De_novo_84, *P. saltuaria* was a 50% sister to *P. oreophila* and a 50% sister to Finnish *P. hyperborea* (Fig. 5e) and in Reference_84 corresponding values were 64.6 and 35.4% (Fig. 5f).

## Discussion

The data analyzed allows us to answer the main questions posed. Firstly, the admixture previously observed among European species of the *saltuaria* group is confirmed by 3 analyses (D-statistics, fastsimcoal2, and PhyloNetworks) but is likely explained by historical rather than current gene flow. Secondly, both ddRADseq and DNA barcodes inform species delimitation but possible introgression of mtDNA can cause complications in the delineation of recently diverged species with complex demographic histories. Thirdly, there is enough genomic evidence to confirm the species status of North American and Greenland populations of currently recognized *P. hyperborea*. Below we discuss these conclusions in light of relevant literature and our findings.

## Taxonomic Status of North American *P. hyperborea*

We demonstrate that North American and Greenland populations of *P. hyperborea* are genetically distinct from European populations to the extent suggesting that they should be regarded as different species. For specimens from Canada, the name *P. luteola* Emerton, 1894 is available (Emerton 1894). The syntype is stored in the Museum of Comparative Zoology and was collected from Lake Louise, Alberta, Canada. Both North American BINs were present

<table><caption>**Table 4.** D-statistics summary for SNP datasets for *n* = 84, *m* = 43 *de novo* (DEN) and reference (REF) datasets. D—D-statistics for the triplet; Z—Z-score for the triplet; P-value—P-value for the triplet; f4-ratio—f4-ratio for the triplet BBAA, ABBA, BABA—number of loci supporting each of the triplet pattern</caption><thead><tr><th rowspan="2">P1</th><th rowspan="2">P2</th><th rowspan="2">P3</th><th colspan="2">D</th><th colspan="2">Z</th><th colspan="2">P</th><th colspan="2">f4</th><th colspan="2">BBAA</th><th colspan="2">ABBA</th><th colspan="2">BABA</th></tr><tr><th>DEN</th><th>REF</th><th>DEN</th><th>REF</th><th>DEN</th><th>REF</th><th>DEN</th><th>REF</th><th>DEN</th><th>REF</th><th>DEN</th><th>REF</th><th>DEN</th><th>REF</th></tr></thead><tbody><tr><td>HypG</td><td>HypF</td><td>Oreo</td><td>0.10</td><td>0.04</td><td>3.96</td><td>1.84</td><td>3.80E-05</td><td>0.033</td><td>0.04</td><td>0.02</td><td>1,315.17</td><td>1,182.98</td><td>658.00</td><td>708.87</td><td>543.15</td><td>657.70</td></tr><tr><td>HypG</td><td>HypF</td><td>Salt</td><td>0.09</td><td>0.05</td><td>3.61</td><td>2.37</td><td>1.53E-04</td><td>0.009</td><td>0.04</td><td>0.03</td><td>1,044.03</td><td>1,015.03</td><td>718.66</td><td>753.99</td><td>602.83</td><td>685.75</td></tr><tr><td>HypF</td><td>Salt</td><td>Oreo</td><td>0.04</td><td>0.07</td><td>1.61</td><td>3.54</td><td>0.054</td><td>2.03E-04</td><td>0.02</td><td>0.04</td><td>1,002.47</td><td>934.09</td><td>730.75</td><td>831.16</td><td>674.02</td><td>718.31</td></tr><tr><td>HypG</td><td>Salt</td><td>Oreo</td><td>0.11</td><td>0.11</td><td>4.16</td><td>4.71</td><td>1.59E-05</td><td>1.23E-06</td><td>0.06</td><td>0.06</td><td>940.31</td><td>856.69</td><td>785.20</td><td>831.13</td><td>624.06</td><td>672.23</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Introgression in *P. saltuaria* species group, *n* = 84, *m* = 43 datasets. *Phyper.FIN*—Finnish *P. hyperborea*, *Phyper.GR*—Greenland *P. hyperborea*; (a, c, e) *de novo* assembly; (b, d, f) reference assembly. (a, b) Fbranch D-statistics summary: dotted lines indicate internal node branches, *fb*—f-branch estimates of admixture; (c, d) fastsimcoal2 scenarios with the highest support: Ne effective population size of corresponding species/population (see fastsimcoal2 section in *Supplementary Material* for details about codes in brackets); *t*<sub>m</sub>—time of population/species merge (split) in generations; *t*<sub>adm</sub>—time of admixture in generations; adm—admixture fraction; (e, f) PhyloNetworks: colored lines and numbers indicate direction and proportion of admixture.</figcaption></figure>
<table><caption>**Table 5.** Fastsomcoal2 model simulation. HypF—*P. hyperborea* from Finland; HypG—*P. hyperborea* from Greenland; S—*P. saltuaria*; O—*P. oreophila*; CA() common ancestor of the pair in brackets; adm—admixture; Likelihood—best maximum estimated likelihood over 100 simulations; AIC—Akaike criterion; Exp(ΔAIC)—relative likelihood of the scenario in comparison with the lowest AIC scenario</caption><thead><tr><th rowspan="2">Scenario</th><th rowspan="2">Description</th><th colspan="3">De novo</th><th colspan="3">Reference</th></tr><tr><th>Likelihood</th><th>AIC</th><th>Exp(ΔAIC)</th><th>Likelihood</th><th>AIC</th><th>Exp(ΔAIC)</th></tr></thead><tbody><tr><td>A</td><td>No admixture</td><td>-2,145.4</td><td>9899.9</td><td>5E-19</td><td>-1,350.3</td><td>6,238.6</td><td>2.29E-21</td></tr><tr><td>B1</td><td>1 adm, HypF-S</td><td>-2,134.2</td><td>9852.6</td><td>9.67E-09</td><td><strong>-1,328.9</strong></td><td><strong>6,143.8</strong></td><td><strong>0.883</strong></td></tr><tr><td>B2</td><td>1 adm, HypF-O</td><td>-2,134.1</td><td>9851.8</td><td>1.43E-08</td><td>-1,334.8</td><td>6,170.9</td><td>1.13E-06</td></tr><tr><td>B3</td><td>1 adm, HypG-S</td><td>-2,155.0</td><td>9948.2</td><td>1.7E-29</td><td>-1,345.1</td><td>6,218.5</td><td>5.29E-17</td></tr><tr><td>B4</td><td>1 adm, HypG-O</td><td>-2,135.6</td><td>9858.8</td><td>4.29E-10</td><td>-1,331.9</td><td>6,157.5</td><td>0.0009</td></tr><tr><td>B5</td><td>1 adm, O-S</td><td>-2,132.5</td><td>9844.4</td><td>5.73E-07</td><td>-1,330.9</td><td>6,152.8</td><td>0.0096</td></tr><tr><td>C1</td><td>2 adm, HypF-S &amp; O</td><td>-2,150.6</td><td>9931.8</td><td>6.05E-26</td><td>-1,339.8</td><td>6,198.2</td><td>1.36E-12</td></tr><tr><td>C2</td><td>2 adm, S-HypF &amp; O</td><td>-2,134.0</td><td>9855.4</td><td>2.35E-09</td><td><strong>-1,328.0</strong></td><td><strong>6,143.5</strong></td><td><strong>1</strong></td></tr><tr><td>C3</td><td>2 adm, O-HypF &amp; S</td><td><strong>-2,125.4</strong></td><td><strong>9815.7</strong></td><td><strong>1</strong></td><td>-1,330.9</td><td>6,157.2</td><td>0.0011</td></tr><tr><td>C4</td><td>2 adm, S-HypF &amp; HypG</td><td>-2,133.2</td><td>9851.8</td><td>1.42E-08</td><td><strong>-1,328.9</strong></td><td><strong>6,148.0</strong></td><td><strong>0.1068</strong></td></tr><tr><td>C5</td><td>2 adm, O-HypF &amp; HypG</td><td>-2,130.9</td><td>9841.0</td><td>3.10E-06</td><td>-1,332.2</td><td>6,163.1</td><td>5.66E-05</td></tr><tr><td>D1</td><td>1 adm, S-CA(HypF &amp; HypG)</td><td>NA</td><td>NA</td><td>NA</td><td>-1,349.6</td><td>6,239.3</td><td>1.55E-21</td></tr><tr><td>D2</td><td>1 adm, O-CA(HypF &amp; HypG)</td><td>-2,139.0</td><td>9874.6</td><td>1.57E-13</td><td>-1,338.0</td><td>6,185.6</td><td>7.28E-10</td></tr><tr><td>D3</td><td>1 adm, CA(O &amp; S)-CA(HypF &amp; HypG)</td><td>-2,143.2</td><td>9893.7</td><td>1.15E-17</td><td>-1,342.2</td><td>6,205.2</td><td>3.96E-14</td></tr></tbody></table>

in the vicinity of the syntype location. There is also a species description by Strand under the name *Lycosa albimontis* Strand, 1909 which was later synonymized with *P. hyperborea*. The spiders were collected in the White Mountains, New Hampshire, USA (Strand 1909). It is likely that specimens from this area belong to the more widespread BIN found in eastern and central Canada and Greenland (BOLD:AAA6658) because the rarer BIN (BOLD:ACE5253) has never been detected east of Alberta.

Taking into account the geographical information, DNA barcoding (BINs) and ddRADseq evidence (ML trees and SNAPP) it could be possible to resurrect both synonyms. A simple solution would be to assign *P. luteola* to the western BIN as it is geographically closer to the syntype location and the name *L. albimontis* to the second BIN. The taxonomic proposal presented above has its complications. First, there could be more species within North American *P. hyperborea* populations that can be detected only through the genetic approach due to the mitonuclear discordance known for the genus (Ivanov et al. 2018). Second, specimens from Canada were intermixed with each other and Greenland specimens in the ML trees calculated for the ddRADseq datasets. We suggest that with the current state of knowledge a more conservative decision is appropriate. Specifically, we propose to assign the name *P. luteola* Emerton, 1894 stat. resurr. to both mtDNA lineages in North America.

Reconstruction of historical migration routes was not the goal of the study. However, several observations are worth mentioning. First, *P. californica*—the second species of the group in North America—is found only in the west (Dondale and Redner 1987, Vogel 2004) as well as the rarer lineage of *P. hyperborea* (BIN BOLD:ACE5253). Second, the STRUCTURE results suggest shared ancestry between the Finnish *P. hyperborea* population and several *P. hyperborea* specimens from mainland Canada in the center and west of the country but not with Greenland specimens (Fig. 4). Third, *p*-distances in COI between European *P. hyperborea* (and thus the whole mixed cluster from the Alps to the European Arctic) and the western BIN in North America are shorter than with the BIN found in both Canada and Greenland (Supplementary Table S5). These observations might

suggest possible historical migratory routes of the *saltuaria* group from Eurasia to North America. Without additional research and especially broader sampling (including Iceland and eastern Greenland) these questions are impossible to answer. Sequencing *P. oksalai* from Russia as well as resequencing of our samples for genetic markers with a suitable mutation rate would help to shed light on this problem.

## Taxonomy of European *Saltuaria* Group Species

The main taxonomic conclusion for European *P. hyperborea* and *P. saltuaria* is that the present taxonomy remains valid despite shared COI. However, there is an outpost population of *P. hyperborea* in the Bohemian Forest in the Czech Republic that was not included in the study and might have an impact on the taxonomic decision. Previously it was shown that specimens from this area were genetically mixed with Scandinavian lineages which in their turn, are mixed with *P. saltuaria* and *P. oreophila* from the Alps (Muster and Berendonk 2006).

The presence of clear differences in DNA barcodes and ddRADseq evidence in *P. oreophila* from the Pyrenees—ML trees, BFD*, STRUCTURE—allows the elevation of this population to species status. Our results agree with previous research done in the group (Muster and Berendonk 2006). However, for formal species descriptions, we would require male representatives to be sequenced at least for COI as well as to be examined morphologically. In addition, it is crucial to test the species boundaries for the Alpine *P. oreophila* against the Pyrenean population with ddRAD or a similar method to clarify the status of both populations.

The taxonomy of European species remains partly unresolved due to the lack of samples from some important areas. The main problem of the current study for European populations is the absence of *P. oreophila* samples from the Alps in the ddRADseq sampling. In addition, other species of the group (*P. evelinae*, *P. drenskii*, *P. oksalai*) could have contributed greatly to the overall assessment of taxonomy. Bearing in mind the aforementioned limitations, we nevertheless suggest that the results of the current species delimitation
analysis will not change irrespective of species or populations. The reason lies in the allopatric distribution that forces researchers to make largely subjective taxonomic decisions as gene flow between populations is interrupted and testing mating compatibility in the laboratory does not reflect the natural situation.

## Historical Admixture

Previous studies of species of the saltuaria group have suggested current migration or a single migration burst that happened ca 20k yr ago between Scandinavia and central Europe but not between the Pyrenees and the rest of the continent (Muster et al. 2009). Our analyses show the presence of admixed nDNA in all tested European species but not in specimens from Canada. The fastsimcoal2 estimates of introgression time between P. hyperborea and P. saltuaria are comparable to previous studies. For the De_novo_84 datasets it was 14–21k generations ago and 24–267k generations ago for the Reference_84 dataset (scenarios B1, C1, C2, and C4, fastsimcoal2 section in Supplementary Material). However, the C3 scenario that has higher support for the de novo dataset does not include such admixture events at all favoring introgression from P. oreophila to P. hyperborea and P. saltuaria (Fig. 5c). On the contrary, scenarios with the highest support for reference dataset always included gene exchange between P. saltuaria and European P. hyperborea (B1, C2, and C4, fastsimcoal2 section, Supplementary Material). Simultaneously, the D-statistics and PhyloNetworks analysis results show a similar picture of gene flow between P. hyperborea and P. saltuaria for both datasets. The estimates of admixed nDNA differ between the methods, ranging from 3% to 50% (Fig. 5). These programs suggest that gene exchange occurred between the 2 species sharing COI after the split between the Finnish and North American P. hyperborea. The f-branch analysis did not detect introgression on internal node branches favoring the same idea. The PhyloNetworks results suggest that P. saltuaria received a third or a half of genomic DNA from the sister lineage of the current P. hyperborea, also after the Europe-North America split. Thus, the scenario that is suggested by our data is that North American species split from European ca 600k generations ago and European P. hyperborea and P. saltuaria had at least 1 secondary contact ca 23k generations ago. The timing is based on the median of all scenarios estimates with fastsimcoal2. The shared mtDNA could have become introgressed during gene exchange events between these 2 species. The latter contradicts a previous suggestion that shared haplotypes between P. saltuaria and P. hyperborea were the result of the retention of ancestral genes in large population (Muster and Berendonk 2006).

Introgression of Scandinavian and Carpathian species with the Pyrenean *P. oreophila* was supported by D-statistics and demographic simulations but not by phylogenetic networks. The amount of shared genetic material between the species ranged from 1 to 10%. The admixture event(s) was suggested to happen after *P. oreophila* split from the common ancestor of *P. hyperborea* and *P. saltuaria*. In this case, however, mtDNA did not become introgressed together with nDNA.

The previous research proposed that the Pyrenean species was trapped at high elevations after several glaciation events that happened during the early Pleistocene (Muster and Berendonk 2006). The fastsimcoal2 partially confirms this hypothesis, estimating the split of *P. oreophila* from the common ancestor of *P. hyperborea* and *P. saltuaria* to have happened over 2 million and ca 1 million generations ago for the *de novo* and reference datasets respectively. The scenarios with the highest support—where introgression with *P. oreophila* is present—suggest relatively recent introgression for the *de novo* dataset (33k and 90k generations ago) while for the

reference dataset the estimate is over 1 million generations ago which violates the order of coalescent events. Therefore, it is difficult to assess the actual timing of the secondary contact. If it was a single event, then gene exchange is likely to have occurred during colder periods in the late Pleistocene (Last Glacial Maximum) or during the Riß glaciation when the cold adapted saltuaria species group expanded its range.

Historical connections between the Pyrenees and the Carpathian are known for plants and animals ranging from virtually no genetic differentiation between regions to various levels of admixture within species (Schmitt et al. 2010). However, to our best knowledge, cases of introgression between deeply diverged species from the Pyrenees and central and northern Europe have never been reported before. Thus, isolated southern species and populations of cold adapted organisms can bear signatures of earlier admixture events and genome-level sequencing can help to reconstruct the demographic dynamics of focal species at deeper time scales. In the case of the saltuaria species group, additional sampling from the Alps, eastern Europe, and Asia is crucial to have a better understanding of the historical dynamics of species and their interactions.

## Conclusions

Our results strongly suggest that *P. hyperborea* populations on both sides of the Atlantic Ocean do not belong to the same species. Therefore, we reinstate *Pardosa luteola* as a valid species name for both mitochondrial lineages currently assigned to *Pardosa hyperborea* that occur in North America. Furthermore, it is likely that the observed split within North America is recurrent in many other species given that similar interglacial dynamics shaped modern boreal fauna in the Nearctic region. We predict that genomics studies will increasingly reveal many presumably Holarctic organisms to be in fact complexes of genetically distinct sister species, with ranges of the species largely following a Palearctic-Nearctic division, although with exceptions, particularly on both sides of the Beringian Strait.

The ddRADseq proved to be a valuable tool for species delimitation in the complex of European saltuaria species. Our study further highlights that more work should be done to disentangle the historical relationships of Arctic-Alpine species of *Pardosa*. Secondary contacts and introgression are likely to be common among allopatric cold adapted species but to fully uncover their history broader sampling and genome sequencing are required. Nevertheless, ddRADseq data analyses confirm previously suggested admixture events and demonstrated unexpected patterns of introgression between highly differentiated species.

## Supplementary Material

Supplementary material is available at *Insect Systematics and Diversity* online.

## Specimen Collection Statement

No specimens were included in this project that required permits of any kind.

## Acknowledgments

We thank Laura Törmälä for invaluable support in the laboratory, 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 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 the Academy of Finland (grants #277984 and #
314702 accorded to MM), the Kvantum Institute of the University of Oulu (a research grant accorded to MM) and the Helsinki Entomological Society (personal grant to VI).

## Author Contributions

Vladislav Ivanov (Conceptualization [Equal], Data curation [Equal], Formal analysis [Lead], Funding acquisition [Supporting], Investigation [Lead], Validation [Equal], Visualization [Equal], Writing—original draft [Lead]), Gergin Blagoev (Data curation [Equal], Investigation [Equal], Writing—original draft [Equal]), Samuel Danflous (Data curation [Equal], Investigation [Equal], Writing—original draft [Equal]), Peter Gajdoš (Data curation [Equal], Investigation [Equal], Writing—original draft [Equal]), Toke Høye (Data curation [Equal], Investigation [Equal], Writing—original draft [Equal]), Kyung Min Lee (Formal analysis [Equal], Visualization [Equal], Writing—original draft [Equal]), Yuri Marusik (Data curation [Equal], Investigation [Equal], Visualization [Supporting], Writing—original draft [Equal]), Cecilie Mielec (Data curation [Equal], Investigation [Equal], Writing—original draft [Equal]), Christoph Muster (Data curation [Equal], Investigation [Equal], Writing—original draft [Equal]), Julien Pétillon (Data curation [Equal], Investigation [Equal], Writing—original draft [Equal]), Jörg Spelda (Data curation [Equal], Investigation [Equal], Writing—original draft [Equal]), and Marko Mutanen (Conceptualization [Equal], Data curation [Equal], Funding acquisition [Lead], Project administration [Equal], Supervision [Equal], Writing—original draft [Equal])

## Data Availability

Raw ddRADseq reads: GenBank SRA, BioProject IDs PRJNA345307 and PRJNA595572. DNA barcodes are publicly available in BOLD, DS-HYPERB dataset.

## References

* Abel C, Schneider JM, Kuntner M, Harms D. Phylogeography of the "cosmopolitan" orb-weaver *Argiope trifasciata* (Araneae: Araneidae). Biol J Linn Soc. 2020;131(1):61–75. https://doi.org/10.1093/biolinnean/blaa078

* Andrews S. FastQC: a quality control tool for high throughput sequence data; 2010.

* Astrin JJ, Höfer H, Spelda J, Holstein J, Bayer S, Hendrich L, Huber BA, Kielhorn KH, Krammer HJ, Lemke M, et al. Towards a DNA barcode reference database for spiders and harvestmen of Germany. PLoS One. 2016;11(9):e0162624–e0162624. https://doi.org/10.1371/journal.pone.0162624

* Bidegaray-Batista L, Ferrández M, Arnedo MA. Winter is coming: miocene and quaternary climatic shifts shaped the diversification of Western-Mediterranean *Harpactocrates* (Araneae, Dysderidae) spiders. Cladistics. 2014;30(4):428–446. https://doi.org/10.1111/cla.12054

* Blagoev GA, DeWaard JR, Ratnasingham S, DeWaard SL, Lu L, Robertson J, Telfer AC, Hebert PDN. Untangling taxonomy: a DNA barcode reference library for Canadian spiders. Mol Ecol Resour. 2016;16:325–341.

* Bonnet T, Leblois R, Rousset F, Crochet P-A. A reassessment of explanations for discordant introgressions of mitochondrial and nuclear genomes. Evolution (NY). 2017;71(9):2140–2158. https://doi.org/10.1111/evo.13296

* Brunke AJ, Salnitska M, Hansen AK, Zmudzinska A, Smetana A, Buffam J, Solodovnikov A. Are subcortical rove beetles truly Holarctic? An integrative taxonomic revision of north temperate *Quedionuchus* (Coleoptera: Staphylinidae: Staphylininae). Org Divers Evol. 2020;20(1):77–116. https://doi.org/10.1007/s13127-019-00422-2

* Chhatre VE, Emerson KJ. StrAuto: automation and parallelization of STRUCTURE analysis. BMC Bioinf. 2017;18(1): 192. https://doi.org/10.1186/s12859-017-1593-0

* Chifman J, Kubatko L. Quartet inference from SNP data under the coalescent model. Bioinformatics. 2014;30(23):3317–3324. https://doi.org/10.1093/bioinformatics/btu530

* Clement M, Snell Q, Walke P, Posada D, Crandall K. TCS: estimating gene genealogies. In: Proceedings of the 16th International Parallel and Distributed Processing Symposium. IEEE; 2002. p. 7.

* DaCosta JM, Sorenson MD. Amplification biases and consistent recovery of loci in a double-digest RAD-seq protocol. PLoS One. 2014;9(9):e106713. https://doi.org/10.1371/journal.pone.0106713

* Dincă V, Lee KM, Vila R, Mutanen M. The conundrum of species delimitation: a genomic perspective on a mitogenetically super-variable butterfly. Proc R Soc B Biol Sci. 2019;286(1911):20191311. https://doi.org/10.1098/rspb.2019.1311

* Divakar PK, Wei X, McCune B, Cubas P, Boluda CG, Leavitt SD, Crespo A, Tchabanenko S, Lumbsch HT. Parallel Miocene dispersal events explain the cosmopolitan distribution of the Hypogymnioid lichens. J Biogeogr. 2019;46(5):945–955. https://doi.org/10.1111/jbi.13554

* Dondale CD, Redner JH. The *atrata*, *cubana*, *ferruginea*, *moesta*, *monticola*, *saltuaria*, and *solituda* groups of the spider genus *Pardosa* in North America (Araneae: Lycosidae). Can Entomol. 1987;119(1):1–19. https://doi.org/10.4039/ent1191-1

* Dufresnes C, Nicieza AG, Litvinchuk SN, Rodrigues N, Jeffries DL, Vences M, Perrin N, Martínez-Solano I. Are glacial refugia hotspots of speciation and cytonuclear discordances? Answers from the genomic phylogeography of Spanish common frogs. Mol Ecol. 2020;29(5): 986–1000. https://doi.org/10.1111/mec.15368

* Dufresnes C, Rodrigues N, Savary R. Slow and steady wins the race: contrasted phylogeographic signatures in two Alpine amphibians. Integr. Zool. 2021;17(1): 181–190. https://doi.org/10.1111/1749-4877.12518

* Durand EY, Patterson N, Reich D, Slatkin M. Testing for ancient admixture between closely related populations. Mol Biol Evol. 2011;28(8):2239–2252. https://doi.org/10.1093/molbev/msr048

* Earl DA, vonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4(2):359–361. https://doi.org/10.1007/s12686-011-9548-7

* Eaton DAR, Overcast I. ipyrad: Interactive assembly and analysis of RADseq datasets. Bioinformatics. 2020;36(8):2592–2594. https://doi.org/10.1093/bioinformatics/btz966

* Elias-Costa AJ, Confalonieri VA, Lanteri AA, Rodriguero MS. Game of clones: Is *Wolbachia* inducing speciation in a weevil with a mixed reproductive mode?. Mol Phylogenet Evol. 2019;133:42–53. https://doi.org/10.1016/j.ympev.2018.12.027

* Emerton JH. Canadian spiders. Trans Connect Acad Arts Sci. 1894;9:400–429.

* Enghoff H. Historical biogeography of the Holarctic: area relationships, ancestral areas, and dispersal of non-marine animals. Cladistics. 1995;11(3):223–263. https://doi.org/10.1111/j.1096-0031.1995.tb00088.x

* Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software structure: a simulation study. Mol Ecol. 2005;14(8):2611–2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x

* Excoffier L, Dupanloup I, Huerta-Sánchez E, Sousa VC, Foll M. Robust demographic inference from genomic and SNP data. PLoS Genet. 2013;9(10):e1003905. https://doi.org/10.1371/journal.pgen.1003905

* Fiorentino V, Manganelli G, Giusti F, Ketmaier V. Recent expansion and relic survival: phylogeography of the land snail genus *Helix* (Mollusca, Gastropoda) from south to north Europe. Mol Phylogenet Evol. 2016;98:358–372. https://doi.org/10.1016/j.ympev.2016.02.017

* Funk DJ, Omland KE. Species-level paraphyly and polyphyly: frequency, causes, and consequences, with Insights from animal mitochondrial DNA. Annu Rev Ecol Syst. 2003;34(1):397–423. https://doi.org/10.1146/annurev.ecolsys.34.011102.132421

* Gilligan TM, Brown JW, Baixeras J. Immigrant Tortricidae: Holarctic versus introduced species in North America. Insects. 2020;11(9):594. https://doi.org/10.3390/insects11090594

* Green RE, Krause J, Briggs AW, Maricic T, Stenzel U, Kircher M, Patterson N, Li H, Zhai W, Fritz MH-Y, et al. A draft sequence of the Neandertal genome. Science. 2010;328:710–722.

* Guindon S, Dufayard J-F, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: Assessing the performance of PhyML 3.0. Syst Biol. 2010;59(3):307–321. https://doi.org/10.1093/sysbio/syq010

* György Z, Tóth EG, Incze N, Molnár B, Höhn M. Intercontinental migration pattern and genetic differentiation of arctic-alpine *Rhodiola rosea* L.: a chloroplast DNA survey. Ecol Evol. 2018;8(23):11508–11521. https://doi.org/10.1002/ece3.4589

* Hewitt G. The genetic legacy of the quaternary ice ages. Nature. 2000;405(6789):907–913. https://doi.org/10.1038/35016000
Hibbets EM, Schumacher KI, Scheppler HB, Boersma PD, Bouzat JL. Genetic evidence of hybridization between *Magellanicus* (*Sphensicus magellanicus*) and Humboldt (*Spheniscus humboldti*) penguins in the wild. *Genetica*. 2020;148(5-6):215-228. https://doi.org/10.1007/s10709-020-00106-2

Hoang DT, Chernomor O, von Haeseler A, Minh BQ, Vinh LS. UFBoot2: improving the ultrafast bootstrap approximation. Mol Biol Evol. 2018;35(2):518-522. https://doi.org/10.1093/molbev/msx281

Hou Y, Nowak MD, Mirré V, Bjorá CS, Brochmann C, Popp M. Thousands of RAD-seq loci fully resolve the phylogeny of the highly disjunct arctic-alpine genus *Diapensia* (Diapensiaceae). PLoS One. 2015;10(10):e0140175. https://doi.org/10.1371/journal.pone.0140175

Hundsdoerfer AK, Kitching IJ. Ancient incomplete lineage sorting of *Hyles* and *Rhodafra* (Lepidoptera: Sphingidae). Org Divers Evol. 2020;20(3):527-536. https://doi.org/10.1007/s13127-020-00445-0

Ivanov V, Lee KM, Mutanen M. Mitonuclear discordance in wolf spiders: genomic evidence for species integrity and introgression. Mol Ecol. 2018;27(7):1681-1695. https://doi.org/10.1111/mec.14564

Ivanov V, Marusik Y, Pétillon J, Mutanen M. Relevance of ddRADseq method for species and population delimitation of closely related and widely distributed wolf spiders (Araneae, Lycosidae). Sci Rep. 2021;11(1):1-14.

Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23(14):1801-1806. https://doi.org/10.1093/bioinformatics/btm233

Johnsen A, Rindal E, Ericson PGP, Zuccon D, Kerr KCR, Stoeckle MY, Lifjeld JT. DNA barcoding of Scandinavian birds reveals divergent lineages in trans-Atlantic species. J Ornithol. 2010;151(3):565-578. https://doi.org/10.1007/s10336-009-0490-3

Khadzhieva UA, Dragan SV, Pospelov IG. Materials to the fauna of spiders (Arachnida: Araneae) in the Ergaki Natural Park. Proc Mordovia State Nat Reserv. 2019;22:186-192.

Kraus O. Why no subspecies in spiders? In: Toft S, Scharff N, editors. Proc. 19th Eur. Colloq. Århus (Denmark): Arachnol. Aarhus University Press; 2002. p. 303-314.

Landry J-F, Nazari V, de Waard JR, Mutanen M, Lopez-Vaamonde C, Huemer P, Hebert PDN. Shared but overlooked: 30 species of Holarctic Microlepidoptera revealed by DNA barcodes and morphology. Zootaxa. 2013;3749-1.

Leaché AD, Fujita MK, Minin VN, Bouckaert RR. Species delimitation using genome-wide SNP Data. Syst Biol. 2014;63(4):534-542. https://doi.org/10.1093/sysbio/syu018

Leduc-Robert G, Maddison WP. Phylogeny with introgression in *Habronattus* jumping spiders (Araneae: Salticidae). BMC Evol Biol. 2018;18(1):24. https://doi.org/10.1186/s12862-018-1137-x

Leigh JW, Bryant D. Popart: full-feature software for haplotype network construction. Methods Ecol Evol. 2015;6(9):1110-1116. https://doi.org/10.1111/2041-210x.12410

Linck E, Epperly K, Van Els P, Spellman GM, Bryson RW, McCormack JE, Canales-Del-Castillo R, Klicka J. Dense geographic and genomic sampling reveals paraphyly and a cryptic lineage in a classic sibling species complex. Syst Biol. 2019;68(6):956-966. https://doi.org/10.1093/sysbio/syz027

Macher JN, Rozenberg A, Pauls SU, Tollrian R, Wagner R, Leese F. Assessing the phylogeographic history of the montane caddisfly *Thremma gallicum* using mitochondrial and restriction-site-associated DNA (RAD) markers. Ecol Evol. 2015;5(3):648-662. https://doi.org/10.1002/ece3.1366

Malinsky M, Matschiner M, Svardal H. Dsuite—fast D-statistics and related admixture evidence from VCF files. Mol Ecol Resour. 2021;21(2):584-595. https://doi.org/10.1111/1755-0998.13265

Malinsky M, Svardal H, Tyers AM, Miska EA, Genner MJ, Turner GF, Durbin R. Whole-genome sequences of Malawi cichlids reveal multiple radiations interconnected by gene flow. Nat Ecol Evol. 2018;2(12):1940-1955. https://doi.org/10.1038/s41559-018-0717-x

Mao X, Rossiter SJ. Genome-wide data reveal discordant mitonuclear introgression in the intermediate horseshoe bat (*Rhinolophus affinis*). Mol Phylogenet Evol. 2020;150:106886. https://doi.org/10.1016/j.ympev.2020.106886

Marusik YM, Hippa H, Koponen S. Spiders (Araneae) from the Altai area, Sothern Siberia. Acta Zool Fenn. 1996;201:11-45.

Mirarab S, Warnow T. ASTRAL-II: Coalescent-based species tree estimation with many hundreds of taxa and thousands of genes. Bioinformatics. 2015;31(12):i44-i52. https://doi.org/10.1093/bioinformatics/btv234

Miura O, Urabe M, Mori H, Chiba S. Ancient drainage networks mediated a large-scale genetic introgression in the East Asian freshwater snails. Ecol Evol. 2020;10(15):8186-8196. https://doi.org/10.1002/ece3.6523

Muster C, Berendonk TU. Divergence and diversity: lessons from an arctic-alpine distribution (*Pardosa saltuaria* group, Lycosidae). Mol Ecol. 2006;15(10):2921-2933. https://doi.org/10.1111/j.1365-294X.2006.02989.x

Muster C, Maddison WP, Uhlmann S, Berendonk TU, Vogler AP. Arctic-Alpine distributions—metapopulations on a continental scale? Am Nat. 2009;173(3):313-326. https://doi.org/10.1086/596534

Mutanen M, Hausmann A, Hebert PDN, Landry J-F, de Waard JR, Huemer P. Allopatry as a Gordian Knot for taxonomists: patterns of DNA barcode divergence in Arctic-Alpine Lepidoptera. PLoS One. 2012;7(10):e47214. https://doi.org/10.1371/journal.pone.0047214

Nentwig W, Blick T, Bosmans R, Gloor D, Hänggi A, Kropf C. Spiders of Europe. Version 12.2020; 2020. https://araneae.nmbe.ch/

Nguyen L-T, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32(1): 268-274. https://doi.org/10.1093/molbev/msu300

Obertegger U, Cieplinski A, Fontaneto D, Papakostas S. Mitonuclear discordance as a confounding factor in the DNA taxonomy of monogonont rotifers. Zool Scr. 2018;47(1):122-132. https://doi.org/10.1111/zsc.12264

Oxford GS, Bolzern A. Molecules v. Morphology—is *Eratigena atrica* (Araneae: Agelenidae) one species or three? Arachnology. 2018;17:337-357.

Pan D, Hülber K, Willner W, Schneeweiss GM. An explicit test of Pleistocene survival in peripheral versus nunatak refugia in two high mountain plant species. Mol Ecol. 2020;29(1):172-183. https://doi.org/10.1111/mec.15316

Paradis E. pegas: an R package for population genetics with an integrated-modular approach. Bioinformatics. 2010;26(3):419-420. https://doi.org/10.1093/bioinformatics/btp696

Paradis E, Claude J, Strimmer K. APE: Analyses of phylogenetics and evolution in R language. Bioinformatics. 2004;20(2):289-290. https://doi.org/10.1093/bioinformatics/btg412

Peterson BK, Weber JN, Kay EH, Fisher HS, Hoekstra HE. Double digest RADseq: an inexpensive method for *de novo* SNP discovery and genotyping in model and non-model species. PLoS One. 2012;7(5):e37135. https://doi.org/10.1371/journal.pone.0037135

Piacentini LN, Ramírez MJ. Hunting the wolf: a molecular phylogeny of the wolf spiders (Araneae, Lycosidae). Mol Phylogenet Evol. 2019;136:227-240. https://doi.org/10.1016/j.ympev.2019.04.004

Pittet L, Fragnière Y, Grünig S, Bétrisey S, Clément B, Gerber E, Ronikier M, Kozlowski G, Parisod C. Genetic structure of the endemic *Papaver occidentale* indicates survival and immigration in the Western Prealps. Alp Bot. 2020;130(2):129-140. https://doi.org/10.1007/s00035-020-00238-3

Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945-959. https://doi.org/10.1093/genetics/155.2.945

Procheş S, Ramdhani S. Eighty-three lineages that took over the world: a first review of terrestrial cosmopolitan tetrapods. J Biogeogr. 2013;40(10):1819-1831. https://doi.org/10.1111/jbi.12125

Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior summarization in Bayesian phylogenetics using tracer 1.7. Syst Biol. 2018;67(5):901-904. https://doi.org/10.1093/sysbio/syy032

Ratnasingham S, Hebert PDN. A DNA-based registry for all animal species: the Barcode Index Number (BIN) system. PLoS One. 2013;8(7):e66213. https://doi.org/10.1371/journal.pone.0066213

Rosenberg N. DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4:137-138.

Sahoo RK, Lohman DJ, Wahlberg N, Müller CJ, Brattström O, Collins SC, Peggie D, Aduse-Poku K, Kodandaramaiah U. Evolution of *Hypolimnas* butterflies (Nymphalidae): out-of-Africa origin and *Wolbachia*-mediated introgression. Mol Phylogenet Evol. 2018;123:50-58. https://doi.org/10.1016/j.ympev.2018.02.001
Sanmartín I, Enghoff H, Ronquist F. Patterns of animal dispersal, vicariance and diversification in the Holarctic. *Biol J Linn Soc.* 2001;73:345–390.

Sanz M, Schönswetter P, Vallès J, Vilatersana R. Glacial survival in and recent long-distance dispersal to the Iberian Mountains: the phylogeographic history of *Artemisia umbelliformis* (Asteraceae). *Bot J Linn Soc.* 2017;183(4):587–599. https://doi.org/10.1093/botlinnean/box002

Sawyer YE, Cook JA. Phylogeographic structure in long-tailed voles (Rodentia: Arvicolinae) belies the complex Pleistocene history of isolation, divergence, and recolonization of Northwest North America's fauna. *Ecol Evol.* 2016;6(18):6633–6647. https://doi.org/10.1002/ece3.2393

Schär S, Talavera G, Espadaler X, Rana JD, Andersen Andersen A, Cover SP, Vila R. Do Holarctic ant species exist? Trans-Beringian dispersal and homoplasy in the formicidae. *J Biogeogr.* 2018;45(8):1917–1928. https://doi.org/10.1111/jbi.13380

Schmitt T. Molecular biogeography of Europe: pleistocene cycles and postglacial trends. *Front Zool.* 2007;4(1):11–13. https://doi.org/10.1186/1742-9994-4-11

Schmitt T. Biogeographical and evolutionary importance of the European high mountain systems. *Front Zool.* 2009;6(1):9–10. https://doi.org/10.1186/1742-9994-6-9

Schmitt T, Louy D, Zimmermann E, Habel JC. Species radiation in the Alps: multiple range shifts caused diversification in Ringlet butterflies in the European high mountains. *Org Divers Evol.* 2016;16(4):791–808. https://doi.org/10.1007/s13127-016-0282-6

Schmitt T, Muster C, Schönswetter P. Are disjunct Alpine and Arctic-Alpine animal and plant species in the western Palearctic really “Relics of a Cold Past?” In: Reli. Species. Berlin (Germany): Springer; 2010. p. 239–252.

Solís-Lemus C, Bastide P, Ané C. PhyloNetworks: a package for phylogenetic networks. *Mol Biol Evol.* 2017;34(12):3292–3298. https://doi.org/10.1093/molbev/msx235

Spelda J, Reip H, Oliveira Biener U, Melzer R. Barcoding Fauna Bavarica: Myriapoda—a contribution to DNA sequence-based identifications of centipedes and millipedes (Chilopoda, Diplopoda). *Zookeys.* 2011:156:123–139. https://doi.org/10.3897/zookeys.156.2176

Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. *Bioinformatics.* 2014;30(9):1312–1313. https://doi.org/10.1093/bioinformatics/btu033

Strand E. Neue oder wenig bekannte amerikanische Lycosiden, aus der Sammlung des vestorbenen Mr Thomas Workman. Z Nat. 1909:81:277–286.

Tseng S-P, Wetterer JK, Suarez AV, Lee C-Y, Yoshimura T, Shoemaker D, Yang C-CS. Genetic diversity and *Wolbachia* infection patterns in a globally distributed invasive ant. *Front Genet.* 2019;10:838. https://doi.org/10.3389/fgene.2019.00838

Vogel BR. A review of the spider genera *Pardosa* and *Acantholycosa* (Araneae, Lycosidae) of the 48 contiguous United States. *J. Arachnol.* 2004;32(1):55–108. https://doi.org/10.1636/h03-8

Wagner ND, Gramlich S, Hörandl E. RAD sequencing resolved phylogenetic relationships in European shrub willows (*Salix* L. subg. Chamaetia and subg. Vetrix) and revealed multiple evolution of dwarf shrubs. *Ecol Evol.* 2018;8(16):8243–8255. https://doi.org/10.1002/ece3.4360

Wang J, Ni Y, Hu W, Yin M. Lineage diversity and gene introgression in freshwater cladoceran crustaceans of the *Chydorus sphaericus* species complex. *Limnol Oceanogr.* 2020;66(1):95–107. https://doi.org/10.1002/lno.11590

Wunderlich J. Seltene und bisher unbekannte Wolfspinnen aus Mitteleuropa und Revision der *Pardosa* saltuaria - Gruppe (Arachnida: Araneae: Lycosidae). Ver Nat Ver Hambg. 1984;27:417–442.

Zhang J, Kobert K, Flouri T, Stamatakis A. PEAR: a fast and accurate illumina paired-end reAd mergeR. *Bioinformatics.* 2014;30(5) 614–620. https://doi.org/10.1093/bioinformatics/btt593