Porous barriers? Assessment of gene flow within and among sympatric long‐eared bat species

Abstract Species are the basic units for measuring biodiversity and for comprehending biological interactions. Yet, their delineation is often contentious, especially in groups that are both diverse and phenotypically conservative. Three cryptic species of long‐eared bats, Plecotus auritus, P. austriacus, and P. macrobullaris, co‐occur over extensive areas of Western Europe. The latter is a fairly recent discovery, questioning the overall diversity of the entire Plecotus complex. Yet, high morphological and acoustic similarities compromise the reliable identification of long‐eared bats in the field. We postulate that such extensive phenotypic overlap, along with the recurrent observation of morphologically intermediate individuals, may hide rampant interspecific hybridization. Based on a geographic sampling centered on areas of sympatry in the Alps and Corsica, we assessed the level of reproductive isolation of these three Plecotus species with mitochondrial and nuclear markers, looking at both inter‐ and intraspecific genetic population structuring. No sign of hybridization was detected between these three species that appear well separated biologically. Genetic structuring of populations, however, reflected different species‐specific responses to environmental connectivity, that is, to the presence of orographic or sea barriers. While the Alpine range and the Ligurian Sea coincided with sharp genetic discontinuities in P. macrobullaris and P. austriacus, the more ubiquitous P. auritus showed no significant population structuration. There were clear phylogeographic discrepancies between microsatellite and mitochondrial markers at the intraspecific level, however, which challenges the reliance on simple barcoding approaches for the delineation of sound conservation units.

Beyond reproductive barriers, species are also confronted to topographical barriers that can influence the genetic architecture of populations. In particular, major landscape features such as mountain ranges or large expanses of open water may have served as important barriers to dispersal in bats (Castella, Ruedi, & Excoffier, 2001;Dàvalos, 2005;Juste et al., 2004). This is particularly true in European biotas which experienced multiple events of range contraction and expansion caused by glacial cycles (Hewitt, 1999). Southern glacial refugia, such as the Iberian, Apennine, and Balkan peninsulas, retained genetic diversity and allowed population differentiation during phases of allopatry, whereas topographical barriers such as the Alps or the Pyrenees constrained recolonization routes from those refugia (Çoraman, Furman, Karataş, & Bilgin, 2013;Hewitt, 2000). Predicting the effect of these barriers on gene flow can be tricky, since co-distributed species with distinct ecological needs may show different responses to the same landscape elements, depending on their ability to cross those (Engler, Balkenhol, Filz, Habel, & Rödder, 2014;Zancolli, Rödel, Steffan-Dewenter, & Storfer, 2014). For instance, large expanses of water can limit drastically the dispersal of some bat species, while they have minor effects on others (García-Mudarra, Ibañez, & Juste, 2009). Likewise, the Alps have been shown to delimit major genetic pools in certain bat species (Ruedi et al., 2008;Wright et al., 2018) when others disperse extensively through this range (Rebelo et al., 2012).

