# The effect of phylogeographic history on species boundaries: a comparative framework in Hyla tree frogs

Christophe Dufresnes, Matthieu Berroneau, Sylvain Dubey, Spartak N Litvinchuk, Nicolas Perrin

► To cite this version:

Christophe Dufresnes, Matthieu Berroneau, Sylvain Dubey, Spartak N Litvinchuk, Nicolas Perrin. The effect of phylogeographic history on species boundaries: a comparative framework in Hyla tree frogs. Scientific Reports, 2020, 10 (1), pp.5502. 10.1038/s41598-020-62382-4. hal-04886987

HAL Id: hal-04886987

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

Submitted on 14 Jan 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>
OPEN

# The effect of phylogeographic history on species boundaries: a comparative framework in *Hyla* tree frogs

Christophe Dufresnes¹,²*, Matthieu Berroneau³, Sylvain Dubey²,⁴,⁵, Spartak N. Litvinchuk⁶,⁷ & Nicolas Perrin²

Because it is indicative of reproductive isolation, the amount of genetic introgression across secondary contact zones is increasingly considered in species delimitation. However, patterns of admixture at range margins can be skewed by the regional dynamics of hybrid zones. In this context, we posit an important role for phylogeographic history: hybrid zones located within glacial refugia (putatively formed during the Late-Pleistocene) should be better defined than those located in post-glacial or introduced ranges (putatively formed during the Holocene and the Anthropocene). We test this hypothesis in a speciation continuum of tree frogs from the Western Palearctic (*Hyla*), featuring ten identified contacts between species spanning Plio-Pleistocene to Miocene divergences. We review the rich phylogeographic literature of this group and examine the overlooked transition between *H. arborea* and *H. molleri* in Western France using a multilocus dataset. Our comparative analysis supports a trend that contacts zones resulting from post-glacial expansions and human translocations feature more extensive introgression than those established within refugial areas. Integrating the biogeographic history of incipient species, i.e. their age since first contact together with their genetic divergence, thus appears timely to draw sound evolutionary and taxonomic inferences from patterns of introgression across hybrid zones.

How and when the continuous process of speciation leads to discrete entities is a central question in evolutionary biology and systematics¹. Reproductive isolation evolves as diverging lineages accumulate genetic changes, ultimately reducing the probability of mating (pre-zygotic barriers) and the fitness of interspecific hybrids (post-zygotic barriers)². Hybridizing taxa thus lay at the boundaries of the biological species concept³ and offer opportunities to understand how much divergence is needed to maintain incipient species apart, despite occasional gene flow⁴.

Since it reflects the balance between dispersal and selection against hybrids, the amount of admixture at phylogeographic transitions serves as an ad hoc measure of reproductive isolation under natural conditions. Although scarce, comparative hybrid zone analyses support a gradual buildup of reproductive isolation with the time diverging lineages have spent in allopatry, where a remarkably short temporal window separates widely admixing and genetically isolated taxa⁵⁻⁷. Characterizing this “grey zone of speciation” is attractive for species delimitation⁷. Hence, assessing the amount of admixture at phylogeographic transitions is increasingly viewed as a key step in integrative taxonomy, notably for controversial species⁸.

However, patterns of introgression at hybrid zones can also be influenced by an array of extrinsic factors that confound with reproductive isolation, e.g. local barriers to gene flow due to landscape features, demographic events of extinction and recolonization, or changes in ecological conditions shifting species transitions⁶,⁹⁻¹². We posit that the time since contact should also play a major role, in the light of the biogeographic history of species. In the Northern hemisphere, the Quaternary climatic oscillations frequently remodeled lineage distributions

¹LASER, College of Biology and the Environment, Nanjing Forestry University, Nanjing, China. ²Department of Ecology & Evolution, University of Lausanne, Lausanne, Switzerland. ³Cistude Nature, Le Haillan, France. ⁴Agrosustain SA, Nyon, Switzerland. ⁵Hintermann & Weber SA, Montreux, Switzerland. ⁶Institute of Cytology, Russian Academy of Sciences, Saint Petersburg, Russia. ⁷Dagestan State University, Makhachkala, Russia. *email: Christophe.Dufresnes@hotmail.fr
and set up the contact zones observed nowadays¹³. Species transitions located in glacial refugia may have existed since or even prior to the last glaciation, and may thus be as old as several hundred thousand of years. In contrast, contact zones formed by post-glacial recolonization post-date the Last Glacial Maximum (LGM, ~21,000 years ago), and are no older than a few thousands of years. Post-glacial expansions should lead to wide, mosaic parapasic distributions as species colonized new, empty ranges, across the climatically instable northern regions, giving rise to multiple local contacts that may not have reached their ecological or genetic equilibrium yet (e.g.¹²). As a result, these transitions would exhibit admixture over large distances, even if the interacting species are partially reproductively isolated. Reciprocally, contact zones that have persisted in glacial refugia had more time to purge hybrid genotypes and should feature more stable transitions. These processes could also affect opportunities for reinforcement, i.e. the evolution of assortative mating as a response to selection against unfit or maladapted hybrids¹⁴,¹⁵. Species could start establishing pre-mating barriers earlier in refugia than in post-glacial ranges, where they are thus expected to experience less gene flow.

Given these contrasting outcomes, did the biogeographic history of secondary contacts significantly affect their amount of admixture, in turn blurring species delimitation assessments? This question is fundamental for both speciation and systematic biology, but because the answer requires replicate transitions between the same species pairs, or at least between pairs of similar divergence, it has never been formally tested. Several unrelated studies yet support remarkable latitudinal differences in gene flow. The rate of hybridization between the collared and pied flycatchers (*Ficedula*) is higher in the Baltic Isles than in Central Europe¹⁶, where both species have been in sympatry for a longer time and evolved character displacement on male plumage color via reinforcement¹⁷. Time since contact was also invoked to explain differences in hybrid zone structures among *Quercus* oaks¹⁸. In house spiders (*Tegenaria*), two sister species co-exist in England and widely hybridize in northern, but not in southern parts of the country¹⁹. These scattered observations match our expectations of wider hybrid zones in recently colonized ranges, for the same species pair.

