Mitochondrial discordance in closely related Theridion spiders (Araneae, Theridiidae), with description of a new species of the T. melanurum group

The incorporation of molecular data into current taxonomic practise has unravelled instances of incongruence among different data sets. Here we report a case of mitochondrial discordance in cobweb spiders of the genus Theridion Walckenaer, 1805 from the Iberian Peninsula. Morphological examination of samples from a country-wide bioinventory initiative revealed the existence of a putative new species and two nominal species belonging to the Theridion melanurum species group. The morphological delineation was supported by the molecular analysis of a nuclear marker but was at odds with the groups circumscribed by a mitochondrial marker. The causes of this discordance remained uncertain, once sample and sequencing errors and the existence of pseudogenes were discarded. The full sorting observed in the alleles of the more slowly evolving nuclear marker ruled out incomplete lineage sorting, while the geographic patterns recovered were difficult to reconciliate with ongoing hybridization. We propose that the apparent incongruence observed is most likely the result of old introgression events in a group with high dispersal abilities. We further speculate that endosymbiont-driven cytoplasmatic incompatibility could be involved in the fixation of mitochondrial haplotypes across species barriers. Additionally, we describe the new species T. promiscuum sp. nov., based on the presence of diagnostic morphological traits, backed up by the nuclear data delimitation. Our study contributes yet another example of the perils of relying on single methods or data sources to summarise the variation generated by multiple processes acting through thousands of years of evolution and supports the key role of biological inventories in improving our knowledge of invertebrate biodiversity.


