Hybrid zone genomics supports candidate species in Iberian Alytes obstetricans

. While estimates of genetic divergence are increasingly used in molecular taxonomy, hybrid zone analyses can provide decisive evidence for evaluating candidate species. Applying a population genomic approach (RAD-sequencing) to a ﬁne-scale transect sampling, we analyzed the transition between two Iberian subspecies of the common midwife toad ( Alytes obstetricans almogavarii and A. o. pertinax ) in Catalonia (northeastern Spain), which putatively diverged since the Plio-Pleistocene. Their hybrid zone was remarkably narrow, with extensive admixture restricted to a single locality (close to Tarragona), and congruent allele frequency clines for the mitochondrial (13 km wide) and the average nuclear genomes (16 km wide). We also ﬁtted clines independently for 89 taxon-diagnostic SNPs: most of them behave like the nuclear background, but a subset (13%) is completely impermeable to gene ﬂow and might be linked to barrier loci involved in hybrid incompatibilities. Assuming that midwife toads are able to disperse in the area of contact, we conclude that these taxa experience partial reproductive isolation and represent incipient species, i.e. Alytes almogavarii and Alytes obstetricans . Interestingly, their evolutionary age and mitochondrial divergence fall below the thresholds proposed in molecular systematics studies, emphasizing the difﬁculty of predicting the outcome of secondary contacts between young lineages entering the grey zone of speciation.