The Western-Palearctic radiation of tree frogs (*Hyla*) provides a unique framework to compare phylogeographic transitions in respect to their late-Quaternary history. The phylogeny and phylogeography of this group have been extensively characterized over the last decade, for all species individually, but also at range margins. W-Palearctic tree frogs diversified from the Miocene (>10 Ma) to the Plio-Pleistocene (3-4 Ma)²⁰ and now include ten species forming replicate transitions across refugial, post-glacial, or even introduced ranges. Here we revisit their patterns of admixture in the light of the history of contacts. We first detail the overlooked hybrid zone between *H. arborea* and *H. molleri* preliminary located in W-France²¹⁻²³, using a population genetic approach. This pair is of high comparative interest since it diverged almost simultaneously (5-6 Ma) with several other species pairs (*H. arborea/orientalis*, *H. arborea/perrini*, *H. savignyi/felixarabica*). Second, we assess how admixture across phylogeographic transitions varies depending on species divergence and relative time since first contact, by combining independently published data from ten contact zones.

## Methods

**Population genetics of the Hyla arborea/molleri hybrid zone.** A total of 251 adult tree frogs (*n* = 231 *H. arborea/molleri*, as well as 20 syntopic *H. meridionalis*) were captured from 25 localities (loc. 4–29) throughout Western France (File S1). DNA was collected with non-invasive buccal swabs, stored dry at –20 °C, and extracted with the BioSprint robotic workstation (Qiagen). The experiments conducted in this research (DNA sampling) were performed in accordance with the regional authorities (Prefectures of Ille-et-Vilaine, Charente and Gironde) under collecting permits, and relevant guidelines were followed. No individuals were harmed or killed, and all were released at their place of capture.

Twenty-six microsatellite loci cross-amplifying in *H. arborea* and *H. molleri* were genotyped in all samples: *Ha*-T3, *Ha*-T52, *Ha*-M2, WHA5–22, *Ha*-T45, *Ha*-T11, *Ha*-T49, *Ha*-T51, *Ha*-T64, *Ha*-T41, *Ha*-T69, *Ha*-T54, *Ha*-T55, *Ha*-T50, *Ha*-T58, *Ha*-T53, *Ha*-T56, *Ha*-T66, *Ha*-T60, *Ha*-T61, *Ha*-T68, *Ha*-T63, *Ha*-T67, WHA1–25, WHA1–67, *Ha*-A11²⁴,²⁵ (see methods therein). Eight loci did not amplify in *H. meridionalis*, in line with the older divergence time. In addition, we re-analyzed microsatellite reads from reference populations of *H. arborea* (from Switzerland, the Netherlands, and France, loc. 1–3, *n* = 80), *H. molleri* (from Spain, loc. 30, *n* = 16) and *H. meridionalis* (from southern France, loc. 31, *n* = 16) from our previous work²⁴⁻²⁶. The new and published genotypes were generated by the same genetic analyzer (an ABI3130) and (re)scored together in GenMapper 4.0 with a single set of bins and panels, to ensure the correspondence of alleles.

The total microsatellite dataset (Western France + reference populations) included 363 individuals. We first performed Bayesian clustering with STRUCTURE²⁷ by running 10 chains of 100'000 iterations, after a burnin of 10'000, for each value of K from 1 to 11. The most likely number of clusters (K) was estimated with STRUCTURE HARVESTER²⁸, concatenated over replicates with CLUMPP²⁹, and graphically displayed with Distruct³⁰. Second, we performed a Principal Component Analysis (PCA) of individual genotypes with the R packages *ade4* and *adegenet*³¹. Third, we computed pairwise genetic distance (*F*<sub>st</sub>) for populations with *n* ≥ 5 using the R package *hierfstat*³².

Mitochondrial DNA from 230 hybrid zone samples was profiled by enzyme restriction with *MseI*, which differentially cuts the three mitotypes (methods in³³). In addition, we included the published mtDNA barcoding data from two additional populations located close to our loc. 15–16 (*n* = 9 *H. arborea/molleri* and 3 *H. meridionalis*;²² see File S1).