Introduction
The incorporation of DNA sequence information in species delimitation and description has become a gold standard in current taxonomic practise. Under an integrative taxonomy framework (Dayrat 2005;Will et al. 2005;Padial et al. 2010) the use of DNA sequence data allows for a finer species delimitation, reveals population structure, simplifies matching life stages or remnants, and facilitates identification for taxa otherwise difficult to diagnose using only morphological charac-ters (e.g. cryptic species) (Pires and Marinoni 2010). The comparative analyses of DNA sequences both from different markers and with phenotypic data have unravelled instances of incongruence between data sets, which may point to relevant biological processes in the origin and maintenance of biological diversity (e.g. Sota and Vogler 2001). At the same time, the reported incongruence may challenge some automated approximations to species delimitation and identification, especially when based on single markers (i.e. DNA barcoding) .
Hybridization and incomplete lineage sorting are among the main processes accounting for the apparent incongruence between different sources of taxonomic evidence, namely mitochondrial DNA, nuclear DNA, or the phenotype. Hybridization is defined as the interbreeding of individuals from different species. When this hybridization process involves repeated backcrossing of the hybrids with the parent species, we talk about introgression. The notion that species boundaries are not impermeable but can be porous to introgression is becoming more and more accepted as new evidence appears (Harrison and Larson 2014). The most iconic case might be the recent evidence provided by ancient DNA of Neanderthal introgression with modern humans (Green et al. 2010). Multiple instances of introgression have also been reported within arthropods. Spiders have not been spared from undergoing hybridization and introgression processes. To cite some examples, several studies have revealed hybrid zones among common funnel weaver spiders of the Eratigena atrica (C.L. Koch, 1843) group (Agelenidae) in Britain and Continental Europe (Croucher et al. 2007;Oxford and Bolzern 2018), and in deeply divergent lineages of the ground-dweller Harpactocrates Simon, 1914 spiders in the Pyrenees (Bidegaray- Batista et al. 2016). In New Zealand, Vink et al. (2008) identified a case of introgression from the Australian black widow Latrodectus hasseltii Thorell, 1870 (Theridiidae) into the endemic L. katipo Powell, 1871, while Lattimore et al. (2011) found evidence of unidirectional introgression within the fishing spider genus Dolomedes Latreille, 1804 (Pisauridae). More recently, Leduc-Robert and Maddison (2018) have revealed that hybridization may have been widespread across phylogenetically distant species of Habronattus F.O. Pickard-Cambridge, 1901 jumping spiders (Salticidae).
Introgression between species can be detected by comparing information from mitochondrial (mtDNA) and nuclear DNA (nucDNA). As mtDNA is maternally inherited (Giles et al. 1980), when two individuals of different species hybridize, the offspring's mtDNA will be 100% identical to that of the maternal species, but the nucDNA will be 50% of each parental species. In successive generations, the nucDNA will homogenize through recombination within the population, while the introgressed mtDNA will remain traceable in the population in following generations, unless lost by genetic drift. The hybridization event will be detectable because some individuals will possess the nucDNA and morphology from the paternal population but the mtDNA from the maternal one.
Another source of incongruence among multiple lines of taxonomic evidence in closely related species is incomplete lineage sorting (ILS). ILS is the process by which, as a result of the segregation of an ancestral polymorphism, i.e. the existence of two or more homologous alleles predating the speciation event, the evolutionary relationships between individuals given by the sequences of a certain gene do not match the species phylogeny (Rogers and Gibbs 2014). The occurrence of ILS is especially pervasive in either recently diverged species or in species with large ancestral effective population sizes (Maddison 1997). Under an ILS scenario, when comparing mtDNA and nucDNA from various species, we would expect mitochondrial gene trees to reflect better the species tree than nuclear gene trees, because, other factors being equal, mtDNA has a smaller effective population size in diploid organisms and is known to have a faster mutation rate than nucDNA (Ballard and Whitlock 2004).
The family Theridiidae is one of the richest and most ecologically diverse spider families, containing 2,516 species grouped in 124 genera (World Spider Catalog 2020). The family is reputed for including, among others, species of medical importance, such as the widow spiders (Latrodectus Walckenaer, 1805), or some of the few examples of sociability within spiders (e.g. Anelosimus Simon, 1891). Theridiids usually construct tangle webs with gumfoots, i.e. sticky droplets on silk threads radiating from mesh retreats or web hubs attached to the substrate, mainly aimed at capturing pedestrian prey (Blackledge et al. 2011). They are occasional to frequent ballooners, and exhibit a great diversity in morphology, ecology, and behaviour, ranging from solitary web-less hunters, such as the species of the genus Euryopis Menge, 1868 (Levi 1954) to elaborate tangled, three-dimensional space webs or even to large colonies of thousands of individuals, like the social spider Anelosimus eximius (Keyserling, 1884).
The genus Theridion Walckenaer, 1805 contains 584 described species, more than one-fifth of all the theridiids found worldwide (World Spider Catalog 2020). This remarkable species diversity is most likely a taxonomic artifact because the genus has been traditionally used as a dumping ground for theridiids with no trace of colulus that do not fit into better diagnosed genera (Forster et al. 1990). Because of its large size, several smaller species groups have been proposed to further divide Theridion (e.g. Levi 1957), some of which have been subsequently erected to genus status. One of these species' complexes, the Theridion melanurum species-group (Wunderlich 2011), contains 10 very similar species, often difficult to tell apart. The group contains, along Theridion melanurum Hahn, 1831, the species T. asopi Vanuytven, 2014, T. bernardi Lecigne, 2017, T. betteni Wiehle, 1960, T. bosniense Wunderlich, 2011, T. cairoense Wunderlich, 2011, T. harmsi Wunderlich, 2011, T. malagaense Wunderlich, 2011, T. musivivum Schmidt, 1956, T. mystaceum L. Koch, 1870, and T. semitinctum Simon 1914 In the general framework of a biological inventory of the spider communities of white oak woodlands of the Spanish National Park Network (Crespo et al. 2018), we had the opportunity to examine a large sample of Theridion specimens and obtained DNA barcoding data for many of these individuals. Comparison between identifications based on morphology and those based on DNA barcoding revealed several instances of incongruence that were subsequently examined and corroborated using the nuclear intron Internal Transcribed Spacer 2 (ITS-2). Here we present the results of the morphological and molecular analyses, describe a new species based on the combination of the available evidence, and further discuss the implications of these results for our understanding of the diversification of species-rich spider lineages.