| Model species
Delimitation of species boundaries in long-eared bats of the genus Plecotus has been challenging for decades. Because of their very similar external morphology, a single species was considered to occur in Europe until a second one was raised to species level by Bauer (1960) and Lanza (1960). More recently, with the use of genetics, no less than three species were further added to the European fauna Mucedda, Kiefer, Pidinchedda, & Veith, 2002;Spitzenberger, Pialek, & Haring, 2001).
In the Alps, where the three species P. auritus, P. austriacus and P. macrobullaris coexist, their discrimination in the field on the basis of external morphology remains difficult (Ashrafi, Bontadina, Kiefer, Pavlinić, & Arlettaz, 2010). Even the examination of cranial characters showed the existence of intermediate specimens with conflicting genetic and morphologic diagnoses . These discrepancies raise the question whether these three Plecotus species are reproductively isolated from each other or not. If interspecific barriers are porous, these inconsistencies could result from hybridization leading to individuals with atypical morphology or introgressed genes. Furthermore, since the degree of reproductive isolation is expected to increase with time since divergence (Edmands, 2002), we hypothesized that hybridization should be more likely between the two more closely related species (P. auritus and P. macrobullaris) than with P. austriacus, which originated from a much older divergence event (Spitzenberger, Strelkov, Winkler, & Haring, 2006).
Since two or more of these morphologically similar species of long-eared bats occur in sympatry over extensive areas in continental Europe and Corsica Courtois, Rist, & Beuneux, 2011;Gilliéron, Schönbächler, Rochet, & Ruedi, 2015;Rutishauser, Bontadina, Braunisch, Ashrafi, & Arlettaz, 2012;Tvrtković, Pavlinić, & Haring, 2005), they also offer the opportunity to compare the barrier effect of major landscape features on gene flow in each species in parallel. Indeed this geographical region includes two main potential topographic barriers to dispersal and gene flow, the Alps mountain range and the Ligurian Sea. Judging from their similar wing morphology, all three species of long-eared bats are considered as poor long-distance flyers (Entwistle, Racey, & Speakman, 1996). Ringing studies indeed show that they are highly sedentary species (Hutterer, Ivanova, Meyer-Cords, & Rodrigues, 2005;Racey & Entwistle, 2003), although little is known about dispersal abilities in P. macrobullaris (Alberdi, Aihartza, et al., 2015a;Ashrafi et al., 2013). Accordingly, genetic structure of populations in these bats is generally pronounced and potentially highly affected by extrinsic factors such as geographic barriers to dispersal (Burland & Worthington-Wilmer, 2001;Racey & Entwistle, 2003).
However, since the general ecology of these three species differs, for instance in terms of altitudinal preferences (Rutishauser et al., 2012), more specific predictions of the effect of topographic barriers can be formulated. P. macrobullaris is a typically alpine-adapted species that is able to breed at altitudes of 2,000 m a.s.l. and higher (Alberdi, Garin, Aizpurua, & Aihartza, 2013;Anonymous, 2014b); the Alps should likely not be a major obstacle to dispersal, while the species might be more reluctant to cross intervening lowland areas. On the contrary, P. austriacus is a typical lowland species (Anonymous, 2014a;Arthur & Lemaire, 2015;Juste et al., 2008), and the reverse situation may likely prevail. Finally, P. auritus occurs across a wide altitudinal range in the Alps (Anonymous, 2014c;Hutson et al., 2008) and thus is expected to be somewhat intermediate regarding the effect of mountain barriers on gene flow. Predictions about the effect of sea channels on gene flow in these three species are more difficult to draw. On the one hand, P. austriacus seems little affected by sea barriers as it is found on most Mediterranean islands and as far as Madeira (Juste et al., 2004;Spitzenberger et al., 2006), whereas P. auritus or P. macrobullaris appear more reluctant to cross open water as they are found only on one or two of these islands Piraccini, 2016). We can thus predict that the Ligurian Sea should be a more effective barrier to gene flow for the latter two than for the former species.
We used a combination of nuclear and mtDNA markers on bats sampled in areas where those species co-occur in strict sympatry to explore (a) their degree of reproductive isolation and (b) to which extent the Alpine range and Ligurian Sea constitute topographic barriers limiting gene flow among populations.

| Sampling design and assignation to major mtDNA lineages
To detect the presence of putative hybrids between the three Plecotus species, we sampled several areas where they have true chances to interbreed, that is, zones of strict sympatry. A total of 349 individuals representing the three target species were gathered in the western parts of the Alps and in Corsica (Supporting information Appendix S1; Supporting Data 1). Ninety-five samples were issued from museum collections, the remaining 254 ones consisted of wing punch biopsies from individuals captured in the field. Samples originated from France (n = 90), Germany (n = 1), Italia (n = 27), and Switzerland (n = 231) and included bats caught in maternity roosts, in hunting grounds, or in transit zones. Only adult bats were genotyped in order to minimize autocorrelation between samples that could arise from mother and pup pairs. All tissues were stored in ethanol at −20°C. Genomic DNA was isolated with the DNeasy Blood & Tissue Kit (Qiagen, Switzerland). As a preliminary taxonomic assignation, individuals were sorted according to their mtDNA lineages, using an array of methods based on short fragments of the 16S gene, as described in Andriollo and Ruedi (2018). This initial mtDNA assignment resulted in 152 samples carrying mitochondrial haplotypes typical of P. auritus, 79 of P. austriacus, and 118 of P. macrobullaris (Supporting Data 1).