Finally, we fitted sigmoid clines to mitochondrial frequency and nuclear admixture coefficients (STRUCTURE's Q) along a north-south transect running along the Atlantic coast (loc. 4–12, 15–16, 21–24, 26–29), using the R package *hzar*³⁴. We performed model selection between two- (width and center) to eight-parameter clines (width, center, parental frequencies, position and steepness of exponential tails), and kept the model with the best AIC score.
**Phylogeography of W-Palearctic Hyla.** In order to provide an up-to-date and comprehensive overview of species distributions, range limits, and the extent of the ten contacts documented so far in W-Palearctic *Hyla*, we combined our new *H. arborea/molleri* population genetic data with those from 24 previously-published studies on the phylogeography and hybrid zones of this group, totaling records from 857 different localities. These included general phylogenetic/phylogeographic accounts<sup>21,23,35</sup>, intra-specific genetic work on *H. meridionalis* and *H. carthaginiensis*<sup>36,37</sup>, *H. savignyi* and *H. felixarabica*<sup>38,39</sup>, *H. sarda*<sup>40,41</sup>, *H. intermedia* and *H. perrini*<sup>42,43</sup>, *H. arborea*<sup>26,44</sup>, *H. molleri*<sup>45,46</sup>, and *H. orientalis*<sup>38,39,47</sup>, as well as targeted surveys of the transitions for *H. meridionalis/arborea-molleri*<sup>22</sup>, *H. savignyi/felixarabica*<sup>48</sup>, *H. arborea/perrini*<sup>49,50</sup>, *H. intermedia/perrini*<sup>20</sup> and *H. arborea/orientalis*<sup>33,51</sup>. Based on the distribution of genetic diversity and ecological niche modelling, these studies also inferred the location of glacial refugia and routes of post-glacial recolonization for each species, providing insights on the relative age of contact zones that we summarized in File S2.

For divergence times, we relied on the fully-resolved phylogeny of<sup>20</sup>, built from genome-wide nuclear sequences (~42 kb), and a molecular clock calibrated with the tree root at ~20 Mya<sup>52,53</sup>, and the Messinian divergence of the Tyrrhenian *H. sarda*<sup>21,35,54</sup>.

*Hyla* species transitions received heterogeneous assessments in terms of methodology (numbers and nature of loci, type of analyses) and sampling schemes. Therefore, rather than a quantitative measure of hybrid zone width (which was available for only three out of the ten contacts), here we distinguished three discrete categories reflecting the type of transition. These categories were: (1) sympatry without gene flow; (2) narrow transitions (<50 km) with admixture (if detected) restricted to few localities at the contact; (3) wide transitions (>50 km) involving widespread admixture far away from the contact.

**Species distribution modelling.** To complement previous knowledge on the putative age of contact zones, we built species distribution models (SDM) to predict the ecological niches and the LGM ranges for all ten *Hyla* taxa. Such analyses have been independently conducted in previous work on some taxa<sup>20,46-48</sup>, but here we aimed to provide comparable outcomes across all of them, as obtained from a single methodology. To this end, analyses were performed with MaxEnt 3.4.1<sup>55</sup>, implementing occurrence filtering, tests for distinct candidate sets of environmental variables, and multiple combinations of model parameters (feature classes, regularization multipliers, and sets of variables), using accessible areas for model calibration and multiple statistical criteria (partial ROC, omission rates, and AICc) for model selection.

A total of 6,363 localities were used, comprising our own records and previously published data (File S3), filtered to avoid spatial autocorrelation and duplication by NicheToolBox<sup>56</sup>. We selected localities of each species to obtain the maximum number of localities that were at least 10 km (0.093°) apart, which resulted in 3,685 presence-only locations for analyses. Altitude and 19 bioclimatic layers (30 arc seconds and 2.5 arc minutes spatial resolutions) representative of the climatic data over ~1950–2000 were extracted from the WorldClim 1.4 database (http://www.worldclim.org). An additional three layers were considered: the aridity index (Global Aridity and Potential Evapo-Transpiration; http://www.cgiar-csi.org/data/global-aridity-and-pet-database), spatial homogeneity of global habitat (EarthEnv; http://www.earthenv.org/texture.html), and the global percent of tree coverage (https://github.com/globalmaps/gm_ve_v1). To consider topography in the model, four landscape layers were calculated with QGIS: aspect, exposition, slope, and terrain roughness index. For the LGM predictions, we used the same 19 bioclimatic layers (2.5 arc minutes spatial resolution) and applied two general atmospheric circulation models: the Community Climate System Model (CCSM; http://www2.cesm.ucar.edu/) and the Model for Interdisciplinary Research on Climate (MIROC<sup>57</sup>). Analyses were conducted with species-specific masks covering to the area of occurrence and adjacent regions.

To eliminate predictor collinearity before generating the model, we calculated Pearson's correlation coefficients for all pairs of bioclimatic variables using ENMTools<sup>58</sup>. For correlated pairs ($|r| > 0.75$), we excluded the variable that appeared less biologically important for *Hyla* tree frogs. Analyses were performed separately for each species, and we ran MaxEnt models with 10 replicate bootstrap resamplings. Model calibration consisted in evaluating models created with distinct regularization multipliers (0.5 to 6 at intervals of 0.5), feature classes (resulted from all combinations of linear, quadratic, product, threshold, and hinge response types), and from two different sets of layers. The first set included altitude, the aridity index, spatial homogeneity, the percent of tree coverage, the four landscape variables, as well as the uncorrelated bioclimatic layers; the second was restricted to the most important layers, i.e. those with a relative contribution above 10% in preliminary analyses. A total of 696 candidate models were thus built for each species. The best parameter settings were selected considering statistical significance (partial ROC), predictive power (omission rates $E = 5\%$), and complexity level (AICc) obtained using the R package *kuenm*<sup>59</sup>. Finally, territories covered by glaciers and seas at LGM times were cut out from the corresponding distributions. These are illustrated in File S4, based on available data<sup>60,61</sup>.

## Results

**Hyla arborea/molleri hybrid zone.** Our data supported widespread nuclear admixture between *H. arborea* and *H. molleri*, spanning from the Garonne basin (loc. 21) to the Massif Central (east of loc. 14) and Normandy (loc. 4) (Fig. 1). The mitochondrial picture was broadly concordant, with mitotypes shared in several localities north of the Garonne (loc. 15–20). Clines were wide (width $w = 98$ km for nuclear loci; $w = 43$ km for mtDNA), remarkably coincident (center $c = 407$ km for nuclear loci; $c = 406$ km for mtDNA), and both bear a long introgression tail on the *H. arborea* side (the selected cline model under the AIC criterion, Fig. 2). There was no trace of mitochondrial or nuclear introgression with the sympatric and syntopic *H. meridionalis* (Fig. 1). The PCA (File S5) and pairwise $F_{ST}$ (File S6) disentangled the three gene pools and confirmed that individuals collected in Western France represent the entire *H. arborea–molleri* spectrum.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Population genetics of the *H. arborea/molleri/meridionalis* contact zone in Western France. Proportion of each mtDNA haplotype (left) and the average admixture coefficient of each population based on 26 microsatellites (right); barplots show individual admixture coefficient (STRUCTURE) for the study area and reference samples; pie charts are proportional to sample sizes. The maps were created in QGIS 3.4 (https://qgis.org).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Cline analysis of mitochondrial frequency (left) and nuclear genome average (right) across the *H. arborea/molleri* hybrid zone along a north-south geographic transect (shown on the inset map); for both datasets, the selected cline model involved asymmetric introgression tails, which were larger on the *H. arborea* side; encircled crosses illustrate sample size for each population; shaded areas represent the 95% confidence interval of fitted clines.</figcaption></figure>

**Species transitions.** The distribution of W-Palearctic *Hyla* taxa is shown in Fig. 3, based on genetic bar-coding from 25 studies (including the present paper). A summary of species interactions and the corresponding literature is provided in Table 1 and File S2, from which we categorized the following. (1) One case of sympatry without gene flow despite occasional hybridization: *H. meridionalis/arborea-molleri*; (2) Five steep parapatric transitions, with introgression (if any) restricted at range margins, and narrow cline width (<50 km): *H. savignyi/felixarabica* in N-Israel; *H. arborea/perrini* in NE-Italy; *H. meridionalis/carthaginiensis* in NE-Algeria; *H. arborea/orientalis* in the Balkans; *H. savignyi/orientalis* in S-Turkey, where evidence for hybridization/admixture is
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Distribution of W-Palearctic tree frogs based on genetically barcoded populations, combining 25 studies. At transitions, populations were considered as a mix (represented in half circles) if they featured both mtDNA haplotypes and/or nuclear genotypes consistent with admixture (STRUCTURE's Q < 0.9 in multilocus analyses; co-segregation of species-diagnostic alleles). The map was created in QGIS 3.4 (https://qgis.org). Photo: *Hyla meridionalis* from southeastern France, by CD.</figcaption></figure>

<table><caption>**Table 1.** Summary of ten secondary contact zones among W-Palearctic tree frogs, in respect to the divergence time between species pairs (Div. in Ma; taken from<sup>20</sup>; nodes are labelled in Fig. 6), the putative epoch of first contact based on previous studies and past species distribution modelling (see File S2), and the type of transition (see File S2, Fig. 3 and the listed references), with the corresponding category in our meta-analysis (a sympatry without gene flow, b steep transition, c wide transition/hybrid swarm); w: hybrid zone width, as measured from cline analyses of transects.</caption><thead><tr><th>Node</th><th>Region</th><th>Species pair</th><th>Div.</th><th>Contact</th><th>Transition</th><th>Reference</th></tr></thead><tbody><tr><td>A</td><td>W-France</td><td><i>H. meridionalis/arborea-molleri</i></td><td>20</td><td>Antiquity</td><td>sympatry without gene flow<sup>a</sup></td><td>22</td></tr><tr><td>B</td><td>S-Anatolia</td><td><i>H. orientalis/savignyi</i></td><td>8.8</td><td>≥LGM</td><td>steep, hybridization/admixture not detected<sup>b</sup></td><td>38,39</td></tr><tr><td>C</td><td>NE-Italy/Slovenia</td><td><i>H. arborea/perrini</i></td><td>6.6</td><td>≥LGM</td><td>steep, faint traces of admixture<sup>b</sup></td><td>49</td></tr><tr><td>C</td><td>W-Switzerland</td><td><i>H. arborea/perrini</i></td><td>6.6</td><td>1950s</td><td>hybrid swarm<sup>c</sup></td><td>50</td></tr><tr><td>D</td><td>Balkans</td><td><i>H. arborea/orientalis</i></td><td>5.1</td><td>≥LGM</td><td>steep (w = 30–32 km)<sup>b</sup></td><td>33</td></tr><tr><td>D</td><td>Poland</td><td><i>H. arborea/orientalis</i></td><td>5.1</td><td>post-glacial</td><td>wide, admixture across a large area (>200 km)<sup>c</sup></td><td>51</td></tr><tr><td>D</td><td>W-France</td><td><i>H. arborea/molleri</i></td><td>5.1</td><td>post-glacial</td><td>wide (w = 98 km)<sup>c</sup></td><td>this study</td></tr><tr><td>E</td><td>Central Italy</td><td><i>H. perrini/intermedia</i></td><td>3.5</td><td>post-glacial</td><td>wide (w = 96 km)<sup>c</sup></td><td>20</td></tr><tr><td>F</td><td>Levant</td><td><i>H. savignyi/felixarabica</i></td><td>6.4</td><td>≥LGM</td><td>steep, geographically restricted admixture<sup>b</sup></td><td>38,48</td></tr><tr><td>G</td><td>E-Algeria</td><td><i>H. meridionalis/carthaginiensis</i></td><td>4.1</td><td>≥LGM</td><td>steep, admixture detected at a single locality<sup>b</sup></td><td>37</td></tr></tbody></table>

presently lacking. (3) Three wide continuous transitions with large nuclear clines (>50 km) and admixture over more than a hundred kilometers: *H. arborea/orientalis* in Poland; *H. arborea/molleri* in SW-France; *H. intermedia/perrini* in Central Italy; the *H. arborea/perrini* hybrid swarm from W-Switzerland was also assigned to this category. Considering the ten contacts altogether, the relationship between the amount of admixture and the divergence time was significant (ordinal logistic regression: LR χ² = 7.72, n = 10, df = 1, P = 0.005).

**LGM distributions and relative ages of the species transitions.** All species distribution models (SDM) received very good fits, indicative of robust predictions (File S7). Projected distributions for each species are shown in Figs. 4 and 5. In Europe, LGM conditions appeared unsuitable across the northern parapatric ranges of *H. arborea* and *H. orientalis*, but were more favorable in the Balkans (Fig. 4). For *H. molleri*, the models identified hospitable glacial ranges along the coasts (Fig. 4). In Italy, the probabilities of occurrence remained high for
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Projected distributions under present and LGM conditions (MIROC and CCSM scenarios) for European *Hyla* species. Color warmth (from white to brown) reflects the probability of occurrence. Dark lines show current species ranges. The maps were created in QGIS 3.4 (https://qgis.org).</figcaption></figure>

*H. perrini*, but not for *H. intermedia*, which appeared scattered throughout the Peninsula (Fig. 4). In N-Africa, the Maghreb coast offered favorable conditions for both *H. meridionalis* and *H. carthaginiensis* (Fig. 5). Finally, the present and LGM predictions were globally similar for the insular *H. sarda*, and the Middle-Eastern *H. savignyi* and *H. felixarabica* (Fig. 5).

Integrating the projected species distributions with the well-documented phylogeography of W-Palearctic *Hyla* (File S2), the wide transitions/hybrid swarms were all presumably formed after the LGM (Holocene or Anthropocene contact): *H. arborea/molleri*, *H. arborea/orientalis* in Poland, *H. intermedia/perrini*, *H.*
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Projected distributions under present and LGM conditions (MIROC and CCSM scenarios) for Tyrrhenian, Middle-Eastern and N-African *Hyla* species. Color warmth (from white to brown) reflects the probability of occurrence. Dark lines show current species ranges. The maps were created in QGIS 3.4 (https://qgis.org).</figcaption></figure>

perrini/arborea in W-Switzerland (Fig. 6). In contrast, the steep transitions all encompass putative glacial refugia (Pleistocene contact): *H. savignyi/orientalis*: *H. arborea/orientalis* in the Balkans; *H. arborea/perrini*; *H. savignyi/felixarabica*; *H. meridionalis/carthaginiensis* (Fig. 6). Finally, the sympatric (but not admixing), early-diverged pair *H. meridionalis/arborea-molleri* result from a human introduction (File S2). Accordingly, the relative age of contact (Holocene/Anthropocene vs Pleistocene) was a significant predictor of these three categories of transition (wide, steep, or sympatry without gene flow), based on an ordinal logistic regression (LR χ² = 3.85, n = 10, df =
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Phylogeny and distribution of Western-Palearctic tree frogs, illustrating post-glacial expansions and the extent of introgression in secondary contact zones in respect to their amount of divergence (capital letters: phylogenetic nodes) and age since contact (subscripts: "young" as post-LGM (y) or "old" as pre-LGM (o)). Species distributions are adapted from<sup>77</sup>. See Table 1 for details on the patterns of admixture and relevant references. The map was created in QGIS 3.4 (https://qgis.org). Photo: *H. arborea* from southeastern France, by CD.</figcaption></figure>

1, *P* = 0.049), and as illustrated in Fig. 7. The effect was also significant when comparing only wide vs steep hybrid zones with a χ² test of independence (χ² = 5.41, *n* = 9, df = 1, *P* = 0.020) – two categories where species pairs do not significantly differ in divergence times (Mann-Whitney-Wilcoxon test, *n* = 9, *W* = 13.5, *P* = 0.45).

**Discussion**

**A wide tree frog transition across W-France.** The widespread introgression between *H. arborea* and *H. molleri* documented across W-France confirms previous studies that reported hybrid specimens in several scattered localities<sup>21-23</sup>. Their genomes are thus mostly compatible, so selection against hybrids is probably weak (see below). The phylogeographic history of these two species suggests a somewhat classic scenario of post-glacial contact. According to our predicted distributions (Fig. 4, see also<sup>46</sup>), the Iberian *H. molleri* could have expanded from some Atlantic refugia into western France (where private mtDNA haplotypes are found<sup>21,46</sup>), while *H. arborea* arrived from the more distant Balkan refugium<sup>26</sup>. Subsequently, these two species introgressed massively along the Atlantic coast, but the dry Garonne valley may nowadays act as an effective bioclimatic barrier preventing northwards dispersal of *H. molleri* (as in other biota<sup>62</sup>). This hypothesis is supported by the asymmetric transition (Fig. 2): *H. arborea* allele frequencies decrease continuously with latitude, up to a sharp drop between the northern (loc. 4–20) and southern sides of the valley (loc. 21–29). Hence, the hybrid zone is probably not stable, and *H. arborea* might eventually fully backcross the hybrid populations over time. This particular scheme of hybrid zone movement is somewhat reminiscent of the *Bufotes viridis/balearicus* green toad transition in Italy, where hybrids became isolated north of the Po River, and cannot escape inflocks of *B. viridis* expanding southward<sup>63</sup>. Here, the *H. arborea/molleri* species pair adds a valuable point to our comparative framework of admixture vs history of secondary contacts.

**Tree frogs admix freely in young hybrid zones, but not in old ones.** The dense phylogeographic framework, including ten contacts from the fully-resolved radiation of W-Palearctic tree frogs, supports that the width of species transitions drastically differs depending on when these were first established (Figs. 6 and 7, Table 1). Lower-Miocene lineages do not introgress (node A), as expected given the deep genetic (and phenotypic) divergence. However, Upper Miocene (nodes B–D and F) and Pliocene (nodes E–G) taxa admix significantly more across post-glacial and introduced ranges than in refugial areas (Figs. 6 and 7). Note that the divergence times reported here (from<sup>20</sup>, see Methods) should be taken with caution, since the genus *Hyla* might be younger than previously assumed<sup>64</sup>.

Because they feature hybrid zones involving the same pair, or between species that coalesced at the same time, two nodes are of particular interest for direct comparisons. First, we highlighted striking differences among the *H. arborea/orientalis* contacts (node D). In their Balkan refugium, mtDNA and nuclear transitions do not exceed 30 km<sup>33</sup>, while they extend over hundreds of kilometers in the Polish post-glacial ranges<sup>51</sup>. The latter mirrors the wide post-glacial hybrid zone documented herein for *H. arborea/molleri*, which share the same amount of divergence (node D). Second, the slightly older *H. arborea* and *H. perrini* (node C) do not admix along their natural margins in NE-Italy/Slovenia<sup>49</sup>, where subrefugia of both species are putatively located<sup>20,26</sup> (see File S2). Nevertheless, this pair widely hybridized in W-Switzerland<sup>50</sup> following introductions of *H. perrini* in the 1950s<sup>65</sup>. All other tree frogs meeting in ecologically stable regions featured steep mitochondrial and nuclear transitions (File S2), i.e. *H. savignyi/orientalis* in southern Anatolia (note B<sup>38,39</sup>), *H. savignyi/felixarabica* in the Levant (node
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Species transitions in respect to divergence time in W-Palearctic tree frogs. Color codes refer to Figs. 3 and 6, and replicate contacts between the same species are specifically labelled. Letters indicate the phylogenetic nodes shown in Fig. 6. Symbols differentiate the contacts according to when they were putatively formed.</figcaption></figure>

F³⁸,⁴⁸, and *H. meridionalis/carthaginiensis* in the E-Maghreb (node G³⁷). In Italy, the *H. intermedia/perrini* contact is presumably post-glacial (see details in File S2), but there the wide transition may also result from the weak reproductive barriers between these recently-diverged taxa (~3.5 Ma²⁰).

Different opportunities for dispersal would hardly explain these contrasting species transitions. European tree frogs are cosmopolitan and some of the most mobile amphibians (up to 12 km movements⁶⁶,⁶⁷). While large rivers and their valleys were identified as contemporary geographic barriers, this involves both refugial and post-glacial transitions, e.g. the Vistula between *H. arborea* and *H. orientalis* in Poland⁵¹; the Garonne between *H. arborea* and *H. molleri* (this study); the Isonzo between *H. arborea* and *H. perrini*⁴⁹; the Jordan Valley between *H. savignyi* and *H. felixarabica*³⁸,⁴⁸. Reciprocally, no geographic barriers were identified for some narrow transitions, e.g. *H. meridionalis/carthaginiensis*³⁷, *H. arborea/orientalis* in the Balkans³³. Ecological barriers also appear unlikely, given that the projected distributions under current conditions remarkably overlapped between parapatric species, notably at their transitions (Figs. 4 and 5).

Several limitations prevent a more quantitative assessment of the relationship between the time since contact and the patterns of introgression across the radiation, which would strengthen our results. First, the heterogeneity of the hybrid zone datasets published on *Hyla* (uneven sampling, diverse types of molecular markers) precludes comparable and spatially-explicit genetic reconstructions of the demographic history of hybridizing populations (e.g.⁶⁸), which could have provided estimates of the ages since first contact. Genomic phylogeographic data (e.g. RAD-seq) and transect sampling for all hybrid systems would have also allowed to compare admixture patterns with more accuracy. Second, our categorization of the relative ages of contact, taken from phylogeographic data and SDM predictions, assumes that species niches did not change since the LGM, and that old hybrid zones remained static through time. Accordingly, it is possible that some contact zones located within southern glacial refugia were established post-LGM following regional expansions from separate micro-refugia, even if they seem connected by suitable LGM conditions in our analyses. Yet, these southern contacts should still be older than the Holocene northern contacts located thousands of kilometers away from their source populations. Third, and while we assume that local geographic barriers played a minor effect on the patterns of introgression (see above), this remains to be thoroughly addressed by habitat analyses, notably in the context of landscape changes during the Anthropocene. In parallel, the relative contribution of pre- vs post-zygotic isolation in limiting admixture across the narrow hybrid zones is unclear; characterizing pre-mating barriers (breeding calls, pheromones) between supposedly cryptic incipient species remains a major gap in the amphibian literature. Accounting for all the factors potentially affecting the age and admixture patterns at secondary contact zones is already colossal for a single transition, and it becomes logistically impossible in a comparative framework involving multiple contacts.

While acknowledging these limitations, the qualitative data in hand thus tend to support our working hypothesis that the contrasted hybrid zones in *Hyla* partly result from the phylogeographic dynamics that led to their establishment. Compared to refugial areas, colonizers at the front of post-glacial expansions should hybridize more with related species⁶⁹, due to limited mate choice, multiple contacts across unsettled parapatric ranges, and few opportunities to evolve reinforcement given the much shorter time since first contact. While the prevalence of reinforcement in the wild is debated¹⁵, this mechanism is known to contribute to speciation in hylid frogs⁷⁰. For the same species pair, different genomic backgrounds (i.e. intraspecific lineages) may also affect patterns of
incompatibilities and admixture; here this could apply to the two *H. arborea/orientalis* contacts, which involve independent refugial lineages<sup>26,51</sup>. Therefore, replicate hybrid zone analyses appear necessary to fully appreciate where two taxa lay on the speciation continuum.

In this context, our findings have profound consequences for species delimitation. Phylogeographic lineages are increasingly considered in taxonomy⁸, and when possible, their specific or subspecific rank may be decided based on whether they are (partially) reproductively isolated, following the universal biological species concept³. Accordingly, wide transitions are interpreted as a lack of reproductive barriers, as opposed to narrow transitions³,⁷¹. An ad hoc way to quantify reproductive isolation from hybrid zones is to calculate selection against hybrids from cline width *w* and dispersal *σ*, with the equations of⁷². Without extrinsic selection (as expected for ecologically-similar cryptic species), the coefficient *s*\* corresponds to the fitness difference between populations from the center and the edge of the hybrid zone (under heterozygote disadvantage), given *w*≈2*σ*/√*s*\*⁷². In *Hyla*, taking a dispersal rate *σ* ~ 0.5 km/year (adapted from⁶⁶, who found average movements of 1.5 km over a three year capture-mark-recapture survey), and average ages at maturity and longevity of 2 and 6 years, respectively (reviewed by⁷³), the refugial contact between *H. arborea* and *H. orientalis* (*w* = 30 km³³,) reflect ten-fold higher estimates (*s*\* = 0.05) than the post-glacial contact between *H. arborea* and *H. molleri* (*s*\* = 0.005, from *w* = 98 km), despite similar genetic divergence. The extent of admixture between incipient species might thus misguide taxonomic conclusions if the history of contacts is not considered. When possible, the grey zone of speciation would be best inferred from refugial areas, where selection had time to operate and species transitions putatively reached an equilibrium. Reciprocally, narrow post-glacial hybrid zones would strongly indicate incipient speciation. To conclude, we recommend that taxonomists interpret introgression patterns at range margins in the light of two key timings of speciation, i.e. the time of divergence and the time since first contact.

Finally, the present results call for caution when assessing the hybridization potential of invasive species. Genetic pollution is a major issue of bioinvasions⁷⁴,⁷⁵, and the risk should not be underestimated even if the local and invasive species rarely hybridize in natural ranges. This is well-illustrated here by the introduction of *H. perrini* north of the Alps, which extensively admixed with local *H. arborea*⁵⁰, despite nearly absent introgression at their natural phylogeographic transition⁴⁹. In the synoptic fishes *Alosa pseufoharengus* and *A. aestivalis*, human-driven habitat fragmentation disrupted their well-established reproductive barriers, leading to whole-hybrid populations⁷⁶. If the diverged genomes of incipient species are partly incompatible, such “forced” hybridization can severely affect the sustainability of local taxa, and thus precipitate their extinctions.

## Conclusions

Combining new data with fifteen years of research on the W-Palearctic radiation of tree frogs (*Hyla*), we emphasized how admixture patterns at species transitions can drastically vary depending on their biogeographic history. For similar evolutionary divergence, species pairs form wider hybrid zone in post-glacial and introduced ranges, compared to regions spanning their putative glacial refugia. We thus conclude that both the divergence time and the time since contact significantly affect the extent of admixture between incipient species. Future studies should focus on identifying the proximate mechanisms, and especially test whether reinforcement evolved in long-term contacts. Informed taxonomic decisions are thus best taken from replicate hybrid zone analyses that account for the late-Quaternary dynamics of phylogeographic transitions.

## Data availability

The data of this study is available as supplementary information.

Received: 2 October 2019; Accepted: 11 March 2020;

Published online: 26 March 2020

## References

1. Avise, J. *Phylogeography: The history and formation of species*. Cambridge, (Harvard University Press) (2000).

2. Coyne, J. A. & Orr, H. A. *Speciation*. (Sinauer Associates) (2004).

3. Mayr, E. *Systematics and the origin of species, from the viewpoint of a zoologist*. (Harvard University Press), (1942).

4. Roux, C. et al. Shedding light on the grey zone of speciation along a continuum of genomic divergence. *PLoS Biol.* **14**, e20000234 (2016).

5. Singhal, S. & Moritz, C. Reproductive isolation between phylogeographic lineages scales with divergence. *Proc. R. Soc. Lond.* **280**, 20132246 (2013).

6. Arntzen, J. W., Wielstra, B. & Wallis, G. P. The modality of nine *Triturus* newt hybrid zones assessed with nuclear, mitochondrial and morphological data. *Biol. J. Linn. Soc.* **113**, 604–622 (2014).

7. Dufresnes, C. et al. Phylogeography of a cryptic speciation continuum in Eurasian spadefoot toads (*Pelobates*). *Mol. Ecol.* **28**, 3257–3270 (2019).

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

9. Barton, N. H. & Hewitt, G. M. Analysis of hybrid zones. *Annu. Rev. Ecol. Evol. Syst.* **16**, 113–148 (1985).

10. Dufkova, P., Macholan, M. & Pialek, J. Inference of selection and stochastic effects in the house mouse hybrid zone. *Evolution* **65**, 993–1010 (2011).

11. Beysard, N. & Heckel, G. Structure and dynamics of hybrid zones at different stages of speciation in the common vole (*Microtus arvalis*). *Mol. Ecol.* **23**, 673–687 (2014).

12. Wielstra, B. et al. A genomic footprint of hybrid zone movement in crested newts. *Evol. Lett.* **1**, 93–101 (2017).

13. Hewitt, G. Quaternary phylogeography: the roots of hybrid zones. *Genetica* **139**, 617–638 (2011).

14. Butlin, R. Speciation by reinforcement. *Trends Ecol. Evol.* **2**, 8–13 (1987).

15. Servedio, M. R. & Noor, M. A. F. The role of reinforcement in speciation: theory and data. *Annu. Rev. Ecol. Evol. Syst.* **34**, 339–364 (2003).

16. Saestre, G.-P., Kral, K., Bures, S. & Ims, R.-A. Dynamics of a clinal hybrid zone and a comparison with island hybrid zones of flycatchers (*Ficedula hypoleuca* and *F. albicollis*). *J. Zool.* **247**, 53–64 (1999).

17. Saestre, G.-P. et al. A sexually selected character displacement in flycatchers reinforces premating isolation. *Nature* **387**, 589–592 (1997).
18. Zeng, Y. F., Lia, W. J., Petit, R. J. & Zang, D. Y. Geographic variation in the structure of oak hybrid zones provides insights into the dynamics of speciation. *Mol. Ecol.* **20**, 4995–5011 (2011).

19. Croucher, P. J., Jones, R. R., Searle, J. B. & Oxford, G. S. Contrasting patterns of hybridization in large house spiders (*Tegenaria atrica* group, Agelenidae). *Evolution* **61**, 1622–1640 (2007).

20. Dufresnes, C. et al. Genomic evidence for cryptic speciation in tree frogs from the Apennine Peninsula, with description of *Hyla perrini* sp. nov. *Front. Ecol. Evol.* **6**, 144 (2018).

21. Stöck, M. et al. Cryptic diversity among Western Palearctic tree frogs: Postglacial range expansion, range limits, and secondary contacts of three European tree frog lineages (*Hyla arborea* group). *Mol. Phylogenet. Evol.* **65**, 1–9 (2012).

22. Drillon, O., Dufresnes, G., Perrin, N., Crochet, P.-A. & Dufresnes, C. Reaching the edge of the speciation continuum: hybridization between three sympatric species of tree frogs (*Hyla*). *Biol. J. Linn. Soc.* **126**, 743–750 (2019).

23. Gvoždík, V. et al. Speciation history and widespread introgression in the European short-call tree frogs (*Hyla arborea* sensu lato, *H. intermedia* and *H. sarda*). *Mol. Phylogenet. Evol.* **83**, 143–155 (2015).

24. Dufresnes, C., Brelsford, A., Béziers, P. & Perrin, N. Stronger transferability but lower variability in transcriptomic- than in anonymous microsatellites: evidence from Hyloid frogs. *Mol. Ecol. Resour.* **14**, 716–725 (2014).

25. Dufresnes, C. et al. Sex-chromosome homomorphy in Palearctic tree frogs results from both turnovers and X-Y recombination. *Mol. Biol. Evol.* **32**, 2328–2337 (2015).

26. Dufresnes, C. et al. Conservation phylogeography: does historical diversity contribute to regional vulnerability in European tree frogs (*Hyla arborea*)? *Mol. Ecol.* **22**, 5669–5684 (2013).

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

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

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

30. Rosenberg, N. A. Distract: a program for the graphical display of population structure. *Mol. Ecol. Resour.* **4**, 137–138 (2004).

31. Jombart, T. adegenet: a R package for the multivariate analysis of genetic markers. *Bioinformatics* **24**, 1403–1405 (2008).

32. Goudet, J. Hierstat, a package for R to compute and test hierarchical F-statistics. *Mol. Ecol. Notes* **5**, 184–186 (2005).

33. Dufresnes, C. et al. Timeframe of speciation inferred from secondary contact zones in the European tree frog radiation (*Hyla arborea* group). *BMC Evol. Biol.* **15**, 155 (2015).

34. Derryberry, E. P., Derryberry, G. E., Maley, J. M. & Brumfield, R. T. HZAR: hybrid zone analysis using an R software package. *Mol. Ecol. Resour.* **14**, 652–663 (2014).

35. Stöck, M. et al. Mitochondrial and nuclear phylogeny of circum-Mediterranean tree frogs from the *Hyla arborea* group. *Mol. Phylogenet. Evol.* **49**, 1019–1024 (2008).

36. Recuero, E., Iraola, A., Rubio, X., Mahordom, A. & García-París, M. Mitochondrial differentiation and biogeography of *Hyla meridionalis* (Anura: Hylidae): an unusual phylogeographical pattern. *J. Biogeogr.* **34**, 1207–1219 (2007).

37. Dufresnes, C. et al. Diversification and speciation in tree frogs from the Maghreb (*Hyla meridionalis* sensu lato), with description of a new African endemic. *Mol. Phylogenet. Evol.* **134**, 291–299 (2019).

38. Gvoždík, V., Moravec, J., Klütsch, C. & Kotlik, P. Phylogeography of the Middle Eastern tree frogs (*Hyla*, Hylidae, Amphibia) as inferred from nuclear and mitochondrial DNA variation, with a description of a new species. *Mol. Phylogenet. Evol.* **55**, 1146–1166 (2010).

39. Gül, S., Kutrup, B. & Özdemir, N. Patterns of distribution of tree frogs in Turkey based on molecular data. *Amphibia-Reptilia* **33**, 95–103 (2012).

40. Bisconti, R., Canestrelli, D., Colangelo, P. & Nascetti, G. Multiple lines of evidence for demographic and range expansion of a temperate species (*Hyla sarda*) during the last glaciation. *Mol. Ecol.* **20**, 5313–5327 (2011).

41. Bisconti, R., Canestrelli, D. & Nascetti, G. Genetic diversity and evolutionary history of the Tyrrhenian treefrog *Hyla sarda* (Anura: Hylidae): adding pieces to the puzzle of Corsica–Sardinia biota. *Biol. J. Linn. Soc.* **103**, 159–167 (2011).

42. Canestrelli, D., Cimmaruta, R. & Nascetti, G. Phylogeography and historical demography of the Italian treefrog, *Hyla intermedia*, reveals multiple refugia, population expansions and secondary contacts within peninsular Italy. *Mol. Ecol.* **16**, 4808–4821 (2007).

43. Canestrelli, D., Verardi, A. & Nascetti, G. Genetic differentiation and history of populations of the Italian treefrog *Hyla intermedia*: lack of concordance between mitochondrial and nuclear makers. *Genetica* **130**, 241–255 (2007).

44. Dufresnes, C. et al. Sex-chromosome differentiation parallels postglacial range expansion in European tree frogs (*Hyla arborea*). *Evolution* **68**, 3445–3456 (2014).

45. Barth, A. et al. Mitochondrial uniformity in populations of the treefrog *Hyla molleri* across the Iberian Peninsula. *Amphibia-Reptilia* **32**, 557–564 (2011).

46. Sánchez-Montes, G., Recuero, E., Marcia Barbosa, A. & Martinez-Solano, I. Complementing the Pleistocene biogeography of European amphibians: Testimony from a southern Atlantic species. *J. Biogeogr.* **46**, 568–583 (2019).

47. Dufresnes, C. et al. Evolutionary melting pots: a biodiversity hotspot shaped by ring diversifications around the Black Sea in the Eastern tree frog (*Hyla orientalis*). *Mol. Ecol.* **25**, 4285–4300 (2016).

48. Dufresnes, C., Mazepa, G., Jablonski, D., Sadek, R. A. & Litvinchuk, S. N. A river runs through it: tree frogs genomics supports the Dead Sea Rift as a rare phylogeographic break. *Biol. J. Linn. Soc.* **128**, 130–137 (2019).

49. Verardi, A., Canestrelli, D. & Nascetti, G. Nuclear and mitochondrial patterns of introgression between the parapatric European tree frogs *Hyla arborea* and *H. intermedia*. *Ann. Zool. Fenn.* **46**, 247–258 (2009).

50. Dufresnes, C., Dubey, S., Ghali, K., Canestrelli, D. & Perrin, N. Introgressive hybridization of threatened European tree frogs (*Hyla arborea*) by introduced *H. intermedia* in Western Switzerland. *Conserv. Genet.* **16**, 1507–1513 (2015).

51. Dufresnes, C. et al. Empirical evidence for large X-effects in animals with undifferentiated sex chromosomes. *Sci. Rep.* **6**, 21029 (2016).

52. Smith, S. A., Stephens, P. R. & Wiens, J. J. Replicate patterns of species richness, historical biogeography and phylogeny in Holarctic treefrogs. *Evolution* **59**, 2433–2450 (2005).

53. Li, J.-T. et al. Amphibians crossing the Bering Land Bridge: evidence from Holarctic treefrogs (*Hyla*, Hylidae, Anura). *Mol. Phylogenet. Evol.* **87**, 80–90 (2015).

54. Kotsakis, T. Osservazioni sui vertebrati quaternary della sardegna. *Boll. Soc. geol. ital.* **99**, 151–165 (1980).

55. Phillips, S., Anderson, R. & Schapire, R. Maximum entropy modeling of species geographic distributions. *Ecol. Model.* **190**, 231–259 (2006).

56. Osorio-Olvera, L., Barve, V., Barve, N., Soberón, J., & Falconi, M. ntbox: From getting biodiversity data to evaluating species distribution models in a friendly GUI environment. https://github.com/luismurao/ntbox (2018).

57. Watanabe, S. et al. MIROC-ESM 2010: model description and basic results of CMIP5-20c3 m experiments. *Geosci. Model Dev.* **4**, 845–872 (2011).

58. Warren, D. L., Glor, R. E. & Turelli, M. ENMTOOLS: a toolbox for comparative studies of environmental niche models. *Ecography* **33**, 607–611 (2010).

59. Cobos, M. E., Peterson, A. T., Barve, N. & Osorio-Olvera, L. kuenm: an R package for detailed development of ecological niche models using Maxent. *PeerJ* **7**, e6281 (2019).
60. Becker, D., Verheul, J., Zickel, M., & Willmes C. LGM paleoenvironment of Europe - Map. CRC 806 Database. https://doi.org/10.5880/SFB806.15 (2015).

61. Zickel, M. Paleocoastline 130m below mean sea level. CRC 806 Database, https://crc806db.uni-koeln.de/layer/show/324/ (2016).

62. Berroneau, M. *Atlas des amphibians et reptiles d'Aquitaine*. (Cistude Nature), (2014).

63. Dufresnes, C. et al. Inferring the degree of incipient speciation in secondary contact zones of closely related lineages of Palearctic green toads (*Bufo viridis* subgroup). *Heredity* **113**, 9–20 (2014).

64. Feng, Y.-J. et al. Phylogenomics reveals rapid, simultaneous diversification of three major clades of Gondwanan frogs at the Cretaceous-Paleogene boundary. *PNAS* **114**, E5864–E5870 (2017).

65. Dubey, S., Lavanchy, G., Thiébaud, J. & Dufresnes, C. Herps without border: a new newt case and a review of transalpine alien introductions in Western Europe. *Amphibia-Reptilia* **40**, 13–27 (2019).

66. Vos, C. C., Ter Braak, C. J. F. & Nieuwenhuizen, W. Incidence function modelling and conservation of the tree frog *Hyla arborea* in the Netherlands. *Ecol. Bull.* **48**, 165–180 (2000).

67. Trochet, A. et al. A database of life-history traits of European amphibians. *Biodivers. Data J.* **2**, e4123 (2014).

68. Peter, B. M. & Slatkin, M. Detecting range expansions from genetic data. *Evolution* **67**, 3274–3289 (2013).

69. Currat, M., Ruedi, M., Petit, R. J. & Excoffier, L. The hidden side of invasions: massive introgression by local genes. *Evolution* **62**, 1908–1920 (2008).

70. Hoskin, C. J., Higgie, M., McDonald, K. R. & Moritz, C. Reinforcement drives rapid allopatric speciation. *Nature* **437**, 1353–1356 (2005).

71. Dufresnes, C. et al. Integrating hybrid zone analyses in species delimitation: lessons from two anuran radiations of the Western Mediterranean. *Heredity* **124**, 423–438 (2020).

72. Barton, N., & Gale, K. S. Genetic analysis of hybrid zones in *Hybrid zones and the evolutionary process* (ed. R., Harrison) 13–45 (Oxford University Press), (1993).

73. Özdemir, N. et al. Variation in body size and age structure among three Turkish populations of the treefrog *Hyla arborea*. *Amphibia-Reptilia* **33**, 25–35 (2012).

74. Petit, R. J. Biological invasion at the gene level. *Divers. Distrib.* **10**, 159–165 (2004).

75. Gaskin, J. T. The role of hybridization in facilitating tree invasion. *AoB Plants* **9**, plw079 (2017).

76. Hasselman, D. et al. Human disturbance causes the formation of a hybrid swarm between two naturally sympatric fish species. *Mol. Ecol.* **23**, 1137–1152 (2014).

77. Dufresnes, C. *Amphibians of Europe, North Africa and the Middle East*. (Bloomsbury), (2019)

## Acknowledgements

We thank J Jacquiéry and O Drillon for providing samples, as well as M Stöck and G Dufresnes for their help on the field. V Gvoždík and two anonymous reviewers provided useful comments on this manuscript. This study was funded by a grant from the Swiss National Science Foundation (n° 31003A_166323) to NP and a PhD fellowship from the Faculty of Biology and Medicine of the University of Lausanne to CD; SLN was supported by a grant from the Russian Foundation for Basic Research (n° RFBR.20-04-00918).

## Author contributions

C.D. and N.P. designed the study. C.D. and M.B. conducted fieldwork. C.D. conducted labwork and analyzed the data with S.D. and S.N.L. C.D. wrote the manuscript, improved by all co-authors.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41598-020-62382-4.

Correspondence and requests for materials should be addressed to C.D.

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

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

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

© The Author(s) 2020