Introduction
The delimitation of cryptic species by phylogeographers is inherently challenged by the continuous and dynamic nature of speciation processes (Avise, 2000;Heethoff, 2018;Struck et al., 2018;Chenuil et al., 2019). In the era of the "molecularisation of taxonomy" (Lee, 2004), evolutionary biologists may never agree on how much reproductive isolation and genetic differentiation is required to consider diverging populations as incipient species (Coyne and Orr, 2004;De Queiroz, 2007;Padial et al., 2010). In particular, hybridizing taxa, as well as allopatric taxa of equivalent divergence, fall within a socalled grey zone of speciation, and thus remain at the heart of taxonomic disputes (Roux et al., 2016).
In the wake of recent attempts to relate the evolutionary divergence of phylogeographic lineages with levels of introgression across hybrid zones (Singhal and Moritz, 2013;Dufresnes et al., 2019a), the concept of "speciation thresholds" might circumvent some practical issues (Fouquet et al., 2007;Vieites et al., 2009). As theoretically and empirically supported, reproductive isolation tends to increase rapidly due to the accumulation of barrier loci after a certain degree of genetic divergence, above which the probability to merge back and reverse the speciation process is considerably reduced (Barton, 1983;Orr, 1995;Gourbiere and Mallet, 2009;Singhal and Moritz, 2013). Hence, if properly characterized and implemented, this checkpoint should facilitate the delimitation of incipient species.
Palearctic amphibians are a good system to assess the performance of such speciation thresholds for integrative systematics, in respect to intrinsic limitations such as the choice of markers (mtDNA, nuclear) and metrics to apply (% of divergence, divergence time). In am-phibians, while some influential studies recommended 3% differentiation at sequences of the mitochondrial 16S gene (Malone and Fontenot, 2008;Vieites et al., 2009), recent comparative surveys suggested a divergence time of 3 My to underlie incipient speciation among European anurans. This age corresponds to the youngest candidate species in several well-studied radiations, while younger lineages usually merge upon secondary contact (Dufresnes et al., 2019a and references therein). However, molecular clocks are specific to the phylogenetic framework considered, and the speciation clock may tick at different paces depending on lineages. Empirical data from additional species groups are thus needed to fully appreciate the progress of reproductive isolation in space and time, and refine the thresholds used in systematics.
Here we focused on the midwife toad Alytes obstetricans, a small anuran that colonized parts of Western Europe from the Iberian Peninsula. This taxon bears the hallmark of the Plio-Pleistocene climatic fluctuations, and previous phylogeographic analyses identified six deeply-diverged mitochondrial lineages distributed across Central Spain (A), N-Spain/W-Europe (B), Galicia and N-Portugal (C), Central Portugal/Sistema Central mountains in Spain (D), the Spanish Pyrenees (E) and Catalonia (F) . Most of these clades feature nuclear differentiation, some corresponding to taxonomic subdivisions (Arntzen and García-París, 1995;Martínez-Solano et al., 2004;Gonçalves et al., 2007;Maia-Carvalho et al., 2014, 2018Gonçalves et al., 2015). The most divergent populations (mtDNA lineages E-F), representing the subspecies A. o. almogavarii, are of particular interest. Dated to the Mio-Pliocene (Martínez-Solano et al., 2004) or to the Plio-Pleistocene , this taxon features the highest allozyme distances to all other A. obstetricans subspecies (Nei's D = 0.17; Arntzen and García-París, 1995; see also García-París, 1995). While its placement on the nuclear tree is not well resolved (Gonçalves et Gonçalves et al., 2015) and amount of 16S divergence (0.99-1.38% in Martínez-Solano et al., 2004) are both below the proposed thresholds of speciation, yet the apparent lack of genetic introgression suggests potential reproductive isolation between A. o. almogavarii and its congeners, which could indicate an unusually young speciation event.
To elucidate whether A. o. almogavarii represents an incipient species, we thoroughly analyzed the transition between Alytes obstetricans, lineages F (almogavarii) and A (pertinax) along the Catalonian coast in northeastern Spain Maia-Carvalho et al., 2018), based on transect sampling and genomic resources (RAD-sequencing). Interestingly, this phylogeographic break is also found in a related amphibian, the common parsley frog Pelodytes punctatus, where subspecies P. p. punctatus and P. p. hespericus introgress over hundreds of kilometers (CD and IMS, unpublished data) despite equivalent divergence time (Díaz-Rodríguez et al., 2015. Assuming similar dispersal capacity in Alytes and Pelodytes (both share mostly terrestrial habits and are often found in syntopy in Mediterranean habitats), we predict widespread admixture across the A. o. almogavarii/pertinax transition if, like P. p. punctatus/hespericus, they did not evolve significant reproductive isolation during their short history of divergence. On the contrary, a narrow contact zone with limited introgression is expected if they reached the threshold of divergence above which reproductive isolation might have quickly accumulated.

Sampling and data generation
We combined 38 samples from a previous study  with 73 specimens that were newly-collected during the summer of 2017 along a geographic transect in Catalonia ( fig. 1, supplementary table S1). Tissue samples consisted of tail tips from subsequently released tadpoles, preserved in 96% ethanol, from which DNA was extracted using the Qiagen Biosprint Robotic workstation.
The 111 samples (18 populations) were barcoded for mitochondrial DNA using custom-designed primers Aly.ND4. F3 (5 -AAACCCCCTGAAGTTTTTCAGG-3 ) and Aly. ND4.R3 (5 -TCTTGGTGGGCAAGAAGGTT-3 ), which amplify ∼500 bp of the mitochondrial gene ND4. These samples were then included in a RAD genomic library built with enzymes SbfI and MseI, following the protocol detailed in Brelsford, Dufresnes and Perrin (2016). The library was sequenced on two Illumina lanes and processed bioinformatically using the STACKS pipeline (v. 1.48), with default parameters (Catchen et al., 2013). We then filtered tags present in all populations (−p 18), in 80% of individuals of each population (−r 0.8), but with a minimum allele frequency of 0.05 (-min_maf 0.05) and a maximum observed heterozygosity of 0.75 (-max_obs_het 0.75) across all individuals. This retained 433 SNPs that were informative across populations, with little missing data (97.4% of matrix completeness), while limiting paralog overmerging.

Genetic analyses
Patterns of genetic structure among Catalonian A. o. almogavarii and A. o. pertinax samples were assessed from the SNP dataset using the clustering algorithm of STRUCTURE (Pritchard et al., 2000). Five replicates from K = 1 to 5 were ran, each including 100 000 iterations after a burnin of 10 000. We also performed a Principal Component Analysis (PCA) on individual allele frequencies with the R package adegenet (Jombart, 2008).
Furthermore, we measured the geographic transition between the two gene pools by fitting sigmoid clines on diagnostic allele frequency data from 17 populations along the transect (all but loc. 5), using the R package hzar (Derryberry et al., 2014). Clines were inferred from three genetic datasets: (1) the average population ancestry (Q) obtained by STRUCTURE with K = 2 (genome-average cline); (2) the mitotype frequencies (mitochondrial cline); (3) the frequencies of 89 SNPs fixed between taxa (all from different RAD tags), i.e. with F st = 1 between loc. 1-3 (A. o. almogavarii) and loc. 18 (A. o. pertinax). We tested clines from two (width w and center c) to six parameters (accounting for introgression tails), and reported those best fitting the data according to the AIC criterion. Finally, we estimated intrinsic selection against hybrids s * (the fitness difference between populations from the center and the edge of the contact zone) from the cline width w, with respect to the dispersal rate per generation σ (p 16 in Barton and Gale, 1993). Since little is known about the dispersal behavior in midwife toads (Tobler, Garner and Schmidt, 2013), here we inferred s * for a range of σ values from 0.1 to 4 km/generation.

Results
Our genomic data (433 SNPs) perfectly differentiated between the gene pools of A. o. almogavarii (coded yellow) and A. o. pertinax (coded blue) along the Catalonian coast ( fig. 1). The best-clustering solution recovered by STRUCTURE involved two clusters (K = 2, K = 1121.1; fig. 1). Accordingly, almost half of the total genetic variance was summarized by the first axis of the PCA and discriminates the two taxa (supplementary fig. S1).
Based on our analyses, these two clusters feature highly restricted genetic introgression. First, STRUCTURE (K = 2) identified a single locality with clear nuclear admixture (loc. 14, where haplotypes from both mtDNA lineages are present), while the neighboring loc. 13 (inhabited by A. o. almogavarii) features two toads with Q scores below 0.95, and two others carrying pertinax mtDNA ( fig. 1). All other samples received ancestry coefficients above 0.95 ( fig. 1). These results are confirmed by the PCA, on which the few introgressed individuals stand out by intermediate values along the first axis (supplementary fig. S1).
Second, the cline analyses inferred narrow transitions for all markers ( fig. 2). Cline models with two parameters (width w and center c) were always preferred (lowest AIC scores), and provided remarkably congruent results for mtDNA (w = 13.5 km, c = 86.1 km), the nuclear genome average (w = 16.2 km, c = 91.1 km) and the 89 diagnostic nuclear loci (w = 20.4 ± 8.5 km and c = 90.2 ± 2.5 km). For most markers, the transition spanned exclusively between our loc. 13 and 15 (located 30 km apart, fig. 1), near the villages of Valls and Montblanc. While the majority of loci (74%) broadly follow the genome average, two subsets feature narrower clines (13% with w tending towards 0 and 12% with w = 2-4 km), as best visualized from the derived s * values ( fig. 2B) and the distribution of cline widths (supplementary fig. S2). Moreover, the hybrid zone may be slightly asymmetric, with the center of the mtDNA cline and of many nuclear SNPs being shifted by a few kilometers on the almogavarii side compared to the genome average cline ( fig. S2).

Discussion
The narrow mitochondrial and nuclear transitions suggest that A. o. pertinax and A. o. almogavarii form a classic tension zone in Catalonia. From the average cline estimates, assuming a dispersal rate of 1.5 km/generation is enough to get s * > 0.03 ( fig. 2B), i.e. the values obtained in well-studied hybrid zones involving cryptic anuran species (e.g. H. arborea/orientalis, Dufresnes et al., 2015;Pelobates fuscus/vespertinus, Dufresnes et al., 2019a; Pelodytes ibericus/atlanticus) (CD and IMS, unpublished data). But it would take up to 3.5 km/generation to reach the strong s * inferred for deeper-diverged, ecomorphologically differentiated taxa (e.g. s * = 0.21 in Bombina variegata/bombina, Barton and Gale, 1993) (fig. 2B). Accounting for a generation time of seven years in Alytes (computed from the average age at sexual maturity and lifespan reported by Márquez, Esteban and Castanet, 1997), 1.5-3.5 km/generation thus corresponds to yearly movements of only 0.2-0.5 km.
Such (relatively short) distances seem realistic for Alytes in the Catalonian landscape. Midwife toads are adapted to xeric environments and accommodate to tiny water points (e.g. droughts), which are visited just a few days per year by the adult males, when they deposit the offspring. Movements up to 1.5 km have been observed in the related Alytes dickhilleni endemic to SE-Iberia (Bosch and González Miras, 2012), a region drier than the Catalonian inlands. Note that habitat fragmentation can significantly limit dispersal, even between nearby subpopulations : Tobler, Garner and Schmidt (2013) reported genetic differentiation between sites located <2 km apart in A. obstetricans from the heavily impacted Swiss plateau. Here however, a stream network runs along the center of the almogavarii/pertinax transition (providing breeding sites, e.g. loc. 14), a landscape element predicted to promote connectivity and persistence among amphibian populations (Grant et al., 2010).
These observations suggest that the limited introgression reflects partial reproductive isolation rather than limited opportunities for dispersal. This is further supported by the heterogeneity of cline widths among our diagnostic markers: the outliers featuring very steep clines (corresponding to strong selective values) could be linked to barrier loci, i.e. involved in hybrid incompatibilities. Theory predicts that most of the genome remains permeable to gene flow in the early stages of incipient speciation, but that post-zygotic isolation arises as barrier loci become recruited and limit introgression across entire genomic regions, due to linkage disequilibrium (Barton, 1983;Orr, 1995;Gourbiere and Mallet, 2009). Surveys across additional almogavarii/pertinax contacts would allow testing whether the same clines remain sharp in replicate transects, thus confirming a selective role. Comparative analyses of amphibian hybrid zones shall provide insights on the genomic architecture of introgression at different stages along the speciation continuum, especially to assess how barrier loci accumulate with genetic divergence.
Our results thus call for a taxonomic revision of Alytes obstetricans. We hereby consider the Catalonian Midwife Toad Alytes almogavarii Arntzen, García-París, 1995 as a distinct incipient species, which was previously considered as a subspecies of the Common Midwife Toad Alytes obstetricans (Laurenti, 1768). The biodiversity of the Spanish and French herpetofauna should thus be re-evaluated to account for this regional endemic, especially for conservation planning and red lists. In France, A. almogavarii appears restricted to the eastern Pyrenees (genetically confirmed in the departments of Aude and Pyrénées-Orientales, and perhaps extending to Ariège; Gonçalves et al., 2015;Maia-Carvalho et al., 2018), so it is extremely localized in the country. Downloaded from Brill.com05/06/2020 12:11:53AM via free access Are Alytes obstetricans and A. almogavarii cryptic species? In the original description, Arntzen and García-París (1995) compared external characters and found high morphological similarity between all Iberian Alytes, with the exception of the early-diverged A. cisternasii and the habitat-specialist A. muletensis (see also Sanchiz, 1984). Based on skull osteology however, A. almogavarii shares some character states with species in subgenus Baleaphryne (including A. muletensis, A. maurus, and A. dickhilleni) rather than with the other A. obstetricans subspecies (Martínez-Solano et al., 2004). Furthermore, bioacoustic investigations suggested inter-population variation in breeding calls between Iberian taxa, including longer notes in A. almogavarii compared to A. obstetricans subspecies (Márquez and Bosch, 1995). To our knowledge, the tadpoles of the different subspecies have not been formally compared. We did observe lighter and more contrasted coloration in tadpoles of A. almogavarii than in those of A. o. pertinax during the fieldwork of this study (IMS & CD, pers. obs.), but this could simply relate to their micro-habitats (streams vs still waters), which are known to influence coloration in Alytes (Polo-Cavia et al., 2016). The new taxonomic status of A. almogavarii should hopefully motivate comparative assessments, especially to document whether the differing breeding calls may induce some pre-mating isolation with A. obstetricans, and thus contribute to the limited hybridization observed.
The tension zone between A. obstetricans and A. almogavarii sharply contrasts with the wide hybrid zone (>150 km) found for the parsley frogs Pelodytes punctatus punctatus/hespericus in Catalonia (Díaz-Rodríguez et al., 2017;CD and IMS, unpublished data), despite putatively similar divergences (∼2.5 My for Alytes; ∼2.8 My for Pelodytes). While massive introgression is the rule rather than the exception for Plio-Pleistocene anuran lineages (e.g. Dufresnes et al., 2019a), the opposite patterns between Alytes and Pelodytes emphasize that the outcomes of secondary contacts are hard to predict for pairs of taxa falling within the grey zone of speciation. First, admixture across hybrid zones depends on a myriad of extrinsic factors, and may thus drastically vary between replicate contacts (e.g. Croucher et al., 2007). For instance, asymmetric introgression across the genome may arise due to hybrid zone movement, a situation increasingly reported in the wild (e.g. van Riemsdijk et al., 2019). Here, some clines appear slightly shifted on the almogavarii side (supplementary fig. S2), but it could also be a methodological bias due to our ∼20 km sampling gap on the pertinax side of the contact (between loc. 14 and 15). Second, the random accumulation of genetic incompatibilities can lead to different levels of postzygotic isolation between independent pairs of lineages, even if they share equivalent divergences (e.g. Arntzen et al., 2017). The latter could be directly assessed in Alytes by investigating the transition between A. almogavarii and A. o. obstetricans (the sister clade of A. o. pertinax), most likely extending along the French Pyrenees . More generally, these aspects are crucial to consider when applying thresholds of speciation: here A. almogavarii falls below both the sequence and time thresholds set up by previous studies (see Introduction), which might explain why it was not previously considered a distinct species. Although thresholds are useful for species delimitation, especially for completely allopatric taxa (e.g. insular species), direct evidences from targeted studies of secondary contacts thus appear essential to make the best informed taxonomic revisions.
Alternatively, in this case, it is conceivable that the mitochondrial inferences  misplaced A. almogavarii in the phylogeny and underestimated its evolutionary age. Relying on transformed-branch length estimates, Martínez-Solano et al. (2004) rather advocated for a Mio-Pliocene split for this clade. Moreover, some nuclear gene trees group A. almogavarii with the Baleaphryne subgenus (Gonçalves et al., 2007; but see Maia-Carvalho et al., 2014), a topology consistent with their os-teological similarities (Martínez-Solano et al., 2004). Could A. almogavarii be much older than previously assumed, but lost its mtDNA during Plio-Pleistocene episodes of hybridization? Complete mitochondrial replacement appears to be common in Palearctic amphibians (e.g. Zieliński et al., 2013), eventually concealing speciation events (the super-cryptic species concept, Dufresnes et al., 2019b). Phylogenomics of the entire Alytes genus should hopefully clarify this situation and bring complementary insights on the evolutionary history of A. almogavarii. Finally, resolving the complex taxonomy of this group will benefit from genomic analyses targeting additional transitions among the polytypic A. obstetricans, and notably within the genetically-structured A. o. obstetricans (Maia-Carvalho et al., 2018).