| Amplification and genotyping of nuclear markers
Individual genotypes were further characterized by 23 autosomal microsatellite loci developed by Burland, Barratt, and Racey (1998) and Razgour et al. (2013). Five multiplex PCR amplifications (Supporting information Appendix S2) were optimized to be performed in a 10-μl reaction volume containing 2-10 ng of genomic DNA, 5 μl HotStarTaq Master Mix, 0.3 μM of forward and reverse primers each, and completed with double distilled water. We used the following cycling protocol on a TC-412 Programmable Thermal Controller (Techne): 35 cycles with 94°C for 30 s, 56°C for 90 s, and 72°C for 60 s. Before the first cycle, a prolonged denaturation step (95°C for 15 min) was included and the last cycle was followed by a 30-min extension at 72°C. Amplicons were sized through fragment analysis by a commercial company (Ecogenics, Switzerland) on an ABI 3730 DNA analyser (Applied BioSystems). Alleles were then scored semi-automatically using the microsatellite plugin in Geneious R10 (Kearse et al., 2012). Twenty-eight individuals (8% of the dataset) were amplified and genotyped independently twice for consistency checks, but no discrepancies were observed in scoring these duplicates. We used the package poppr (Kamvar, Tabima, & Grünwald, 2014) in the R environment (R Core Team, 2017) to (a) discard individuals with more than 25% missing genotypes, (b) discard loci with more than 40% missing data in each intraspecific subset, (c) test for Hardy-Weinberg equilibrium (hereafter HWE), and (d) calculate genetic diversity indices. Since our sampling was geographically uneven, some isolated individuals could not be aggregated into biologically meaningful populations for HWE tests.
These statistical tests were thus restricted to the better sampled populations, and loci with consistent departure across these populations (i.e., showing sign of null alleles) were discarded from further analysis. The software coancestry (Wang, 2011) was used to estimate the coefficients of identity Δ 7 and Δ 8 described by Jacquard (1972) through the triadic likelihood method (Wang, 2007). These inbreeding coefficients can identify putative monozygotic twins (Δ 7 = 1; Δ 8 = 0) or parent-offspring pairs (Δ 7 = 0; Δ 8 = 1) that could bias allelic frequency estimations and cause departure from expectations inherent to population genetics models (Wang, 2011).

| Interspecific analyses
In order to quantify levels of interspecific admixture, we analyzed the overall dataset with structure v. 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). This Bayesian clustering method uses a MCMC algorithm to assign individuals in a preset number of groups (K) by minimizing departure from HWE and linkage disequilibrium (LD), with no prior about individual group-membership. Ten independent chains, consisting of 100,000 generations after 50,000 generations of burn-in, were run for each K ranging from 1 to 10. A model with independent allele frequencies among populations was selected.
The assumptions of population equilibrium models (HWE and LD) on which Bayesian clustering approaches rely (Pritchard et al., 2000) are often violated with uneven sample sizes and may bias assignments (Puechmaille, 2016). Since uneven sampling prevailed in our study, we also employed multivariate analyses that are free from these assumptions to describe the diversity of individual genotypes. These methods also perform better in detecting genetic clines as they are not forcing individual assignment into discrete clusters (Patterson, Price, & Reich, 2006). A principal component analysis (PCA) implemented in adegenet (Jombart, 2008) was carried out on the scaled allele frequencies, with missing data (<6% of genotypes) replaced by mean allele frequencies. When applicable, the elbow criterion (Ketchen & Shook, 1996) was used on the scree plot of eigenvalues to characterize the number of principal components best explaining the structure of the data.