Material and methods
Specimens were collected using semi-quantitative methods as part of a larger project that aimed to understand the diversity patterns of Iberian spider communities (Crespo et al. 2018). The sampling design included 16 × 1 ha plots distributed in white oak forests across six national parks of the Iberian Peninsula. Sampling was conducted using the COBRA protocol (Cardoso 2009) in May and June of 2013 and 2014. The sampling protocol combined timed direct capture, beating, and sweeping with 48 pitfall traps per plot, which were active for two weeks. The distribution of the sampling sites is shown in Figure 1. Detailed collection and genetic information of all specimens used in this work can be found in Suppl. material 1.
Specimens were sorted and identified under a ZEISS Stemi 2000 stereomicroscope. We took photographs using a Leica DFC 450 camera attached to a Leica MZ 16A stereomicroscope, with the software Leica Application Suite v. 4.4. Both the male palp and the female epigyne were excised with the help of entomological needles to facilitate observation under the scope. The muscle tissue of the epigyne was further removed with the needles and digested using immersion in potassium hydroxide (KOH) at a 30% concentration. For SEM examination, palps were excised and cleaned ultrasonically for 1 min and then transferred to 100% ethanol overnight. Palps were submitted to critical point drying, glued to flat-headed rivets and gold sputter coated. Imaging was conducted with the help of a Quanta 200 environmental SEM. We used the Araneae: Spiders of Europe online identification tool to identify most of the species found in our samples (Nentwig et al. 2019). Type specimens were deposited at Museu de Ciències Naturals de Barcelona, Spain (CMCNB) and additional voucher specimens at the Centre de Recursos de Biodiversitat Animal, University of Barcelona (CRBA). Description of the new species followed the format of Vanuytven (2014). All measurements are in millimetres.
For molecular analyses, we aimed at including both sexes and all the sites were each species was collected. We extracted DNA from two to four legs from each specimen using REDExtract-N-Amp™ Tissue PCR Kit Protocol from Sigma-Aldrich, following the manufacturer's protocol, performed in 96-well plates. We amplified fragments of the animal DNA barcode gene cytochrome c oxidase subunit I (COI) and the nuclear Internal Transcribed Spacer 2 (ITS2). Primers used for amplification are shown in Table 1. We performed polymerase chain reactions in 96-well plates using 8 µL REDExtract-N-Amp™ PCR ReadyMix from Sigma-Aldrich, primers forward and reverse, 4 µL of diluted DNA and ultrapure, distilled water up to a total reaction volume of 20 µL. PCR condi-tions were as follows: initial denaturing step at 95 °C for 5 min, 35 amplification cycles (94 °C for 30 s, 45 °C for 35 s, 72 °C for 45 s) and a final step at 72 °C for 5 min. PCR products were cycle-sequenced in both directions with the sample amplification primers at Macrogen Inc. (Seoul, South Korea).
We edited and manipulated all sequences using Geneious v. 10.2 (Kearse et al. 2012). For the COI sequences alignment was trivial since there was no evidence of insertions or deletions. Conversely, fragments of the ITS2 sequences showed between and within individual length polymorphisms. Single individual length polymorphisms were manually phased by editing original chromatograms using a sliding window approach. Phased ITS2 sequences were aligned using the automatic algorithm implemented in the online version v. 7.017 of MAFFT (Katoh and Standley 2013) using the G-INS-i strategy. Gaps in the alignment were treated as separate absence/presence characters in subsequent analyses using FastGap v. 1.2 (Borchsenius 2009).
We inferred a Maximum Likelihood tree for the COI with the program RAxML (Stamatakis 2014), using 20 runs of random addition of taxa, using only unique sequences and assessing support by means of MRE bootstrapping. We used the best codon partition and corresponding model retrieved by PartitionFinder v. 2.1 (Lanfear et al. 2012) using the AICc criterion, which was the general time-reversible nucleotide substitution model with Gamma-distributed among-site rate variation and invariant sites (GTR+G+I), and with a separate partition for each codon position. Additionally, we constructed haplotype and allele networks for the COI and ITS2, respectively, using the TCS algorithm (Clement et al. 2002) as implemented in the program PopArt v. 1.7 (Leigh and Bryant 2015). We calculated intraspecific and interspecific p-distances for the COI and the ITS2 for all species using MEGA v. 7.0.26 (Kumar et al. 2016).
We conducted additional species delineation using the barcode index number (BIN) method (Ratnasingham and Hebert 2013) (2 males). In addition, we found two morphotypes that did not fit into any of the nominal species descriptions, tentatively referred to as T. sp. 6 (seven males, five females) and T. sp. 15 (one male). The former was found to be a new species to science, hereby described.