| Detection of hybrids
In order to estimate the detection power of interspecific hybrids provided by our multilocus genotyping approach, we applied maximum-likelihood methods to estimate the ancestry index S (the amount of genetic information inherited from a given parent) and the interclass heterozygosity index H I (the proportion of loci with one allele of each parental species) on the basis of prior estimates of parental allele frequencies using the R package HIest (Fitzpatrick, 2012). We calculated parental allele frequencies on the basis of individuals classified in the interspecific structure analysis with Qvalues higher than 0.96, which were considered to represent pure parental individuals. Hybrid genotypes were then generated in silico from these parental pools by simulating random interspecific mating using the package adegenet. S and H I were estimated a posteriori for parental individuals, for simulated offspring and for individuals with uncertain parentage (i.e., those with Q-values < 0.96, which were not attributed to parental populations). Since this descriptive method does not try to sort individuals into predefined classes, the power analysis relied on critical evaluation of S and H I values overlap between different hybrid categories.

| Intraspecific analyses
We ran structure analyses for the three intraspecific datasets as described above, but using the correlated model for allele frequencies.
Since the second-order statistics ΔK cannot correctly identify the best K when K = 1, the Evanno's criterion was used in combination with the maximal logarithm probability of the data when no structure was suspected. PCA were also performed for each species independently. Intraspecific analyses of molecular variance (AMOVA) were performed with arlequin, version 3.5.2.2 (Excoffier & Lischer, 2010). In order to allow comparisons of diversity indices between the three species, only the 15 loci that consistently amplified in all of them were considered. Individuals were grouped into populations according to their geographic locations. Pairwise F-statistics were calculated as the number of different alleles between microsatellite genotypes, and statistical significance evaluated using a nonparametric test with 10,000 permutations. 95% confidence intervals for fixation indices were computed with 10,000 bootstraps using the R package hierfstat (Goudet, 2005). Isolation-by-distance (IBD) was tested with a Mantel test (n = 9999 iterations), conducted between the pairwise genetic distance matrix and the matrix of Euclidian geographical distances between individuals, as implemented in adegenet (Jombart, 2008). Pairwise genetic distances were estimated as the Euclidean Reynolds' distance (Re) (Reynolds, Weir, & Cockerham, 1983), while geographic distances were calculated with the R package geosphere (Hijmans, Williams, & Vennes, 2016), as orthodromic distances (in km) between geographical coordinates.

| Multilocus genotype datasets
A minimal number of 15 loci that were reliably amplifying in all three species were kept for the interspecific dataset. For the intraspecific datasets, we retained 15 to 19 loci depending on the species. Among the discarded loci, three (Paus03, Paus11, and Paus19) could not be reliably called in all three species because of amplification artifacts. Likewise, loci Paus01, Paus07, and Paus09 in P. auritus, Paur06 in P. austriacus and Paus04 and Paus07 in P. macrobullaris did not amplify well in those species. Finally, two loci in P. auritus (Paus13 and Paus16) and a single one in P. macrobullaris (Paus15) showed significant and consistent departure from HWE for several populations likely due to the presence of null alleles, and were thus also discarded for these species. Other signs of disequilibrium concerned only one or two populations per species and other loci were thus kept for intraspecific analyses (Supporting information Appendix S3). Inbreeding analyses carried with coancestry did not identify twins in the dataset (all estimated Δ 7 lower than 0.56). Moreover, the individuals exhibiting the highest kinship scores originated from different colonies or regions and were thus unlikely mother-pup pairs. Only two possible pairs of individuals having parent-offspring characteristics (Δ 7 = 0; Δ 8 > 0.96) were identified, but one concerned two lactating females of P. auritus sampled in the same breeding colony in Isère, and the other were two P. austriacus from distinct collection sites sampled at five years interval. Again they could likely represent related animals issued from successive years of reproduction, but since these cases were marginal, all the 349 individuals were kept in the dataset for subsequent analyses. The datasets overall contained less than 6% missing data, the mean number of alleles per locus ranged from 7.9 to 17.3 across species, with a mean heterozygosity ranging from 0.5 to 0.8 (Supporting information Appendix S4). The complete multilocus genotype dataset for the 349 individuals is provided in Supporting Data 2.

| Interspecific analyses
For the global dataset, the Evanno's criterion calculated from Bayesian clustering analyses based on 15 loci indicated K = 3 groups was the most likely structure (Supporting information Appendix S5A). Clustering of individuals into these three groups was fully concordant with the initial species assignation based on mitochondrial identifications, and corresponded, respectively, to 152 P. auritus, 79 P. austriacus and 118 P. macrobullaris. No sign of recent admixture was detected as all but one Q-values were greater than 0.96. This was notably the case for all P. auritus and P. macrobullaris individuals originating from breeding colonies located in the same building in Valais, or for other samples taken in close geographic proximity.
A single individual sampled in Piedmont (Italy) was assigned to the P. auritus cluster with a slightly lower Q-value (0.89), and to P. macrobullaris with a Q-value of 0.11. Increasing the K value did not reveal any additional meaningful cluster, but rather assigned isolated individuals to a mixture of clusters (Supporting information Appendix S6). Thus, a posteriori nuclear-based assignations to the three groups and initial mtDNA lineage identifications showed no cytonuclear discordance (Supporting information Appendix S6). The multivariate analysis indicated that the two first axes of the PCA carried about 10% of cumulated variance and represented the main structure of the dataset. The first axis perfectly discriminated all P. austriacus individuals from the others, while the second axis discriminated all P. auritus from P. macrobullaris, with no individual set in a genetically intermediate position (Figure 1a).
Since no sign of admixture or introgression was detected in P. austriacus with any of the methods used, power analyses of hybrid detection were only performed for the other two, most closely related species, namely P. auritus and P. macrobullaris. All estimates of S in pure parental forms were either very close to 0 or to 1, whereas H I was close to zero in all of these individuals (Figure 2). For 300 simulated F1 offspring individuals, S and H I ranged from 0.34 to 0.58 and from 0.68 to 1, respectively. There was no overlap between values observed for parental and simulated F1 individuals, providing