Molecular analyses
We obtained the COI sequences of 117 specimens ( The ML tree obtained for the COI (Fig. 2) supported the close relationships between the haplotypes of the T. melanurum group but failed to recover monophyly of the different species. Of the four main lineages recovered, albeit some of them with low support, one included all the haplotypes of T. harmsi, along with single haplotypes of T. mystaceum and T. sp. 6, the second included sequences of T. mystaceum, T. sp. 6, and T. betteni, the third two divergent haplotypes of T. mystaceum and the fourth haplotypes of T. mystaceum, T. sp. 6, and T. melanurum. The average uncorrected genetic distances within the T. melanurum group was 2.8% (SE 0.004). Each of the remaining Theridion species sampled were recovered as monophyletic and clearly divergent from the haplotypes of the T. melanurum group (9% or above average uncorrected genetic distance), including the single haplotype of the morphotype T. sp. 15.
The COI sequences were split into nine different BINs, one unique and the remaining including additional specimen available in BOLD. The species delineation of the T. melanurum complex yielded either BINs with mixed species (Fig. 2) or BINs with single specimens, in all cases identified as T. mystaceum. Conversely, in the remaining Theridion species, the BINs mirrored the morphological identifications.
The COI haplotype network ( Fig. 3) was largely congruent with the ML tree but provided a more detailed view on the relationships and geographic distribution of the T. melanurum group haplotypes. Neither a taxonomic nor a geographic structure was readily apparent in the T. melanurum group network. All T. harmsi haplotypes were highly similar, separated by two or less mutations, and the most common haplotype was widespread all over the sampled parks. The closest haplotype to this group was found in individuals belonging to T. mystaceum and T. sp. 6 from different parks. The bulk of T. mystaceum haplotypes were separated by five or less mutations, and also included haplotypes found in T. sp. 6, closely related to a haplotype of T. mystaceum form the same park. One of the haplotypes of this cluster was shared between two parks but also found in Germany and in a population of T. betteni from Britain. A second cluster of T. mystaceum haplotypes, separated by six or less mutations, included a haplotype shared with T. sp. 6 from a different park, and a very closely related haplotype found in T. melanurum from Switzerland.
In contrast to the mitochondrial network, the nuclear ITS2 network perfectly matched the morphology-based species delimitations (Fig. 4). The alleles sampled from each species were exclusive and each other's closest relatives. The highest allele diversity was found in T. harmsi, whose most common alleles were found across most of the parks. The remaining species yielded two or fewer alleles.
The ITS2 maximum uncorrected intraspecific divergences ranged between 0 and 4.5% (Table 2). Minimum interspecific divergences within the species complex ranged from 2.6% (between T. mystaceum and T. sp. 6) to 4.3% (between T. harmsi and T. sp. 6), as seen in Table  3. These values fall within or above the ranges observed within and among closely related spider species across several families (see Agnarsson 2010).     The combination of diagnostic morphological features both in males and females, similar to those found among other species within the T. melanaurum group, and the additional support of the nuclear data delimitation led us to propose that T. sp. 6 morphotype actually corresponds to a new species. Etymology. The specific name promiscuum is derived from the Latin word for "intermingling", referring to the morphological similarity with other species of the T. melanurum group. It also makes reference to the possible introgression between this species and other closely related species here reported.

Diagnosis.
Males of T. promiscuum sp. nov. can be distinguished from other species in T. melanurum group by a bent and twisted embolus, forming a marked angle (Fig. 5e); embolus shorter in T. promiscuum sp. nov. than in sympatric T. mystaceum (Fig. 6e) and T. harmsi (Fig. 7e); tip of the embolus pointing distally in T. promiscuum sp. nov. (Fig. 5c) and pointing retrolaterally in remaining species (Figs 6c, 7c); base, prolateral part of the embolus more sclerotized in T. promiscuum sp. nov. (Fig. 5e) than in remaining species of the T. melanurum group (Figs 6e, 7e); presence of small barbs in embolus base distal border is unique to T. promiscuum sp. nov. (only visible in SEM; Fig. 8a, b). Females of T. promiscuum sp. nov. can be distinguished from T. mystaceum and T. harmsi by the lower position of copulatory duct openings in atrium in T. promiscuum sp. nov. (Fig. 5f), which are central in T. harmsi (Fig. 7f) and close to or under lateral walls of atrium in T. mystaceum (Fig. 6f); copulatory ducts shorter, fewer coils in T. promiscuum sp. nov. (Fig.  5g) than in T. mystaceum (Fig. 6g) and T. harmsi (Fig. 7g).

Genetic information.
We obtained a 568 bp sequence of the COI mitochondrial gene for the holotype (stored in Genbank under the code MT215600) and five paratypes (MZB 2017-3710, MZB 2017-3712, MZB 2017-3713, MZB 2017-3714, and MZB 2017 with GenBank codes MT215603, MT215604, MT215602, MT215606 and MT215601, respectively). We also obtained sequences of the ITS2 of up to 469 bp, including up to 26 bp of the 28S gene and up to 69 bp of the 5.8S gene, for the holotype (GenBank code MT117179) and the five paratypes mentioned before (codes MT117182, MT117181, MT117180, MT117183 and MT117184). Distribution and natural history. The holotype and paratypes were collected in an oak forest of Quercus pubescens Willd. in Aigüestortes i Estany de Sant Maurici National Park, located in the southern slopes of the Catalan Pyrenees. The specimens were captured either by beating, sweeping or direct sampling, but not by pitfall traps, which suggests that this species is found in the vegetation at a certain height above the ground. Only the three specimens captured by beating method were captured during the day, whereas the sweeping and direct sampling ones were captured at night, which indicates that this species is mostly active at night.
Other taxa analysed. The specimen tentatively referred to as T. sp. 15 is morphologically similar to the species T. cinereum Thorell, 1875, T. petraeum L. Koch, 1872, T. furfuraceum Simon, 1914, T. pyrenaeum Denis, 1944, and T. wiehlei Schenkel, 1938. However, the male palp of the single specimen available showed slight differences in the shape of the median apophysis, slightly tilted, and the more pronounced basal curvature of the shorter embolus, which refrained us from assigning it to any of the former species. The DNA barcodes, on the other hand, identified the specimen unambiguously as T. cinereum.

Discussion
The genus Theridion is already one of the largest genera within spiders, yet new species continue to be described yearly, with no evidence of reaching a plateau. Surprisingly, even the well-known European fauna is still contributing new Theridion species. The newest, T. bernardi, also belonging to the T. melanurum group, was described from Portugal in 2017 (Lecigne 2017). From this prospective, the discovery of yet another new species of Theridion from the T. melanurum group does not greatly improve our knowledge of this diverse yet poorly defined genus. However, the identification of a clear incongruence between some of the sources providing evidence for the evolutionary independence of this new species presents some challenges and opens interesting evolutionary questions. This study contributes yet another example of the perils of relying on single methods or data sources to summarise the variation generated by multiple processes acting through thousands of years of evolution. Operational approaches such as DNA barcoding may greatly help to accelerate biological inventories and provide useful information for refining species delimitations and speed up identification. However, the multiple factors involved in the generation and evolution of species -e.g. speciation mechanisms, species ontogeny, gene flow -require the integration of multiple sources of evidence to provide a sound hypothesis of species.  Our data clearly shows that while the phenotype and the nuclear data nicely delimit species boundaries in the studied species, the mitochondrial data suggest mixture across members of the T. melanurum group. The most obvious explanation for the apparent incongruence would be the human error in the manipulation of specimens or in the lab procedures. We discarded this source of error by re-extracting and resequencing specimens involved at least twice, independently. On the other hand, it is well known that mitochondrial pseudogenes, either as duplications within the mitochondrial genome or as copies inserted into the nuclear genome (e.g. NUMTS) may compromise the recovery of species boundaries using DNA barcoding approaches (Song et al. 2008) and, therefore, be in conflict with other sources of evidence. We examined all mitochondrial sequences generated in our study in search for some evidence of the presence of pseudogenes or NUMTs in our data set, conducting a careful examination of sequence characteristics, including indels, in-frame stop codons and nucleotide composition. We found neither changes in reading frames nor in the base composition-all sequences showed the characteristic AT bias observed in many arthropod mitogenomes and all passed the composition chi-squared test.
Among the biological processes that may account for inconsistencies among different molecular markers, the incomplete sorting of ancestral polymorphism seems unlikely in our case. At least in spiders, the higher mutation rates and smaller population sizes of mitochondrial markers should ensure sorting in those markers previous to the nuclear ones, an opposite pattern of the one here reported. Although slower evolutionary rates in mitochondrial DNA have been found in some tetrapods (Nabholz et al. 2009) and cnidarians (Shearer et al. 2002), they have never been reported in spiders. Hybridization is the main biological culprit for mitochondrial discordance. However, the individuals sharing the same mitochondrial haplotype in our study were located hundreds of kilometres apart, which casts some doubts on ongoing gene flow as a main explanation for the incongruence pattern recovered. It could be argued, however, that, like many other spiders, Theridion is able to drift long distances by ballooning, i.e. silk mediated airborne dispersal, which partially explains the large distribution ranges of many species. However, the extremely high dispersal rates required to explain the reported patterns of mixture in mtDNA should most likely also affect the nuclear and phenotypic data, which does not seem to be the case. Alternatively, the present day mixing of haplotypes within the species of the T. melanurum complex investigated could have been the result of older, hybridization events, possibly at a younger stage of the species ontogeny. The high dispersal rates of the species would have subsequently blurred the geographical signal of the hybridization processes.
Although at the moment we lack any solid evidence, there are also some chances that the mitochondrial incongruence observed could have been maintained by the involvement of endosymbiont driven cytoplasmic incompatibility. It is well known that endosymbiotic bacteria may affect the patterns of mitochondrial variation in invertebrates, which may compromise the inferences made on host evolution from these patterns (Hurst and Jiggins 2005). Although less well characterised than in insects, spiders are also known to host endosymbiotic bacteria (Goodacre et al. 2006;Zhang et al. 2018). The infection frequency may be as high as seventy percent of the individuals and specimens may be infected by multiple bacteria species and strain types (White et al. 2019). Studies on the diversity and structuring of Wolbachia endosymbionts in the funnel web spider genus Agelenopsis Giebel, 1869 suggest that Wolbachia strains induce cytoplasmic incompatibility, which may explain fixation of mtDNA haplotypes across wide geographic distances in these spiders (Baldo et al. 2008). Similarly, endosymbiont driven fixation of mtDNA haplotypes has been invoked to explain the mito-nuclear discordance recently reported in wolf spider species (Ivanov et al. 2018). We propose that cytoplasmatic incompatibility driven by endosymbionts could have fixed mitochondrial haplotype mixtures across species in the T. melanurum complex originated from ancestral gene flow or ILS. Genotyping the endosymbiont strains hosted by those individuals from different species bearing similar haplotypes will have to be conducted to test our hypothesis.
The genus Theridion is a classic example of catch-all genus, a poorly defined group to include species with no trace of colulus that do not match other more precisely diagnosed genera. The redefinition of the genus would require an exhaustive systematic revision of a thorough sample of species, which is beyond the scope of the present paper. However, our results do support the existence of complexes of species within the genus. Specifically, molecular data (Fig. 2) confirms the delimitation of the T. melanurum species group, given the short branch lengths within it, in comparison to the branches that separate the remaining species. The Theridion species that were included in the tree and that were previously thought to be-long to this species complex (T. mystaceum, T. harmsi, T. melanurum, and T. betteni) are confirmed to be part of it. Additionally, the new species T. promiscuum sp. nov. is supported as a new member of the T. melanurum group. Species complexes, i.e. closely related species with limited morphological variation, are ideally suited to unravel the processes underlying speciation and the factors driving phenotypic change. Assuming that similarity is the result of recent divergence, there are better chances that the patterns of variation exhibited by species in those complexes are directly linked to the causal agents of the diversification of the group, and that they have not been blurred by subsequent processes. Future studies may take advantage of the melanurum and other species complexes within the species rich genus Theridion to gain better insights on the origins and maintenance of spider diversity in an age of major biodiversity loss.

Conclusions
We investigated here one of the few cases reported of mitochondrial discordance within spiders. A countrywide bioinventorying project revealed the existence of a new, morphologically diagnosable spider species within the T. melanurum group. Subsequent mtDNA barcode screening of specimens, however, identified instances of haplotype mixing across closely related species. Molecular information from a nuclear marker, on the other hand, supported the morphological delimitations, including the new specific status. The lack of geographic structure in the shared haplotypes and the lack of sorting in the fastest evolving gene suggests that mechanisms other than ongoing gene flow and deep coalescence are responsible for the observed patterns. We propose that cytoplasmatic incompatibility mediated by endosymbionts may have been instrumental in generating mito-nuclear discordance, probably originated from old introgression events. Finally, this study highlights the important role that bioinventories play in improving our knowledge of biodiversity, especially in a time when fieldwork studies that gather new data are becoming less popular than those using pre-existing data, like modelling or meta-analysis studies (Ríos-Saldaña et al. 2018). ra Baldo kindly contributed insights into the interpretation of the results. We thank Cor Vik and Nadine Dupérré for their comments. This work was funded by the Organismo Autónomo de Parques Nacionales (OAPN #485/2012). M.D. thanks the University of Barcelona for the APIF PhD fellowship supporting this research. Additional support was provided by 2017SGR73 from the Catalan Government.