| Intraspecific analyses
For P. auritus, the highest ΔK was observed for K = 6 (Supporting

| Complete reproductive isolation among species
Despite several situations of strict sympatry among the three sampled species, and their conservative morphology, all conducted genetic analyses indicated that reproductive barriers between those species are not porous. Indeed, no recent hybrid was detected by the nuclear markers, including between individuals from distinct species sharing the same roost. A single individual from north-eastern Italy exhibited a slightly lower cluster assignment (Q = 0.89) due to the presence in this bat of otherwise uncommon alleles in P. auritus.
Our simulation analyses suggest that this intermediate assignment is unlikely due to past introgression of alleles from P. macrobullaris, but might rather reflect the existence of a slightly diverging local population of P. auritus. More samples from this area would be required to tackle this question more precisely. Additionally to the observed absence of contemporary nuclear gene flow among species, mitochondrial identifications of all individuals perfectly matched the species assignation inferred by microsatellite analyses, as each group was characterized by typical nuclear genotypes and carried speciesspecific mitochondrial lineage. This complete cytonuclear concordance indicates a lack of historical introgression between the three taxa. We thus confirm that the use of simple mitochondrial barcodes (typically COI) is perfectly suited to identify these long-eared bat species in Western Europe despite the existence of morphologically intermediate individuals . This concordance of genetic markers and lack of interspecific gene flow also validates previous genetic assignations that were exclusively based on mitochondrial markers (Ashrafi, Beck, Rutishauser, Arlettaz, & Bontadina, 2011;Mattei-Roesli, 2010). This clear-cut situation of strong barriers to gene flow among the three Plecotus contrasts with the situation of other pairs of cryptic species such as Myotis bats which hybridize and exchange both nuclear and mitochondrial genes (Berthier et al., 2006;Morales & Carstens, 2018;Morales et al., 2017).

| Porous geographical barriers for local populations
The Alps have been a strong barrier to gene flow for many organisms while they recolonized suitable habitats northwards after glaciations (Hewitt, 1999). Despite their potential for overtaking topographic obstacles, most species of bats tested so far recolonized Central and Northern Europe from the Southeast (Balkans) or Southwest (Iberia) but rarely from the southern Italian peninsula (the Apennines), probably due to the presence of the alpine range (Alberdi, Gilbert, et al., 2015b;Flanders et al., 2009;Petit, Excoffier, & Mayer, 1999;Razgour et al., 2013;Wright et al., 2018).
Only Barbastella barbastellus seems to have recolonized northern parts of its range from Italy (Rebelo et al., 2012), while fairly limited transalpine migration has been noted in few species such as M. myotis and M. blythii (Castella et al., 2001). For the lowland-adapted P. austriacus, Razgour et al. (2013) reported a strong differentiation of Italian versus Central and Western European populations, with directional gene flow from the Italian Peninsula toward Western Europe, following a species' range constriction and subsequent expansion due to climate change . Although a strong effect of the Alps as a blocking orographic barrier to gene flow is observed in this species (Razgour et al., 2014), limited gene F I G U R E 3 Bivariate plot of the relationship between genetic and geographical distances for the three species datasets used: (a) P. auritus, (b) P. austriacus, and (c) P. macrobullaris. The black line represents the linear curve best fitting the data. For readability purposes, a two-dimensional kernel density estimation heat map is superimposed to the scatterplot, with red for highest densities  Patthey, 2014;Ruedi et al., 1989). Clearly these areas provide mating opportunities for individuals issued from vast catchment areas (Rivers, Butlin, & Altringham, 2006), and thus potentially connects populations with high gene flow even across the Alps, as suggested by our analyses. However, the mating behavior is largely unknown in P. macrobullaris, and further knowledge is required to better understand factors promoting gene flow in this species.
Regarding the effect of open sea as a barrier to gene flow for longeared bats, our predictions arising from the patterns of island occupation in the Mediterranean is also largely corroborated by the results of nuclear genetic structure analyses. Except for P. auritus, which apparently does not occur on Corsica (Courtois et al., 2011), the effect of the Ligurian Sea is minimal for the potentially good colonizing species (P. austriacus), but more significant for the highland-adapted, apparently poor colonizer P. macrobullaris. Indeed, the few individuals of P. austriacus from mainland Italy appear to be genetically identical to animals from Corsica and differ from those sampled elsewhere in continental Europe (Supporting information Appendix S9 and Figure 1c).
This differentiation, however, also coincides with the Alpine range and might in part be due to an isolation-by-distance effect. As P. austriacus is very rare or absent from the intervening, northern regions of Italy immediately adjacent to the Alps, a proper sampling design to test these effects seems difficult to implement in the field. Clearly, populations sampled in south-eastern France, that is, in the area bordering the Ligurian Sea along the occidental limits of the Alps, could bring useful insights on gene flow and migration routes. Elsewhere, in the Italian Peninsula and Sardinia, no study involving nuclear variation have been conducted so far, but the limited information issued from mitochondrial lineages (Galimberti et al., 2012;Razgour et al., 2013) suggests that P. austriacus from this region are not particularly differentiated.
Regarding the effect of open sea barriers, an opposite situation prevails for P. macrobullaris, as all individuals from Corsica form a distinct cluster relative to those from the continent (Supporting information Appendix S10 and Figure 1d). The only exception is an animal from Porto-Vecchio in south-eastern Corsica, which exhibits a multilocus genotype identical to that of continental samples. This exceptional similarity is not due to an analytical error, as its genotype was extracted and characterized twice independently without visible inconsistency. This animal might therefore represent a recent migrant from the continent, indicating that the Ligurian Sea is not an absolute barrier for P. macrobullaris, although such occasional gene flow is not sufficient to counter the local differentiation of Corsican populations.
Our analyses of nuclear DNA show no genetic breaks between populations bearing these lineages, suggesting that individuals are interbreeding regardless of their mtDNA lineage. These mtDNA variants are better considered as deep conspecific lineages (Padial, Miralles, Riva, & Vences, 2010), which mirror the situation found in other bats from this region (Andriollo, Naciri, & Ruedi, 2015). Likewise, without evidence from biparentally inherited markers, it seems premature to consider the Iberian lineage associated to the taxon begognae (Ibáñez, et al., 2006) as specifically distinct on the sole basis of mitochondrial data (Mayer, Dietz, & Kiefer, 2007). The apparent contrast between strong population structure inferred from mitochondrial markers versus weak differentiation inferred with nuclear markers is not uncommon in bats and may be a consequence of a stronger female philopatry compared to male migration/dispersal (Burland & Worthington-Wilmer, 2001). Long-eared bats fit this pattern, as shown in P. auritus, where the strong philopatry of females (which transmit the mitochondrial lineages to the next generations) evidenced in the breeding colonies is counterbalanced by extensive male-mediated gene flow (Burland, Barratt, Beaumont, & Racey, 1999;Burland, Barratt, Nichols, & Racey, 2001). This increased malemediated gene flow is particularly evident at the autumn and spring swarming sites, where females from distant colonies may mate with unrelated males (Furmankiewicz & Altringham, 2007).
The other two species of long-eared bats, P. austriacus and P. macrobullaris, do not exhibit any major phylogeographic pattern at mtDNA markers in Central Europe, although major distinct lineages exist in Iberia for P. austriacus (Razgour et al., 2013)  in order to preserve distinct potentials for evolution, each of these gene pools should be considered as evolutionarily significant units (Moritz, 1994), although again no mitochondrial differentiation was observed at this geographic scale (Kiefer, 2007).

| CO N C LU S I O N S
Our genetic analyses, combining mitochondrial and nuclear markers, showed no evidence of interspecific admixture between the three long-eared bat species in the Alpine and adjacent areas, where their distributions overlap extensively. These cryptic species therefore mate assortatively based on other clues than external morphology.
Furthermore, despite striking phenotypic similarities, species-specific responses to topographic barriers could be evidenced, which correlated to distinct local altitudinal preferences of the species.
The Alps, which act as an orographic barrier to gene flow in several European bats, were surprisingly porous for two of the long-eared bat species studied (P. auritus and P. macrobullaris), but more effective in the third (P. austriacus). Similarly, the 100-km-wide channel of the Ligurian Sea also appeared to be a variably porous barrier for the two species sampled in Corsica (P. austriacus and P. macrobullaris). Overall, at the geographic scale envisioned here, these allegedly poor dispersers seem to cope well with these topographic barriers, but to different degrees. These differences led to idiosyncratic phylogeographic histories in the three species, resulting in distinct contemporary genetic structuring of populations. These intraspecific structures inferred from nuclear data largely contradict those suggested by previous analyses based on mitochondrial markers. Although three major mtDNA lineages were documented in P. auritus, all belong to a single nuclear genetic pool and should therefore be managed as a single unit of conservation. Conversely, a single mitochondrial lineage characterized each of P. austriacus and P. macrobullaris, yet microsatellites data segregated two welldefined ESUs within both these species. This highlights again the importance of defining priority conservation units on the basis of multiple markers. Simple barcode approaches are certainly useful for inferring phylogeographic patterns, but not particularly suited for estimating levels of gene flow among lineages.

ACK N OWLED G M ENTS
We warmly thank Orly Razgour who provided essential information on microsatellites genotyping, and Andrea Galimberti, David Progin and Roberto Toffoli who shared critical tissue samples. We also gratefully acknowledge the multiple colleagues for guiding or helping us during fieldwork, including (

E T H I C A L S TAT E M E N T
Most tissue samples were issued from museum collections or from biopsy samples (3 mm in diameter) achieved from monitoring programs. Additional biopsy samples specifically collected for this study were performed under appropriate licenses and ethical approval ("autorizzazione cantonale TI-05-2016" in Ticino, Switzerland; "Programme Régional de Conservation des Chiroptères en Corse" in Corsica, France; "arrêté préfectoral n o 38-2016-06-28-005" in Isère, France). Capture methods included mist nets, harp traps, and hand-net (Kunz & Kurta, 1988;Kunz, Hodgkison, & Weise, 2009).
Animals were retained for a maximum of 20 min in cotton bags, before being released on the spot of capture.

CO N FLI C T O F I NTE R E S T
None declared.