# Spatial genetic structure in European marbled newts revealed with target enrichment by sequence capture

Christos Kazilas, Christophe Dufresnes, James France, Konstantinos

Kalaentzis, Iñigo Martínez-Solano, Manon C de Visser, Jan W Arntzen, Ben

Wielstra

► To cite this version:

Christos Kazilas, Christophe Dufresnes, James France, Konstantinos Kalaentzis, Iñigo Martínez-Solano, et al. Spatial genetic structure in European marbled newts revealed with target enrichment by sequence capture. Molecular Phylogenetics and Evolution, 2024, 194, pp.108043. 10.1016/j.ympev.2024.108043. hal-04547713

HAL Id: hal-04547713

https://hal.science/hal-04547713

Submitted on 16 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

# Spatial genetic structure in European marbled newts revealed with target enrichment by sequence capture

Christos Kazilas<sup>a, b, *</sup>, Christophe Dufresnes<sup>c, d</sup>, James France<sup>a, b</sup>, Konstantinos Kalaentzis<sup>a, b</sup>,
Iñigo Martínez-Solano<sup>e</sup>, Manon C. de Visser<sup>a, b</sup>, Jan W. Arntzen<sup>a, b, 1</sup>, Ben Wielstra<sup>a, b, 1</sup>

## ARTICLE INFO

### Keywords:

Amphibians
Iberian Peninsula
Phylogenomics
Salamander
*Triturus*

## ABSTRACT

European marbled newts come in two species that have abutting ranges. The northern species, *Triturus marmoratus*, is found in France and the northern part of the Iberian Peninsula, whereas the southern species, *T. pygmaeus*, is found in the southwestern corner of the Iberian Peninsula. We study the intraspecific genetic differentiation of the group because morphological data show geographical variation and because the Iberian Peninsula is a recognized center of speciation and intraspecific genetic diversity for all kinds of organisms, amphibians included. We use target enrichment by sequence capture to generate c. 7 k nuclear DNA markers. We observe limited genetic exchange between the species, which confirms their distinctiveness. Both species show substantial genetic structuring that is only in part mirrored by morphological variation. Genetically differentiated groups are found in the south (*T. marmoratus*) and west (*T. pygmaeus*) of the species ranges. Our observations highlight the position of the Iberian Peninsula as a hotspot for genetic differentiation.

## 1. Introduction

A fundamental goal of evolutionary biogeography is to determine how long-term environmental processes influence the distribution and diversification patterns of taxa, especially within groups of closely related species composed of multiple evolutionary lineages (Avise, 2000). The Mediterranean region has been the subject of intensive phylogeographic surveys over the past decades, presumably triggered by the high amount of biodiversity and endemism that it hosts (Cuttelod et al., 2009; Médail and Quezél, 1999; Myers et al., 2000). The species richness of the Mediterranean Basin is thought to be the product of its ecological heterogeneity and complex geological history. As a consequence, the entire area is recognized as one of the world's major biodiversity hotspots (Médail and Myers, 2004).

In shallower times, the Mediterranean region also functioned as a glacial refugium, which allowed populations of temperate taxa to survive locally during the cold spells of the Quaternary Ice Age, when

environmental conditions further north became intolerable (Hewitt, 2011a). Populations in glacial refugia are typically characterized by higher genetic variation (Hewitt, 2000), both due to a more stable demographic history and the cyclical dynamics of lineage allopatry and fusion. Furthermore, glacial refugia appear more complex than initially thought, and harbored multiple, fragmented areas of suitable habitat: a concept called “refugia-within-refugia” (Abellán and Svenning, 2014; Gómez and Lunt, 2007). Accordingly, phylogeographical studies of Mediterranean taxa often reveal distinct evolutionary lineages. At contact zones these lineages may show wide zones of genetic intergradation (suggesting genetic isolation is limited) to narrow hybrid zones with little to no introgression (suggesting distinct species are involved) (Hewitt, 2011b).

Amphibians have played an important role in exposing the patterns and processes underlying Mediterranean biodiversity (Ehl et al., 2019; Oosterbroek and Arntzen, 1992). As terrestrial ectotherms, amphibians are particularly sensitive to past and present climatic fluctuations

<footnote><sup>a</sup> Naturalis Biodiversity Center, P.O. Box 9517, 2300 RA Leiden, The Netherlands</footnote>

<footnote><sup>b</sup> Institute of Biology Leiden, Leiden University, P.O. Box 9505, 2300 RA Leiden, The Netherlands</footnote>

<footnote><sup>c</sup> LASER, College of Biology and the Environment, Nanjing Forestry University, Nanjing, People's Republic of China</footnote>

<footnote><sup>d</sup> Institut de Systématique, Evolution, Biodiversité (ISYEB), Muséum national d'Histoire naturelle, CNRS, Sorbonne Université, EPHE, Université des Antilles, Paris, France</footnote>

<footnote><sup>e</sup> Museo Nacional de Ciencias Naturales, MNCN-CSIC, c/ José Gutiérrez Abascal 2, 28006 Madrid, Spain</footnote>

<footnote>* Corresponding author at: Zoodochou Pigis 15, 38222 Volos, Greece.</footnote>

<footnote>E-mail address: ckazilas@gmail.com (C. Kazilas).</footnote>

<footnote n="1">¹ these authors contributed equally.</footnote>
<figure><img src="image_2.png" /></figure>

# Spatial genetic structure in European marbled newts revealed with target enrichment by sequence capture

Christos Kazilas<sup>a, b, *</sup>, Christophe Dufresnes<sup>c, d</sup>, James France<sup>a, b</sup>, Konstantinos Kalaentzis<sup>a, b</sup>,
Iñigo Martínez-Solano<sup>e</sup>, Manon C. de Visser<sup>a, b</sup>, Jan W. Arntzen<sup>a, b, 1</sup>, Ben Wielstra<sup>a, b, 1</sup>

## ARTICLE INFO

### Keywords:

Amphibians
Iberian Peninsula
Phylogenomics
Salamander
*Triturus*

## ABSTRACT

European marbled newts come in two species that have abutting ranges. The northern species, *Triturus marmoratus*, is found in France and the northern part of the Iberian Peninsula, whereas the southern species, *T. pygmaeus*, is found in the southwestern corner of the Iberian Peninsula. We study the intraspecific genetic differentiation of the group because morphological data show geographical variation and because the Iberian Peninsula is a recognized center of speciation and intraspecific genetic diversity for all kinds of organisms, amphibians included. We use target enrichment by sequence capture to generate c. 7 k nuclear DNA markers. We observe limited genetic exchange between the species, which confirms their distinctiveness. Both species show substantial genetic structuring that is only in part mirrored by morphological variation. Genetically differentiated groups are found in the south (*T. marmoratus*) and west (*T. pygmaeus*) of the species ranges. Our observations highlight the position of the Iberian Peninsula as a hotspot for genetic differentiation.

## 1. Introduction

A fundamental goal of evolutionary biogeography is to determine how long-term environmental processes influence the distribution and diversification patterns of taxa, especially within groups of closely related species composed of multiple evolutionary lineages (Avise, 2000). The Mediterranean region has been the subject of intensive phylogeographic surveys over the past decades, presumably triggered by the high amount of biodiversity and endemism that it hosts (Cuttelod et al., 2009; Médail and Quezél, 1999; Myers et al., 2000). The species richness of the Mediterranean Basin is thought to be the product of its ecological heterogeneity and complex geological history. As a consequence, the entire area is recognized as one of the world's major biodiversity hotspots (Médail and Myers, 2004).

In shallower times, the Mediterranean region also functioned as a glacial refugium, which allowed populations of temperate taxa to survive locally during the cold spells of the Quaternary Ice Age, when

environmental conditions further north became intolerable (Hewitt, 2011a). Populations in glacial refugia are typically characterized by higher genetic variation (Hewitt, 2000), both due to a more stable demographic history and the cyclical dynamics of lineage allopatry and fusion. Furthermore, glacial refugia appear more complex than initially thought, and harbored multiple, fragmented areas of suitable habitat: a concept called “refugia-within-refugia” (Abellán and Svenning, 2014; Gómez and Lunt, 2007). Accordingly, phylogeographical studies of Mediterranean taxa often reveal distinct evolutionary lineages. At contact zones these lineages may show wide zones of genetic intergradation (suggesting genetic isolation is limited) to narrow hybrid zones with little to no introgression (suggesting distinct species are involved) (Hewitt, 2011b).

Amphibians have played an important role in exposing the patterns and processes underlying Mediterranean biodiversity (Ehl et al., 2019; Oosterbroek and Arntzen, 1992). As terrestrial ectotherms, amphibians are particularly sensitive to past and present climatic fluctuations

<footnote><sup>a</sup> Naturalis Biodiversity Center, P.O. Box 9517, 2300 RA Leiden, The Netherlands</footnote>

<footnote><sup>b</sup> Institute of Biology Leiden, Leiden University, P.O. Box 9505, 2300 RA Leiden, The Netherlands</footnote>

<footnote><sup>c</sup> LASER, College of Biology and the Environment, Nanjing Forestry University, Nanjing, People's Republic of China</footnote>

<footnote><sup>d</sup> Institut de Systématique, Evolution, Biodiversité (ISYEB), Muséum national d'Histoire naturelle, CNRS, Sorbonne Université, EPHE, Université des Antilles, Paris, France</footnote>

<footnote><sup>e</sup> Museo Nacional de Ciencias Naturales, MNCN-CSIC, c/ José Gutiérrez Abascal 2, 28006 Madrid, Spain</footnote>

<footnote>* Corresponding author at: Zoodochou Pigis 15, 38222 Volos, Greece.</footnote>

<footnote>E-mail address: ckazilas@gmail.com (C. Kazilas).</footnote>

<footnote n="1">¹ these authors contributed equally.</footnote>
(Zeisset and Beebee, 2008) as well as topographic barriers to dispersal. Mediterranean amphibians regularly represent groups of closely related species that are themselves characterized by high intraspecific genetic variability, with evolutionary lineages of varying ages (Ambu et al., 2023; Dufresnes et al., 2020a; Fijarczyk et al., 2011; Martinez-Solano et al., 2006; Pabijan et al., 2015; Recuero et al., 2012; Vences et al., 2013; Wielstra et al., 2013). Typically, relatively young evolutionary lineages show wide zones of genetic admixture between them, while introgression between relatively old evolutionary lineages is restricted (Arntzen et al., 2014; Dufresnes et al., 2023; Kalaentzis et al., 2023; Wielstra et al., 2018, 2017). Many of these evolutionary lineages have recently been recognized or confirmed as distinct (cryptic) species (Díaz-Rodríguez et al., 2017; Dufresnes et al., 2020b, 2019; Sequeira et al., 2020; van Riemsdijk et al., 2022; Wielstra and Arntzen, 2016).

Despite recent advances to identify shared biogeographic drivers at the community level (Ehl et al., 2019), there is still little consensus as to what processes contributed to the diversification of Iberian herpetofauna. Beside the effect of the Pleistocene glaciations (which were presumably more pronounced in Iberia compared to other southern European refugia; Macaluso et al., 2023), the Iberian Peninsula was affected by suites of geological events since the Middle-Miocene, including the formation of saline megalakes (the Betic salinity crisis; Lonergan and White, 1997), the accretion of the Betic mountains in southwestern Spain (Lonergan and White, 1997), as well as several marine transgression promoting insularity in southern Iberia (Achalhi et al., 2016; Martin, 2014). These events presumably promoted in situ

diversification in several Iberian species, while also mediating exchanges with northwestern Africa through land-bridges up until the Late Pliocene (Booth-Rea et al., 2018).

Until recently, mitochondrial DNA (mtDNA) has been the marker of choice in phylogeographic studies, allowing researchers to build testable hypotheses on the evolutionary history of species (Beheregaray, 2008). However, studies based entirely on mtDNA sequence data, or on only a few nuclear genetic markers, may lead to an underestimation of the true genetic diversity (e.g. Irwin et al., 2009; Milá et al., 2010), give misleading results due to introgression (e.g. Ballard and Whitlock, 2004), or lead to erroneous conclusions due to selection by maternally inherited symbionts (Hurst and Jiggins, 2005), potentially creating taxonomic chaos (Dufresnes and Jablonski, 2022). To avoid such limitations, recent phylogeographic surveys employ genome-wide, multi-locus datasets (Dufresnes et al., 2019, 2017; Gernandt et al., 2018; Huang et al., 2019; van Riemsdijk et al., 2022).

Generating genome-wide datasets can be problematic, especially when working with organisms that are characterized by massive genomes such as salamanders (Newman and Austin, 2016; Sun et al., 2012). Fortunately, whole genome data are usually not necessary, because sufficient information can be extracted from subsets of independent genetic markers (Jones and Good, 2016). Various genome reduction techniques are currently available to avoid the high-costs and computational workload of whole genome datasets (Davey et al., 2011). A particularly efficient approach is target enrichment by sequence capture, which allows for DNA regions of interest to be selected and

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Map showing the distribution and study localities of marbled newts in western Europe. The green and brown shades represent the ranges of *Triturus marmoratus* and *T. pygmaeus* (after Arntzen et al., 2007). Round symbols represent studied localities and are colored according to the genetic subgroup they belong, with northern (dark blue) and southern (light blue) T. marmoratus and western (yellow) and eastern (orange) T. pygmaeus. Localities where interspecific introgression was suggested in the Admixture analysis are marked with a white spot. Numbers correspond to individuals listed in Supplementary Table 1.</figcaption></figure>
(Zeisset and Beebee, 2008) as well as topographic barriers to dispersal. Mediterranean amphibians regularly represent groups of closely related species that are themselves characterized by high intraspecific genetic variability, with evolutionary lineages of varying ages (Ambu et al., 2023; Dufresnes et al., 2020a; Fijarczyk et al., 2011; Martinez-Solano et al., 2006; Pabijan et al., 2015; Recuero et al., 2012; Vences et al., 2013; Wielstra et al., 2013). Typically, relatively young evolutionary lineages show wide zones of genetic admixture between them, while introgression between relatively old evolutionary lineages is restricted (Arntzen et al., 2014; Dufresnes et al., 2023; Kalaentzis et al., 2023; Wielstra et al., 2018, 2017). Many of these evolutionary lineages have recently been recognized or confirmed as distinct (cryptic) species (Díaz-Rodríguez et al., 2017; Dufresnes et al., 2020b, 2019; Sequeira et al., 2020; van Riemsdijk et al., 2022; Wielstra and Arntzen, 2016).

Despite recent advances to identify shared biogeographic drivers at the community level (Ehl et al., 2019), there is still little consensus as to what processes contributed to the diversification of Iberian herpetofauna. Beside the effect of the Pleistocene glaciations (which were presumably more pronounced in Iberia compared to other southern European refugia; Macaluso et al., 2023), the Iberian Peninsula was affected by suites of geological events since the Middle-Miocene, including the formation of saline megalakes (the Betic salinity crisis; Lonergan and White, 1997), the accretion of the Betic mountains in southwestern Spain (Lonergan and White, 1997), as well as several marine transgression promoting insularity in southern Iberia (Achalhi et al., 2016; Martin, 2014). These events presumably promoted in situ

diversification in several Iberian species, while also mediating exchanges with northwestern Africa through land-bridges up until the Late Pliocene (Booth-Rea et al., 2018).

Until recently, mitochondrial DNA (mtDNA) has been the marker of choice in phylogeographic studies, allowing researchers to build testable hypotheses on the evolutionary history of species (Beheregaray, 2008). However, studies based entirely on mtDNA sequence data, or on only a few nuclear genetic markers, may lead to an underestimation of the true genetic diversity (e.g. Irwin et al., 2009; Milá et al., 2010), give misleading results due to introgression (e.g. Ballard and Whitlock, 2004), or lead to erroneous conclusions due to selection by maternally inherited symbionts (Hurst and Jiggins, 2005), potentially creating taxonomic chaos (Dufresnes and Jablonski, 2022). To avoid such limitations, recent phylogeographic surveys employ genome-wide, multi-locus datasets (Dufresnes et al., 2019, 2017; Gernandt et al., 2018; Huang et al., 2019; van Riemsdijk et al., 2022).

Generating genome-wide datasets can be problematic, especially when working with organisms that are characterized by massive genomes such as salamanders (Newman and Austin, 2016; Sun et al., 2012). Fortunately, whole genome data are usually not necessary, because sufficient information can be extracted from subsets of independent genetic markers (Jones and Good, 2016). Various genome reduction techniques are currently available to avoid the high-costs and computational workload of whole genome datasets (Davey et al., 2011). A particularly efficient approach is target enrichment by sequence capture, which allows for DNA regions of interest to be selected and

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Map showing the distribution and study localities of marbled newts in western Europe. The green and brown shades represent the ranges of *Triturus marmoratus* and *T. pygmaeus* (after Arntzen et al., 2007). Round symbols represent studied localities and are colored according to the genetic subgroup they belong, with northern (dark blue) and southern (light blue) T. marmoratus and western (yellow) and eastern (orange) T. pygmaeus. Localities where interspecific introgression was suggested in the Admixture analysis are marked with a white spot. Numbers correspond to individuals listed in Supplementary Table 1.</figcaption></figure>
sequenced by hybridization to specific probes (Andermann et al., 2020; Hutter et al., 2022; Jones and Good, 2016; McCartney-Melstad et al., 2016).

The marbled newt, *Triturus marmoratus* (Latreille, 1800), and the pygmy marbled newt, *T. pygmaeus* (Wolterstorff, 1905), are two closely related salamander species distributed in southwestern France and the Iberian Peninsula (Fig. 1). These embody the discovery of distinct evolutionary lineages that meet at a narrow hybrid zone and were eventually recognized as species, which themselves turn out to be subjected to extensive intraspecific genetic structuring. Early research on marbled newts relied on morphology, osteology, chromosome banding and allozymes (Busack et al., 1988; Dorda and Esteban, 1986; García-París et al., 1993), whereas later studies incorporated mtDNA sequence data (García-París et al., 2001). While initially recognized as subspecies, given the significant morphological and phylogenetic divergence their species status was advocated by García-París et al. (2001), which is the current consensus (Arntzen, 2023; Frost et al., 2006). The sampling scheme of 21st century studies that employed (several up to several dozens of) nuclear markers was geographically restricted, particularly to focussing on the parapatric range of the two taxa, revealing the local and restricted occurrence of hybridization (Arntzen, 2018; Arntzen et al., 2020; Espregueira Themudo and Arntzen, 2007a; López-Delgado et al., 2020). While morphological information suggests the presence of distinct intraspecific structuring (Arntzen, 2023), intraspecific phylo-geography has so far only been studied from the perspective of a single marker, namely mitochondrial DNA, and suggests shallow differentia-tion without a strong geographical component (Wielstra et al., 2013).

A range-wide, genome-wide study of marbled newts is required to determine the geographical extent of genetic admixture between the two marbled newt species, as well as the geographical genetic structuring within these species. We use target enrichment by sequence capture to collect genome-wide data for individuals sampled across the marbled newt range. This allows us to determine the extent of introgression be-tween the marbled newt species and to delineate intraspecific genetic structure. We explain our findings in the context of palaeogeological and paleoclimatic history of the Iberian Peninsula, taking into account comparative phylogeographical patterns in co-distributed Iberian amphibians.

# 2. Materials & methods

## 2.1. Sampling scheme and DNA extraction

We studied 59 marbled newt individuals, of which 28 were *T. marmoratus* and 31 were *T. pygmaeus* (Supplementary Table 1). Samples spanned the entire range of both species (Fig. 1). We included individuals from close to the species contact zone to determine the depth of introgressive hybridization (Arntzen et al., 2020, 2014; López-Delgado et al., 2020). Whole genomic DNA was extracted with the Promega Wizard™ Genomic DNA Purification Kit (Promega, Madison, WI, USA), according to the salt-based extraction protocol of Sambrook and Russell (2001).

## 2.2. Laboratory methods and data preparation

The NEBNext Ultra™ II FS DNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA) was used for library preparations, by following the manufacturer's protocol with all volumes divided by four. Additionally, the fragmentation time was optimized for *Triturus* samples (6:30 min). The isolated DNA was first sheared enzymatically, NEB adapters were then ligated and cleaved with USER enzyme (New England Biolabs, Ipswich, MA, USA). Next, size selection was performed, using NucleoMag™ magnetic beads (Macherey-Nagel, Düren, Germany) targeting an insert size of 300 bp. Unique combinations of i5 and i7 index primers (Illumina Inc., San Diego, CA, USA) were assigned to each sample via PCR amplification. Finally, the products were purified with

magnetic beads, and the concentration and fragment size of each sample was assessed with the Agilent 2200 TapeStation System (Agilent Technologies, Santa Clara, CA, USA).

Libraries were equimolarly pooled into batches of 16 samples, with a total mass of 4000 ng (250 ng per sample). For target enrichment the MyBaits-II kit (Arbor Bioscience, Ann Arbor, MI, USA; Ref# 170210-32) was used, targeting 7102 exon sequences for a total target length of c. 2.3 million bp (for details see Wielstra et al., 2019). The pooled libraries were first incubated with *Triturus*-derived C0t-1 repetitive sequence blockers for 30 min, this replaces Blocks C and O from the kit (see Wielstra et al., 2019). Then the libraries were hybridized for 30 h at 65 °C with the *Triturus*-based RNA probes designed by Wielstra et al. (2019). Following hybridization, the biotinylated RNA probes bound to the target DNA were captured with streptavidin coated magnetic beads, then washed to remove off-target DNA. The bead-bound libraries were then amplified in 14 cycles of PCR and then purified with NucleoMag™ magnetic beads. The Agilent 2200 TapeStation system was used to measure the size distribution and concentration of the pools. For each pool 16 Gb of 150 bp paired-end sequencing was performed on the Illumina NovaSeq 6000 platform (Illumina Inc., San Diego, CA, USA) by BaseClear (Leiden, the Netherlands).

## 2.3. Processing of sequence-capture data

A data-analysis pipeline was implemented to analyze the acquired sequences for each target region. First the raw reads were trimmed to remove adapters and reads of poor quality using Skewer (Jiang et al., 2014). We checked the quality of the reads before and after trimming using FastQC (Andrews, 2010). Next the reads were mapped to a reference (the *Triturus* reciprocal best blast hit reference assembly from Wielstra et al., 2019) with BWA_MEM (Li, 2013), and read groups were added and duplicates were removed with Picard (https://broadinstitute.github.io/picard/). Subsequently, Haplotype Caller and Genotype GVCFs were used in GATK (McKenna et al., 2010) to call variants, produce genomic variant call format (gVCF) files per individual, create a multi-sample VCF file by combining the individual files, and finalize this merged file by performing joint genotyping. A Hardy Weinberg Equi-librium (HWE) filter for heterozygote excess was applied to exclude putative paralogous loci from the analysis. BCFtools (Danecek et al., 2021; Li, 2011) was used to estimate the heterozygote excess p-values per sample, and sites with a score below 0.05 were removed from the dataset. Finally, additional filtering was performed on the VCF file as explained below to remove low-quality data.

The pipeline was run on two datasets. The first dataset contained 59 marbled and pygmy newt samples and was used for the Admixture, principal component analysis (PCA) and NewHybrids (see below). The second dataset was used for the phylogenetic analysis and included 52 marbled and pygmy marbled newt samples, because seven putative hybrids were removed (see Results), while three northern crested newt samples (*T. cristatus* (Laurenti,1768)) taken from Wielstra et al. (2019) were added as an outgroup, and this dataset was used in the RAxML and BEAST analyses (see below).

## 2.4. Admixture analysis

To assess the genetic differentiation within the marbled newt species group, individuals of the first dataset (without outgroup) were clustered according to their SNP genotypes using the software Admixture (Alexander et al., 2009). Input consisted of a VCF file containing one random SNP per marker and sites with no missing data, resulting in 5,951 SNPs. For the number of subpopulations (K), values from 1 to 20 were tested and the analysis was run for 25 replicates. Replicates for each K value were combined with CLUMPAK (Kopelman et al., 2015). Afterwards, the optimal K value was determined by Admixture's cross validation error rate. The results were visualized using the R packages reshape2 (Wickham, 2007) and ggplot2 (Wickham, 2016).
sequenced by hybridization to specific probes (Andermann et al., 2020; Hutter et al., 2022; Jones and Good, 2016; McCartney-Melstad et al., 2016).

The marbled newt, *Triturus marmoratus* (Latreille, 1800), and the pygmy marbled newt, *T. pygmaeus* (Wolterstorff, 1905), are two closely related salamander species distributed in southwestern France and the Iberian Peninsula (Fig. 1). These embody the discovery of distinct evolutionary lineages that meet at a narrow hybrid zone and were eventually recognized as species, which themselves turn out to be subjected to extensive intraspecific genetic structuring. Early research on marbled newts relied on morphology, osteology, chromosome banding and allozymes (Busack et al., 1988; Dorda and Esteban, 1986; García-París et al., 1993), whereas later studies incorporated mtDNA sequence data (García-París et al., 2001). While initially recognized as subspecies, given the significant morphological and phylogenetic divergence their species status was advocated by García-París et al. (2001), which is the current consensus (Arntzen, 2023; Frost et al., 2006). The sampling scheme of 21st century studies that employed (several up to several dozens of) nuclear markers was geographically restricted, particularly to focussing on the parapatric range of the two taxa, revealing the local and restricted occurrence of hybridization (Arntzen, 2018; Arntzen et al., 2020; Espregueira Themudo and Arntzen, 2007a; López-Delgado et al., 2020). While morphological information suggests the presence of distinct intraspecific structuring (Arntzen, 2023), intraspecific phylo-geography has so far only been studied from the perspective of a single marker, namely mitochondrial DNA, and suggests shallow differentia-tion without a strong geographical component (Wielstra et al., 2013).

A range-wide, genome-wide study of marbled newts is required to determine the geographical extent of genetic admixture between the two marbled newt species, as well as the geographical genetic structuring within these species. We use target enrichment by sequence capture to collect genome-wide data for individuals sampled across the marbled newt range. This allows us to determine the extent of introgression be-tween the marbled newt species and to delineate intraspecific genetic structure. We explain our findings in the context of palaeogeological and paleoclimatic history of the Iberian Peninsula, taking into account comparative phylogeographical patterns in co-distributed Iberian amphibians.

# 2. Materials & methods

## 2.1. Sampling scheme and DNA extraction

We studied 59 marbled newt individuals, of which 28 were *T. marmoratus* and 31 were *T. pygmaeus* (Supplementary Table 1). Samples spanned the entire range of both species (Fig. 1). We included individuals from close to the species contact zone to determine the depth of introgressive hybridization (Arntzen et al., 2020, 2014; López-Delgado et al., 2020). Whole genomic DNA was extracted with the Promega Wizard™ Genomic DNA Purification Kit (Promega, Madison, WI, USA), according to the salt-based extraction protocol of Sambrook and Russell (2001).

## 2.2. Laboratory methods and data preparation

The NEBNext Ultra™ II FS DNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA) was used for library preparations, by following the manufacturer's protocol with all volumes divided by four. Additionally, the fragmentation time was optimized for *Triturus* samples (6:30 min). The isolated DNA was first sheared enzymatically, NEB adapters were then ligated and cleaved with USER enzyme (New England Biolabs, Ipswich, MA, USA). Next, size selection was performed, using NucleoMag™ magnetic beads (Macherey-Nagel, Düren, Germany) targeting an insert size of 300 bp. Unique combinations of i5 and i7 index primers (Illumina Inc., San Diego, CA, USA) were assigned to each sample via PCR amplification. Finally, the products were purified with

magnetic beads, and the concentration and fragment size of each sample was assessed with the Agilent 2200 TapeStation System (Agilent Technologies, Santa Clara, CA, USA).

Libraries were equimolarly pooled into batches of 16 samples, with a total mass of 4000 ng (250 ng per sample). For target enrichment the MyBaits-II kit (Arbor Bioscience, Ann Arbor, MI, USA; Ref# 170210-32) was used, targeting 7102 exon sequences for a total target length of c. 2.3 million bp (for details see Wielstra et al., 2019). The pooled libraries were first incubated with *Triturus*-derived C0t-1 repetitive sequence blockers for 30 min, this replaces Blocks C and O from the kit (see Wielstra et al., 2019). Then the libraries were hybridized for 30 h at 65 °C with the *Triturus*-based RNA probes designed by Wielstra et al. (2019). Following hybridization, the biotinylated RNA probes bound to the target DNA were captured with streptavidin coated magnetic beads, then washed to remove off-target DNA. The bead-bound libraries were then amplified in 14 cycles of PCR and then purified with NucleoMag™ magnetic beads. The Agilent 2200 TapeStation system was used to measure the size distribution and concentration of the pools. For each pool 16 Gb of 150 bp paired-end sequencing was performed on the Illumina NovaSeq 6000 platform (Illumina Inc., San Diego, CA, USA) by BaseClear (Leiden, the Netherlands).

## 2.3. Processing of sequence-capture data

A data-analysis pipeline was implemented to analyze the acquired sequences for each target region. First the raw reads were trimmed to remove adapters and reads of poor quality using Skewer (Jiang et al., 2014). We checked the quality of the reads before and after trimming using FastQC (Andrews, 2010). Next the reads were mapped to a reference (the *Triturus* reciprocal best blast hit reference assembly from Wielstra et al., 2019) with BWA_MEM (Li, 2013), and read groups were added and duplicates were removed with Picard (https://broadinstitute.github.io/picard/). Subsequently, Haplotype Caller and Genotype GVCFs were used in GATK (McKenna et al., 2010) to call variants, produce genomic variant call format (gVCF) files per individual, create a multi-sample VCF file by combining the individual files, and finalize this merged file by performing joint genotyping. A Hardy Weinberg Equi-librium (HWE) filter for heterozygote excess was applied to exclude putative paralogous loci from the analysis. BCFtools (Danecek et al., 2021; Li, 2011) was used to estimate the heterozygote excess p-values per sample, and sites with a score below 0.05 were removed from the dataset. Finally, additional filtering was performed on the VCF file as explained below to remove low-quality data.

The pipeline was run on two datasets. The first dataset contained 59 marbled and pygmy newt samples and was used for the Admixture, principal component analysis (PCA) and NewHybrids (see below). The second dataset was used for the phylogenetic analysis and included 52 marbled and pygmy marbled newt samples, because seven putative hybrids were removed (see Results), while three northern crested newt samples (*T. cristatus* (Laurenti,1768)) taken from Wielstra et al. (2019) were added as an outgroup, and this dataset was used in the RAxML and BEAST analyses (see below).

## 2.4. Admixture analysis

To assess the genetic differentiation within the marbled newt species group, individuals of the first dataset (without outgroup) were clustered according to their SNP genotypes using the software Admixture (Alexander et al., 2009). Input consisted of a VCF file containing one random SNP per marker and sites with no missing data, resulting in 5,951 SNPs. For the number of subpopulations (K), values from 1 to 20 were tested and the analysis was run for 25 replicates. Replicates for each K value were combined with CLUMPAK (Kopelman et al., 2015). Afterwards, the optimal K value was determined by Admixture's cross validation error rate. The results were visualized using the R packages reshape2 (Wickham, 2007) and ggplot2 (Wickham, 2016).
## 2.5. Principal component analysis

A principal component analysis was carried out using the sequence alignment VCF file from the first dataset (without outgroup) to estimate the genetic population structure of the marbled newt species group. The snpgdsPCA function of the SNPRelate R package (Zheng et al., 2012) was used under the default parameters to conduct PCA. We used SNPs with a minimum quality score of 20 and no more than 50 % missing data, while additional pruning was performed by the R script to remove monomorphic sites, after which 105,753 SNPs were kept.

## 2.6. NewHybrids analysis

NewHybrids (Anderson and Thompson, 2002) is a Bayesian model-based clustering method that computes the posterior probability for each individual to belong to predefined, distinct parental and hybrid classes over a genealogical depth set by the user. For seven presumed hybrid individuals identified in the Admixture analysis (see Results, marked with a white spot in Fig. 1), we investigated the depth of hybridization with the program NewHybrids. The one SNP per target VCF file generated for Admixture was converted into NewHybrids format using PGD Spider v.2.1.1.5 (Lischer and Excoffier, 2012). All non-hybrids were marked as either *T. marmoratus* or *T. pygmaeus* by invoking the z-option. Under the 'two generations' boundary condition, individuals are allocated to up to six classes (two parental, F1 and F2 hybrids, and backcross hybrids in either direction). Under the 'three generations' setting, another nine discernable hybrid categories are considered. The program was run for 10,000 iterations under default conditions. Because of capacity limitations, the dataset was reduced to the 345 most informative loci (the largest number with which we could run NewHybrids). The criterion used was Cohen's kappa (with 1 - kappa acting as a measure of disconcordance, i.e. diagnosticity).

## 2.7. Concatenated maximum likelihood phylogeny

The second dataset (excluding admixed marbled newts and including the crested newt outgroup) was subjected to a concatenated maximum likelihood analysis with RAxML (Stamatakis, 2014). We used SNPs with a minimum quality score of 20 and no more than 50 % missing data. The filtered VCF file from the data-analysis pipeline was converted into PHYLIP format using PGD Spider v.2.1.1.5 (Lischer and Excoffier, 2012). Next, the 'ascbias.py' Python script (https://github.com/btmartin721/raxml_ascbias) was used to exclude invariant sites from the analysis. The phylogenetic tree was built based on an alignment of 122,989 SNPs from 6,884 target loci. The GTR + GAMMA substitution model was used for the analysis, with 100 bootstrap replicates for branch-support estimation, followed by a correction for ascertainment bias (Lewis' ascertainment correction; Lewis, 2001). Finally, the resulting phylogenetic tree was visualized with FigTree 1.4.3 (https://www.tree.bio.ed.ac.uk/software/figtree).

## 2.8. Molecular dating

In order to obtain a dated phylogeny, we conducted a concatenated Bayesian analysis in BEAST 2.7 (Bouckaert et al., 2014). This was based on SNPs with a minimum quality score of 20 and no more than 50 % missing data. We initially converted the filtered VCF file into NEXUS format with Mesquite 3.81 (Maddison and Maddison, 2023) and next to XML in BEAUti 2.7 (part of the BEAST package). A fossil dated at 24 Mya was used as a minimum estimate for the most recent common ancestor of the genus *Triturus* (Steinfartz et al., 2007) and given a lognormally distributed prior with an offset of 24 and a mean and standard deviation of 1.0. We applied the GTR + I + G model of sequence evolution, the strict clock model and the Yule speciation model. We conducted a 25 million generation run with a sampling frequency of 0.0001 on the Cipres Science Gateway (Miller et al., 2010). We confirmed that sample

sizes for parameters were > 100 (and in most cases ≫ 200) in TRACER v.1.7 (Rambaut et al., 2018), discarded the first half of the sampled trees as burn-in and calculated an annotated Maximum Clade Credibility tree in TreeAnnotator 2.7 (part of the BEAST package).

## 2.9. mtDNA analysis

We sequenced 658 bp of an mtDNA gene (ND4) with primer pair KARF4-KAR1 following Wielstra et al. (2013). Sanger sequencing was outsourced to BaseClear B.V. and sequences were edited and trimmed in Geneious Prime 2022.1.1. We conducted calibrated Bayesian inference in BEAST v.2.1, using the same settings as above, with the mtDNA partitioned by codon position.

## 3. Results

A mean of 10,958,752 read pairs per sample was generated, with a standard deviation of 3,383,913. An average of 30.03 % raw reads (s.d. 11.98 %) were mapped. The percentage of duplicate reads, which were filtered out, accounted for an average of 44.20 % of mapped reads (s.d. 11.28 %).

In the Admixture analysis, K = 4 yields the best fit for the entire marbled newt dataset (Fig. 2). These four groups are geographically coherent, with two *T. marmoratus* groups in the north versus the south of the species' range, and two *T. pygmaeus* groups in the west versus the east of the species' range (Fig. 1). An Admixture analysis for K = 2 reveals limited genetic admixture between the two species based on individual ancestry fractions ($Q_m$). Three *T. marmoratus* individuals were classified as admixed (0.92 > $Q_m$ > 0.93), as well as four *T. pygmaeus* individuals (0.02 < $Q_m$ < 0.17) (Fig. 2).

In the principal component analysis, the first axis (PC1) accounts for 18.8 % of the total variance and separates *T. marmoratus* from *T. pygmaeus* (Fig. 3). The second principal component (PC2) accounts for 3.27 % of the total variance and splits each species into two genetic groups. The results are overall similar to those of the Admixture analysis, with a northern and a southern group in *T. marmoratus* and a western and an eastern group in *T. pygmaeus*, but note that no intermediates between *T. marmoratus* – *T. pygmaeus* are found on PC1, and a single intermediate is found between the *T. marmoratus* subgroups on PC2 (individual P1093).

Under the 'two generations' constraint in NewHybrids, six out of seven of the presumed hybrids are allocated to one or the other parental species, whereas individual P7803 is classified as a backcross hybrid in direction of *T. pygmaeus*, all with probabilities > 0.99 (details are in Supplementary Table 2). Conversely, under the 'three generations' boundary condition, all individuals are classified as third generation hybrids, with the exception of individuals P7814 and P7824, which remain classified as *T. pygmaeus*.

The concatenated phylogenetic analysis with RAxML, in which the seven genetically admixed individuals were removed, resolves both marbled newt species as monophyletic with high bootstrap support (Fig. 4). Each species is divided into two major clades with high branch support, in general agreement with the Admixture and PCA analyses. The BEAST analysis dates the split between *T. marmoratus* and *T. pygmaeus* to the late Miocene at c. 5 Ma and the basal splits in both *T. marmoratus* and *T. pygmaeus* to the early Pleistocene at c. 2 Ma (Fig. 5a). Note that the position of individual 1093, identified as inter-mediary between the two groups within *T. marmoratus*, deviates between the RAxML and the BEAST analysis.

We obtained mtDNA sequences for all but two individuals. We do not observe interspecific introgression of mtDNA, i.e. we do not see individuals that are predominantly *T. marmoratus* based on nuclear DNA with *T. pygmaeus* mtDNA or vice versa. All mtDNA haplotypes present in marbled newts identified as genetically admixed based on nuclear DNA are also present in individuals that are identified as genetically pure. The two species are represented by distinct mtDNA clades for which the split
## 2.5. Principal component analysis

A principal component analysis was carried out using the sequence alignment VCF file from the first dataset (without outgroup) to estimate the genetic population structure of the marbled newt species group. The snpgdsPCA function of the SNPRelate R package (Zheng et al., 2012) was used under the default parameters to conduct PCA. We used SNPs with a minimum quality score of 20 and no more than 50 % missing data, while additional pruning was performed by the R script to remove monomorphic sites, after which 105,753 SNPs were kept.

## 2.6. NewHybrids analysis

NewHybrids (Anderson and Thompson, 2002) is a Bayesian model-based clustering method that computes the posterior probability for each individual to belong to predefined, distinct parental and hybrid classes over a genealogical depth set by the user. For seven presumed hybrid individuals identified in the Admixture analysis (see Results, marked with a white spot in Fig. 1), we investigated the depth of hybridization with the program NewHybrids. The one SNP per target VCF file generated for Admixture was converted into NewHybrids format using PGD Spider v.2.1.1.5 (Lischer and Excoffier, 2012). All non-hybrids were marked as either *T. marmoratus* or *T. pygmaeus* by invoking the z-option. Under the 'two generations' boundary condition, individuals are allocated to up to six classes (two parental, F1 and F2 hybrids, and backcross hybrids in either direction). Under the 'three generations' setting, another nine discernable hybrid categories are considered. The program was run for 10,000 iterations under default conditions. Because of capacity limitations, the dataset was reduced to the 345 most informative loci (the largest number with which we could run NewHybrids). The criterion used was Cohen's kappa (with 1 - kappa acting as a measure of disconcordance, i.e. diagnosticity).

## 2.7. Concatenated maximum likelihood phylogeny

The second dataset (excluding admixed marbled newts and including the crested newt outgroup) was subjected to a concatenated maximum likelihood analysis with RAxML (Stamatakis, 2014). We used SNPs with a minimum quality score of 20 and no more than 50 % missing data. The filtered VCF file from the data-analysis pipeline was converted into PHYLIP format using PGD Spider v.2.1.1.5 (Lischer and Excoffier, 2012). Next, the 'ascbias.py' Python script (https://github.com/btmartin721/raxml_ascbias) was used to exclude invariant sites from the analysis. The phylogenetic tree was built based on an alignment of 122,989 SNPs from 6,884 target loci. The GTR + GAMMA substitution model was used for the analysis, with 100 bootstrap replicates for branch-support estimation, followed by a correction for ascertainment bias (Lewis' ascertainment correction; Lewis, 2001). Finally, the resulting phylogenetic tree was visualized with FigTree 1.4.3 (https://www.tree.bio.ed.ac.uk/software/figtree).

## 2.8. Molecular dating

In order to obtain a dated phylogeny, we conducted a concatenated Bayesian analysis in BEAST 2.7 (Bouckaert et al., 2014). This was based on SNPs with a minimum quality score of 20 and no more than 50 % missing data. We initially converted the filtered VCF file into NEXUS format with Mesquite 3.81 (Maddison and Maddison, 2023) and next to XML in BEAUti 2.7 (part of the BEAST package). A fossil dated at 24 Mya was used as a minimum estimate for the most recent common ancestor of the genus *Triturus* (Steinfartz et al., 2007) and given a lognormally distributed prior with an offset of 24 and a mean and standard deviation of 1.0. We applied the GTR + I + G model of sequence evolution, the strict clock model and the Yule speciation model. We conducted a 25 million generation run with a sampling frequency of 0.0001 on the Cipres Science Gateway (Miller et al., 2010). We confirmed that sample

sizes for parameters were > 100 (and in most cases ≫ 200) in TRACER v.1.7 (Rambaut et al., 2018), discarded the first half of the sampled trees as burn-in and calculated an annotated Maximum Clade Credibility tree in TreeAnnotator 2.7 (part of the BEAST package).

## 2.9. mtDNA analysis

We sequenced 658 bp of an mtDNA gene (ND4) with primer pair KARF4-KAR1 following Wielstra et al. (2013). Sanger sequencing was outsourced to BaseClear B.V. and sequences were edited and trimmed in Geneious Prime 2022.1.1. We conducted calibrated Bayesian inference in BEAST v.2.1, using the same settings as above, with the mtDNA partitioned by codon position.

## 3. Results

A mean of 10,958,752 read pairs per sample was generated, with a standard deviation of 3,383,913. An average of 30.03 % raw reads (s.d. 11.98 %) were mapped. The percentage of duplicate reads, which were filtered out, accounted for an average of 44.20 % of mapped reads (s.d. 11.28 %).

In the Admixture analysis, K = 4 yields the best fit for the entire marbled newt dataset (Fig. 2). These four groups are geographically coherent, with two *T. marmoratus* groups in the north versus the south of the species' range, and two *T. pygmaeus* groups in the west versus the east of the species' range (Fig. 1). An Admixture analysis for K = 2 reveals limited genetic admixture between the two species based on individual ancestry fractions ($Q_m$). Three *T. marmoratus* individuals were classified as admixed (0.92 > $Q_m$ > 0.93), as well as four *T. pygmaeus* individuals (0.02 < $Q_m$ < 0.17) (Fig. 2).

In the principal component analysis, the first axis (PC1) accounts for 18.8 % of the total variance and separates *T. marmoratus* from *T. pygmaeus* (Fig. 3). The second principal component (PC2) accounts for 3.27 % of the total variance and splits each species into two genetic groups. The results are overall similar to those of the Admixture analysis, with a northern and a southern group in *T. marmoratus* and a western and an eastern group in *T. pygmaeus*, but note that no intermediates between *T. marmoratus* – *T. pygmaeus* are found on PC1, and a single intermediate is found between the *T. marmoratus* subgroups on PC2 (individual P1093).

Under the 'two generations' constraint in NewHybrids, six out of seven of the presumed hybrids are allocated to one or the other parental species, whereas individual P7803 is classified as a backcross hybrid in direction of *T. pygmaeus*, all with probabilities > 0.99 (details are in Supplementary Table 2). Conversely, under the 'three generations' boundary condition, all individuals are classified as third generation hybrids, with the exception of individuals P7814 and P7824, which remain classified as *T. pygmaeus*.

The concatenated phylogenetic analysis with RAxML, in which the seven genetically admixed individuals were removed, resolves both marbled newt species as monophyletic with high bootstrap support (Fig. 4). Each species is divided into two major clades with high branch support, in general agreement with the Admixture and PCA analyses. The BEAST analysis dates the split between *T. marmoratus* and *T. pygmaeus* to the late Miocene at c. 5 Ma and the basal splits in both *T. marmoratus* and *T. pygmaeus* to the early Pleistocene at c. 2 Ma (Fig. 5a). Note that the position of individual 1093, identified as inter-mediary between the two groups within *T. marmoratus*, deviates between the RAxML and the BEAST analysis.

We obtained mtDNA sequences for all but two individuals. We do not observe interspecific introgression of mtDNA, i.e. we do not see individuals that are predominantly *T. marmoratus* based on nuclear DNA with *T. pygmaeus* mtDNA or vice versa. All mtDNA haplotypes present in marbled newts identified as genetically admixed based on nuclear DNA are also present in individuals that are identified as genetically pure. The two species are represented by distinct mtDNA clades for which the split
<figure><img src="image_3.png" /><figcaption>Fig. 2. Individual genetic admixture proportions bar plots for marbled newts estimated with Admixture based on 5,951 nuclear DNA markers (a single SNP per marker). The best-fitting model had a cross-validation error of 0.1 and four ancestry components (K = 4, bottom panel). The model with two ancestry components (K = 2, top panel) was created to estimate admixture between the species with *Triturus marmoratus* in blue and *T. pygmaeus* in orange, with seven admixed individuals shown by an asterisk. Sample numbers correspond to Supplementary Table 1.</figcaption></figure>

is again dated to c. 5 Ma (Fig. 5b). The intraspecific groups determined by the nuclear DNA analysis are not recovered in the mtDNA BEAST analysis.

# 4. Discussion

## 4.1. Near-absence of gene flow between marbled newt species

The separation of *T. marmoratus* and *T. pygmaeus* in all the analyses supports their taxonomic distinction, and the little admixture between them is consistent with the species level (Espregueira Themudo and Arntzen, 2007b; Wielstra et al., 2019). Admixture is geographically restricted: the seven admixed individuals all originate from near the parapatric range border of the two species (Figures 1 & 2), where introgressive hybridization has been previously documented (Arntzen, 2018; Arntzen et al., 2020). All but two of the admixed individuals are identified as deeper generation backcrosses. One of the admixed individuals (individual P2888) represents a *T. marmoratus* enclave that is enveloped by *T. pygmaeus* (Espregueira Themudo and Arntzen, 2007a), resulting from a near-complete species replacement (Arntzen et al., 2020; López-Delgado et al., 2020). This population managing to maintain its genetic integrity is an indication that introgression is selected against, in spite of the apparent opportunity for hybridization and gene flow.

Molecular dating based on either nuclear DNA or mtDNA suggests a split between the two marbled newt species around 5 Ma (Fig. 5). This dating roughly coincides with the diversification of other Iberian amphibians, including the *Pelodytes* genus (Díaz-Rodríguez et al., 2017; Dufresnes et al., 2020c), the *Rana temporaria* complex (Dufresnes et al., 2020b), or the emergence of *Baleaphryne* clade of *Alytes* (Ambu et al., 2023). However, it is, with respect to the distributions of these species groups, hard to pinpoint a single paleographic event responsible for their simultaneous divergences. Here, we see two plausible scenarios for the split of *T. marmoratus* and *T. pygmaeus*. On the one hand, it could be

the consequence of the Atlantic-Mediterranean seaways that momentarily re-opened the Betic Straits in the Betic Cordillera by the end of the Miocene (Martin, 2014), subsequently isolating Southern Iberia from the rest of the peninsula. On the other hand, the ancestor of both our *Triturus* species may have historically colonized North Africa during pre-Messinian times, e.g., by Late Miocene land-bridges across the Alboran volcanic-arc (which opened as early as 7Mya; Booth-Rea et al., 2018). The ancestor would have diverged following the refilling of the Mediterranean Sea as the Strait of Gibraltar re-opened, allowing the neoformed *T. pygmaeus* to make it back to Iberia (and subsequently disappearing from North Africa). Post-Messinian faunal dispersal events are also known from another newt, *Pleurodeles waltl* (Gutiérrez-Rodríguez et al., 2017a). In particular, this scenario was proposed to explain the unexpected phylogenetic relationship between the Betic and Moroccan endemic midwife toad species *Alytes dickhilleni* and *A. maurus* (Ambu et al., 2023). Whatever the initial drivers of differentiation, the subsequent divergence of *T. pygmaeus* is likely to have been promoted by the isolation of Southern Iberia throughout much of the Pliocene and Pleistocene, as environmental connections with west and east Iberia became modulated by climatic fluctuations. Consequently, the region is inhabited by many regional endemic amphibians (e.g., *Pelodytes ibericus*, *Alytes dickhilleni*, *Discoglossus galganoi* jeanneae, *Salamandra salamandra longirostris*).

## 4.2. Genetic structure within *T. marmoratus*

Our analyses reveal two subgroups within *T. marmoratus* with a split dated around 2 Ma, during the Quaternary Ice Age, and show that most genetic diversity of the species is in the Iberian Peninsula. The northern subgroup covers a wide area, ranging from the west of France to the north of Portugal. The southern subgroup has a narrow distribution across north and central Portugal and adjacent Spain. Evidence of genetic admixture between the subgroups is found in the north of Portugal (individual P1093) and north of the Central System (Figs. 2 and 3). On
<figure><img src="image_4.png" /><figcaption>Fig. 2. Individual genetic admixture proportions bar plots for marbled newts estimated with Admixture based on 5,951 nuclear DNA markers (a single SNP per marker). The best-fitting model had a cross-validation error of 0.1 and four ancestry components (K = 4, bottom panel). The model with two ancestry components (K = 2, top panel) was created to estimate admixture between the species with *Triturus marmoratus* in blue and *T. pygmaeus* in orange, with seven admixed individuals shown by an asterisk. Sample numbers correspond to Supplementary Table 1.</figcaption></figure>

is again dated to c. 5 Ma (Fig. 5b). The intraspecific groups determined by the nuclear DNA analysis are not recovered in the mtDNA BEAST analysis.

# 4. Discussion

## 4.1. Near-absence of gene flow between marbled newt species

The separation of *T. marmoratus* and *T. pygmaeus* in all the analyses supports their taxonomic distinction, and the little admixture between them is consistent with the species level (Espregueira Themudo and Arntzen, 2007b; Wielstra et al., 2019). Admixture is geographically restricted: the seven admixed individuals all originate from near the parapatric range border of the two species (Figures 1 & 2), where introgressive hybridization has been previously documented (Arntzen, 2018; Arntzen et al., 2020). All but two of the admixed individuals are identified as deeper generation backcrosses. One of the admixed individuals (individual P2888) represents a *T. marmoratus* enclave that is enveloped by *T. pygmaeus* (Espregueira Themudo and Arntzen, 2007a), resulting from a near-complete species replacement (Arntzen et al., 2020; López-Delgado et al., 2020). This population managing to maintain its genetic integrity is an indication that introgression is selected against, in spite of the apparent opportunity for hybridization and gene flow.

Molecular dating based on either nuclear DNA or mtDNA suggests a split between the two marbled newt species around 5 Ma (Fig. 5). This dating roughly coincides with the diversification of other Iberian amphibians, including the *Pelodytes* genus (Díaz-Rodríguez et al., 2017; Dufresnes et al., 2020c), the *Rana temporaria* complex (Dufresnes et al., 2020b), or the emergence of *Baleaphryne* clade of *Alytes* (Ambu et al., 2023). However, it is, with respect to the distributions of these species groups, hard to pinpoint a single paleographic event responsible for their simultaneous divergences. Here, we see two plausible scenarios for the split of *T. marmoratus* and *T. pygmaeus*. On the one hand, it could be

the consequence of the Atlantic-Mediterranean seaways that momentarily re-opened the Betic Straits in the Betic Cordillera by the end of the Miocene (Martin, 2014), subsequently isolating Southern Iberia from the rest of the peninsula. On the other hand, the ancestor of both our *Triturus* species may have historically colonized North Africa during pre-Messinian times, e.g., by Late Miocene land-bridges across the Alboran volcanic-arc (which opened as early as 7Mya; Booth-Rea et al., 2018). The ancestor would have diverged following the refilling of the Mediterranean Sea as the Strait of Gibraltar re-opened, allowing the neoformed *T. pygmaeus* to make it back to Iberia (and subsequently disappearing from North Africa). Post-Messinian faunal dispersal events are also known from another newt, *Pleurodeles waltl* (Gutiérrez-Rodríguez et al., 2017a). In particular, this scenario was proposed to explain the unexpected phylogenetic relationship between the Betic and Moroccan endemic midwife toad species *Alytes dickhilleni* and *A. maurus* (Ambu et al., 2023). Whatever the initial drivers of differentiation, the subsequent divergence of *T. pygmaeus* is likely to have been promoted by the isolation of Southern Iberia throughout much of the Pliocene and Pleistocene, as environmental connections with west and east Iberia became modulated by climatic fluctuations. Consequently, the region is inhabited by many regional endemic amphibians (e.g., *Pelodytes ibericus*, *Alytes dickhilleni*, *Discoglossus galganoi* jeanneae, *Salamandra salamandra longirostris*).

## 4.2. Genetic structure within *T. marmoratus*

Our analyses reveal two subgroups within *T. marmoratus* with a split dated around 2 Ma, during the Quaternary Ice Age, and show that most genetic diversity of the species is in the Iberian Peninsula. The northern subgroup covers a wide area, ranging from the west of France to the north of Portugal. The southern subgroup has a narrow distribution across north and central Portugal and adjacent Spain. Evidence of genetic admixture between the subgroups is found in the north of Portugal (individual P1093) and north of the Central System (Figs. 2 and 3). On
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Principal component analysis of the marbled newts *Triturus marmoratus* (shown in light and dark blue) and *T. pygmaeus* (yellow and orange). Admixed samples are highlighted with a white spot (see Fig. 2). Sample numbers correspond to Supplementary Table 1.</figcaption></figure>

balance, these observations suggest that the genetic transition of the subgroups is likely to be fairly wide.

Divergence of populations from central and northern Iberia may be the consequence of major rivers flowing in a West-East orientation acting as a barrier to dispersal, initiating incipient speciation. Similar patterns of differentiation are observed in other amphibians, such as *Chioglossa lusitanica* (Sequeira et al., 2022), *Lissotriton boscai* (Martinez-Solano et al., 2006), *Alytes obstetricans* (Dufresnes and Hernandez, 2021; Gonçalves et al., 2015; Maia-Carvalho et al., 2018), and *Discoglossus galganoi* (Martinez-Solano, 2004). Remarkably, the phylogeographic transition between the two *T. marmoratus* lineages coincides with the contact area of two phylogeographic lineages attributed to the midwife toad *A. o. boscai* (Maia-Carvalho et al., 2018). Amphibian dispersal in this area may be restricted by the Douro River, the highest flow river of the Iberian Peninsula, which forms a > 250 m wide brackish estuary tens of kilometers from its mouth.

*Triturus marmoratus* samples from northern Spain show genetic similarity with populations in the west of France. This observation suggests either a large glacial refugium encompassing northern Spain and France, or the post-glacial expansion from northern Spain into France. A Spanish origin and an Atlantic corridor of colonization of western Europe has been suggested for many amphibians (e.g. *Hyla molleri*, Sánchez-Montes et al., 2019; and *Lissotriton helveticus*, Recuero and García-París, 2011). However, it is unclear whether the *Triturus*

refugium also encompassed the French side of the Pyrenean Mountains, as could be the case in *Alytes o. obstetricans* (Ambu et al., 2023), and *Bufo spinosus* (Arntzen et al., 2017), and as suggested for other species (e.g., *Epidalea calamita*, Rowe et al., 2006). Moreover, a Mediterranean corridor cannot be excluded. For instance, the widespread *Pelobates cultripes* colonized France from the Mediterranean coast, and then expanded westward inland to reach the Atlantic coast through the Aude and Garonne valleys, which offered suitable conditions during the Holocene (Bailon, 2003; Gutiérrez-Rodríguez et al., 2017b). Likewise, *Pelodytes p. punctatus* expanded from Catalonia to most of Western France, where it often co-occurs with *T. marmoratus* (Díaz-Rodríguez et al., 2015; Dufresnes, 2019). Either way, the lack of genetic structure and diversity between France and Spain (Fig. 4), and the predicted absence of the species in model-based glacial distributions (Wielstra et al., 2013), clearly support out-of-Iberia post-glacial colonization.

### 4.3. Genetic structure within *T. pygmaeus*

Our analyses also revealed the presence of two subgroups within *T. pygmaeus*, again with a Quaternary Ice Age split dated around 2 Ma, with a western one restricted to the western part of central Portugal and an eastern one covering the remainder of the species' range. Intraspecific introgression was observed in the Admixture analysis not only in samples from near the contact zone of the two subgroups, but all the way
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Principal component analysis of the marbled newts *Triturus marmoratus* (shown in light and dark blue) and *T. pygmaeus* (yellow and orange). Admixed samples are highlighted with a white spot (see Fig. 2). Sample numbers correspond to Supplementary Table 1.</figcaption></figure>

balance, these observations suggest that the genetic transition of the subgroups is likely to be fairly wide.

Divergence of populations from central and northern Iberia may be the consequence of major rivers flowing in a West-East orientation acting as a barrier to dispersal, initiating incipient speciation. Similar patterns of differentiation are observed in other amphibians, such as *Chioglossa lusitanica* (Sequeira et al., 2022), *Lissotriton boscai* (Martinez-Solano et al., 2006), *Alytes obstetricans* (Dufresnes and Hernandez, 2021; Gonçalves et al., 2015; Maia-Carvalho et al., 2018), and *Discoglossus galganoi* (Martinez-Solano, 2004). Remarkably, the phylogeographic transition between the two *T. marmoratus* lineages coincides with the contact area of two phylogeographic lineages attributed to the midwife toad *A. o. boscai* (Maia-Carvalho et al., 2018). Amphibian dispersal in this area may be restricted by the Douro River, the highest flow river of the Iberian Peninsula, which forms a > 250 m wide brackish estuary tens of kilometers from its mouth.

*Triturus marmoratus* samples from northern Spain show genetic similarity with populations in the west of France. This observation suggests either a large glacial refugium encompassing northern Spain and France, or the post-glacial expansion from northern Spain into France. A Spanish origin and an Atlantic corridor of colonization of western Europe has been suggested for many amphibians (e.g. *Hyla molleri*, Sánchez-Montes et al., 2019; and *Lissotriton helveticus*, Recuero and García-París, 2011). However, it is unclear whether the *Triturus*

refugium also encompassed the French side of the Pyrenean Mountains, as could be the case in *Alytes o. obstetricans* (Ambu et al., 2023), and *Bufo spinosus* (Arntzen et al., 2017), and as suggested for other species (e.g., *Epidalea calamita*, Rowe et al., 2006). Moreover, a Mediterranean corridor cannot be excluded. For instance, the widespread *Pelobates cultripes* colonized France from the Mediterranean coast, and then expanded westward inland to reach the Atlantic coast through the Aude and Garonne valleys, which offered suitable conditions during the Holocene (Bailon, 2003; Gutiérrez-Rodríguez et al., 2017b). Likewise, *Pelodytes p. punctatus* expanded from Catalonia to most of Western France, where it often co-occurs with *T. marmoratus* (Díaz-Rodríguez et al., 2015; Dufresnes, 2019). Either way, the lack of genetic structure and diversity between France and Spain (Fig. 4), and the predicted absence of the species in model-based glacial distributions (Wielstra et al., 2013), clearly support out-of-Iberia post-glacial colonization.

### 4.3. Genetic structure within *T. pygmaeus*

Our analyses also revealed the presence of two subgroups within *T. pygmaeus*, again with a Quaternary Ice Age split dated around 2 Ma, with a western one restricted to the western part of central Portugal and an eastern one covering the remainder of the species' range. Intraspecific introgression was observed in the Admixture analysis not only in samples from near the contact zone of the two subgroups, but all the way
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Phylogenetic tree for marbled newts built with RAxML based on an alignment of 122,989 SNPs representing 6,884 nuclear DNA markers. *Triturus marmoratus* individuals are shown in light and dark blue and *T. pygmaeus* individuals in yellow and orange. Outgroup individuals of *T. cristatus* are shown in grey. Some long branches were truncated, as shown by double slashes (//). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

from central to southern Portugal (Fig. 2). A similar pattern of genetic differentiation is observed in other amphibian species in this region (e.g. *Pelodytes atlanticus*, Díaz-Rodríguez et al., 2017; and *Lissotriton boscai*, Peñalver-Alcázar et al., 2021; Sequeira et al., 2020; Teixeira et al., 2015). The limited data available so far suggest that the Tagus river may be separating the two groups. Alternatively, regional uplifts linked to Pleistocene tectonic activities in southwestern Iberia (Figueiredo et al., 2014), combined with Quaternary sea level changes, may have locally promoted population isolation.

There is little differentiation in the eastern part of the range, pre-
sumably due to postglacial expansion from a southwestern refugium to
the central parts of the Iberian Peninsula (Arntzen, 2018; Arntzen and
Espregueira Themudo, 2008; Espregueira Themudo et al., 2012; Espre-
gueira Themudo and Arntzen, 2007a; Wielstra et al., 2013). Such phy-
logeographic homogeneity is also tangible from other amphibians,
including *Alytes cisternasii* (Ambu et al., 2023; Gonçalves et al., 2009).

Climate deterioration during glacial spells from the perspective of *T. pygmaeus* in the eastern part of its range is known from climate reconstructions (Hijmans et al., 2005), and reduced habitat suitability has been predicted from species distribution modelling (Wielstra et al., 2013).

Individuals from the southwestern corner of Portugal (P1082, 5016 and P7402) form a sister clade relationship with the southeastern remainder of the group, from which they are separated by the Guadiana river. Our results show that individuals from the Betic region (P7643, P7651, P7663 and P7698) are nested within the eastern subgroup of *T. pygmaeus*, along with an individual from Doñana National Park (P4199 from north of the Guadalquivir, see Arntzen, 2023 for a possible explanation) and two individuals (P7809 and P7722) from the source area of that river. Genetically distinct intraspecific groups in the Betic region have been observed in other amphibians (*Salamandra salamandra longirostris*, Antunes et al., 2018; Pleurodeles waltl, Gutiérrez-Rodríguez et al., 2017a; and *Discoglossus galganoi* jeanneae, Dufresnes et al., 2020c). Southern *T. pygmaeus* are also morphologically distinct, differing in body size and colouration patterning, and have recently been described as subspecifically distinct (Arntzen, 2023, 2018; García-París et al., 1993). However, the phylogenetic position of two individuals from north of the Guadalquivir (P6751 and P6537) suggests that, as currently defined, the nominotypical subspecies of *T. pygmaeus* (encompassing the Betic *T. pygmaeus*) and *T. pygmaeus lusitanicus* (encompassing the remainder of the range) are not reciprocally monophyletic, suggesting the taxonomy of *T. pygmaeus* requires further study.

**5. Conclusion**

Our study provides new insights into the phylogeography of an Iberian amphibian species pair, the marbled newts, by employing a large number of nuclear markers obtained through target enrichment by sequence capture. Although the costs per individual with this method still limit sampling density (both in terms of the number of localities and the number of individuals per locality), we show target enrichment by sequence capture to be a robust approach for next-generation phylo- geography. The observed lack of introgression corroborates species status of *T. marmoratus* and *T. pygmaeus* and both species show sub- stantial genetic structuring. Comparative phylogeography of the marbled newts and other Iberian amphibians highlights the role of the Iberian Peninsula as a hotspot for inter- and intraspecific differentiation.

CRediT authorship contribution statement

**Christos Kazilas:** Resources, Writing – review & editing. **Christophe Dufresnes:** Resources, Writing – review & editing. **James France:** Investigation, Methodology, Writing – review & editing. **Konstantinos Kalaentzis:** Formal analysis, Writing – review & editing. **Inigo Martínez-Solano:** Resources, Writing – review & editing. **Manon C. de Visser:** Investigation, Methodology, Writing – review & editing. **Jan W. Arntzen:** Conceptualization, Formal analysis, Methodology, Resources, Writing – original draft. **Ben Wielstra:** Conceptualization, Formal analysis, Methodology, Resources, Supervision, Writing – original draft, Writing – review & editing.

**Declaration of competing interest**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Data availability**

The Illumina sequencing reads generated for this study have been
submitted to the NCBI Sequence Read Archive (SRA) and can be
retrieved through BioProject PRJNA1081190. All scripts utilized can be
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Phylogenetic tree for marbled newts built with RAxML based on an alignment of 122,989 SNPs representing 6,884 nuclear DNA markers. *Triturus marmoratus* individuals are shown in light and dark blue and *T. pygmaeus* individuals in yellow and orange. Outgroup individuals of *T. cristatus* are shown in grey. Some long branches were truncated, as shown by double slashes (//). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

from central to southern Portugal (Fig. 2). A similar pattern of genetic differentiation is observed in other amphibian species in this region (e.g. *Pelodytes atlanticus*, Díaz-Rodríguez et al., 2017; and *Lissotriton boscai*, Peñalver-Alcázar et al., 2021; Sequeira et al., 2020; Teixeira et al., 2015). The limited data available so far suggest that the Tagus river may be separating the two groups. Alternatively, regional uplifts linked to Pleistocene tectonic activities in southwestern Iberia (Figueiredo et al., 2014), combined with Quaternary sea level changes, may have locally promoted population isolation.

There is little differentiation in the eastern part of the range, pre-
sumably due to postglacial expansion from a southwestern refugium to
the central parts of the Iberian Peninsula (Arntzen, 2018; Arntzen and
Espregueira Themudo, 2008; Espregueira Themudo et al., 2012; Espre-
gueira Themudo and Arntzen, 2007a; Wielstra et al., 2013). Such phy-
logeographic homogeneity is also tangible from other amphibians,
including *Alytes cisternasii* (Ambu et al., 2023; Gonçalves et al., 2009).

Climate deterioration during glacial spells from the perspective of *T. pygmaeus* in the eastern part of its range is known from climate reconstructions (Hijmans et al., 2005), and reduced habitat suitability has been predicted from species distribution modelling (Wielstra et al., 2013).

Individuals from the southwestern corner of Portugal (P1082, 5016 and P7402) form a sister clade relationship with the southeastern remainder of the group, from which they are separated by the Guadiana river. Our results show that individuals from the Betic region (P7643, P7651, P7663 and P7698) are nested within the eastern subgroup of *T. pygmaeus*, along with an individual from Doñana National Park (P4199 from north of the Guadalquivir, see Arntzen, 2023 for a possible explanation) and two individuals (P7809 and P7722) from the source area of that river. Genetically distinct intraspecific groups in the Betic region have been observed in other amphibians (*Salamandra salamandra longirostris*, Antunes et al., 2018; Pleurodeles waltl, Gutiérrez-Rodríguez et al., 2017a; and *Discoglossus galganoi* jeanneae, Dufresnes et al., 2020c). Southern *T. pygmaeus* are also morphologically distinct, differing in body size and colouration patterning, and have recently been described as subspecifically distinct (Arntzen, 2023, 2018; García-París et al., 1993). However, the phylogenetic position of two individuals from north of the Guadalquivir (P6751 and P6537) suggests that, as currently defined, the nominotypical subspecies of *T. pygmaeus* (encompassing the Betic *T. pygmaeus*) and *T. pygmaeus lusitanicus* (encompassing the remainder of the range) are not reciprocally monophyletic, suggesting the taxonomy of *T. pygmaeus* requires further study.

**5. Conclusion**

Our study provides new insights into the phylogeography of an Iberian amphibian species pair, the marbled newts, by employing a large number of nuclear markers obtained through target enrichment by sequence capture. Although the costs per individual with this method still limit sampling density (both in terms of the number of localities and the number of individuals per locality), we show target enrichment by sequence capture to be a robust approach for next-generation phylo- geography. The observed lack of introgression corroborates species status of *T. marmoratus* and *T. pygmaeus* and both species show sub- stantial genetic structuring. Comparative phylogeography of the marbled newts and other Iberian amphibians highlights the role of the Iberian Peninsula as a hotspot for inter- and intraspecific differentiation.

CRediT authorship contribution statement

**Christos Kazilas:** Resources, Writing – review & editing. **Christophe Dufresnes:** Resources, Writing – review & editing. **James France:** Investigation, Methodology, Writing – review & editing. **Konstantinos Kalaentzis:** Formal analysis, Writing – review & editing. **Inigo Martínez-Solano:** Resources, Writing – review & editing. **Manon C. de Visser:** Investigation, Methodology, Writing – review & editing. **Jan W. Arntzen:** Conceptualization, Formal analysis, Methodology, Resources, Writing – original draft. **Ben Wielstra:** Conceptualization, Formal analysis, Methodology, Resources, Supervision, Writing – original draft, Writing – review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

The Illumina sequencing reads generated for this study have been
submitted to the NCBI Sequence Read Archive (SRA) and can be
retrieved through BioProject PRJNA1081190. All scripts utilized can be
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Time-calibrated phylogenies created with BEAST for 216,192 bp of nuDNA (a) and 658 bp of mtDNA (b). *Triturus marmoratus* individuals are shown in light and dark blue and *T. pygmaeus* individuals in yellow and orange in (a) and mtDNA haplotypes in (b) are colored according to the nuDNA-based composition of the individuals they were recovered from. The outgroup is not shown. Sample numbers and haplotype codes correspond to Supplementary Table 1.</figcaption></figure>

found in the following GitHub repository: https://github.com/Wielstra-Lab/marbled_newts.

## Acknowledgements

Licenses to collect were provided as follows: France - Prefecture de la Mayenne, by permit number 2003-A-2007; Spain - various provinces by permits numbers CN0010/12/ACA, CN03/0085, CN04/0269, CN10/0030, DGMEN/SEN/avp_12_015_aut, DNP 27/2008, E.P-107/04 (MG), IS/pa/EPCYL/129/2012 and SGYB/FOA/AFR/CFS, Consejerías de Medio Ambiente de Castilla la Mancha, Castilla y León, Comunidad de Madrid (Grant PID2020-116289GB-100, FEDER/Ministerio de Ciencia, Innovación y Universidades-Agencia Estatal de Investigación, Spain), and Portugal - Instituto da Conservação da Natureza, by letters dated 26/10/1998, 19/4/2000 and 19/3/2002 and by permit numbers 397/2007/CAPT, 102/2010/CAPT, 103/2010/CAPT and 107/2012/CAPT. Sampling in the Doñana National Park was carried out on the basis of the approved project proposal 'Phylogeography of marbled newts (*Triturus marmoratus* and *T. pygmaeus*)', with permit number 27/2008. We thank Alberto Gil Chamorro from the Junta de Extremadura, Consejería de Medio Ambiente y Rural, Políticas Agrarias y Territorio, Dirección General de Medio Ambiente, Dirección Banco Genético de la Bio-diversidad de Extremadura for sample numbers BGBD-EX2773D (7961) and BGBD-EX2867D (7959). Other tissue samples were made available by Salvador Carranza, Daniel Fernandez, Julia Lopez-Delgado and Ernesto Recuero. KK was the recipient of an Onassis Foundation scholarship.

## Appendix A. Supplementary material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.ympev.2024.108043.

## References

* Abellán, P., Svenning, J.C., 2014. Refugia within refugia - patterns in endemism and genetic divergence are linked to Late Quaternary climate stability in the Iberian Peninsula. Biol. J. Linn. Soc. 113, 13–28. https://doi.org/10.1111/bij.12309.

* Achalhi, M., Münch, P., Cornée, J.-J., Azdimousa, A., Melinte-Dobrinescu, M., Quillévéré, F., Drinia, H., Fauquette, S., Jiménez-Moreno, G., Merzeraud, G., Moussa, A.B., El Kharim, Y., Feddi, N., 2016. The late Miocene Mediterranean-Atlantic connections through the North Rifian Corridor: new insights from the Boudinar and Arbaa Taourirt basins (northeastern Rif, Morocco). Palaeogeogr. Palaeoclimatol. Palaeoecol. 459, 131–152. https://doi.org/10.1016/j.palaeo.2016.06.040.

* Alexander, D.H., Novembre, J., Lange, K., 2009. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664. https://doi.org/10.1101/gr.094052.109.

* Ambu, J., Martínez-Solano, Í., Suchan, T., Hernandez, A., Wielstra, B., Crochet, P.A., Dufresnes, C., 2023. Genomic phylogeography illuminates deep cyto-nuclear discordances in midwife toads (*Alytes*). Mol. Phylogenet. Evol. 183, 107783 https://doi.org/10.1016/j.ympev.2023.107783.

* Andermann, T., Torres Jiménez, M.F., Matos-Maraví, P., Batista, R., Blanco-Pastor, J.L., Gustafsson, A.L.S., Kistler, L., Liberal, I.M., Oxelman, B., Bacon, C.D., Antonelli, A., 2020. A guide to carrying out a phylogenomic target sequence capture project. Front. Genet. 10, 1–20. https://doi.org/10.3389/fgene.2019.01407.

* Anderson, E.C., Thompson, E.A., 2002. A model-based method for identifying species hybrids using multilocus genetic data. Genetics 160, 1217–1229. https://doi.org/10.1093/genetics/160.3.1217.

* Andrews, S., 2010. FastQC. A quality control tool for high throughput sequence data [WWW Document]. URL http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed 7.7.23).

* Antunes, B., Lourenço, A., Caeiro-Dias, G., Dinis, M., Gonçalves, H., Martínez-Solano, I., Tarroso, P., Velo-Antón, G., 2018. Combining phylogeography and landscape genetics to infer the evolutionary history of a short-range Mediterranean relic,
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Time-calibrated phylogenies created with BEAST for 216,192 bp of nuDNA (a) and 658 bp of mtDNA (b). *Triturus marmoratus* individuals are shown in light and dark blue and *T. pygmaeus* individuals in yellow and orange in (a) and mtDNA haplotypes in (b) are colored according to the nuDNA-based composition of the individuals they were recovered from. The outgroup is not shown. Sample numbers and haplotype codes correspond to Supplementary Table 1.</figcaption></figure>

found in the following GitHub repository: https://github.com/Wielstra-Lab/marbled_newts.

## Acknowledgements

Licenses to collect were provided as follows: France - Prefecture de la Mayenne, by permit number 2003-A-2007; Spain - various provinces by permits numbers CN0010/12/ACA, CN03/0085, CN04/0269, CN10/0030, DGMEN/SEN/avp_12_015_aut, DNP 27/2008, E.P-107/04 (MG), IS/pa/EPCYL/129/2012 and SGYB/FOA/AFR/CFS, Consejerías de Medio Ambiente de Castilla la Mancha, Castilla y León, Comunidad de Madrid (Grant PID2020-116289GB-100, FEDER/Ministerio de Ciencia, Innovación y Universidades-Agencia Estatal de Investigación, Spain), and Portugal - Instituto da Conservação da Natureza, by letters dated 26/10/1998, 19/4/2000 and 19/3/2002 and by permit numbers 397/2007/CAPT, 102/2010/CAPT, 103/2010/CAPT and 107/2012/CAPT. Sampling in the Doñana National Park was carried out on the basis of the approved project proposal 'Phylogeography of marbled newts (*Triturus marmoratus* and *T. pygmaeus*)', with permit number 27/2008. We thank Alberto Gil Chamorro from the Junta de Extremadura, Consejería de Medio Ambiente y Rural, Políticas Agrarias y Territorio, Dirección General de Medio Ambiente, Dirección Banco Genético de la Bio-diversidad de Extremadura for sample numbers BGBD-EX2773D (7961) and BGBD-EX2867D (7959). Other tissue samples were made available by Salvador Carranza, Daniel Fernandez, Julia Lopez-Delgado and Ernesto Recuero. KK was the recipient of an Onassis Foundation scholarship.

## Appendix A. Supplementary material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.ympev.2024.108043.

## References

* Abellán, P., Svenning, J.C., 2014. Refugia within refugia - patterns in endemism and genetic divergence are linked to Late Quaternary climate stability in the Iberian Peninsula. Biol. J. Linn. Soc. 113, 13–28. https://doi.org/10.1111/bij.12309.

* Achalhi, M., Münch, P., Cornée, J.-J., Azdimousa, A., Melinte-Dobrinescu, M., Quillévéré, F., Drinia, H., Fauquette, S., Jiménez-Moreno, G., Merzeraud, G., Moussa, A.B., El Kharim, Y., Feddi, N., 2016. The late Miocene Mediterranean-Atlantic connections through the North Rifian Corridor: new insights from the Boudinar and Arbaa Taourirt basins (northeastern Rif, Morocco). Palaeogeogr. Palaeoclimatol. Palaeoecol. 459, 131–152. https://doi.org/10.1016/j.palaeo.2016.06.040.

* Alexander, D.H., Novembre, J., Lange, K., 2009. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664. https://doi.org/10.1101/gr.094052.109.

* Ambu, J., Martínez-Solano, Í., Suchan, T., Hernandez, A., Wielstra, B., Crochet, P.A., Dufresnes, C., 2023. Genomic phylogeography illuminates deep cyto-nuclear discordances in midwife toads (*Alytes*). Mol. Phylogenet. Evol. 183, 107783 https://doi.org/10.1016/j.ympev.2023.107783.

* Andermann, T., Torres Jiménez, M.F., Matos-Maraví, P., Batista, R., Blanco-Pastor, J.L., Gustafsson, A.L.S., Kistler, L., Liberal, I.M., Oxelman, B., Bacon, C.D., Antonelli, A., 2020. A guide to carrying out a phylogenomic target sequence capture project. Front. Genet. 10, 1–20. https://doi.org/10.3389/fgene.2019.01407.

* Anderson, E.C., Thompson, E.A., 2002. A model-based method for identifying species hybrids using multilocus genetic data. Genetics 160, 1217–1229. https://doi.org/10.1093/genetics/160.3.1217.

* Andrews, S., 2010. FastQC. A quality control tool for high throughput sequence data [WWW Document]. URL http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed 7.7.23).

* Antunes, B., Lourenço, A., Caeiro-Dias, G., Dinis, M., Gonçalves, H., Martínez-Solano, I., Tarroso, P., Velo-Antón, G., 2018. Combining phylogeography and landscape genetics to infer the evolutionary history of a short-range Mediterranean relic,
Salamandra salamandra longirostris. Conserv. Genet. 19, 1411–1424. https://doi.org/10.1007/s10592-018-1107-10.

Arntzen, J.W., 2018. Morphological and molecular characters to describe a marbled newt hybrid zone in the Iberian peninsula. Contrib. to Zool. 87, 167–185. https://doi.org/10.1163/18759866-08703003.

Arntzen, J.W., 2023. Morphological and genetic diversification of Old-World marbled news, with the description of a new and 'not-at-all-cryptic' subspecies from the Iberian Peninsula (Triturus, Salamandridae). Contrib. to Zool. 1–26 https://doi.org/10.1163/18759866-bja10055.

Arntzen, J.W., Espregueira Themudo, G., 2008. Environmental parameters that determine species geographical range limits as a matter of time and space. J. Biogr. 35, 1177–1186. https://doi.org/10.1111/j.1365-2699.2007.01875.x.

Arntzen, J.W., Espregueira Themudo, G., Wielstra, B., 2007. The phylogeny of crested news (Triturus cristatus superspecies): nuclear and mitochondrial genetic characters suggest a hard polytypy, in line with the paleogeography of the centre of origin. Contrib. to Zool. 76, 261–278. https://doi.org/10.1163/18759866-07604005.

Arntzen, J.W., Wielstra, B., Wallis, G.P., 2014. The modality of nine Triturus newt hybrid zones assessed with nuclear, mitochondrial and morphological data. Biol. J. Linn. Soc. 113, 604–622. https://doi.org/10.1111/bij.12358.

Arntzen, J.W., de Vries, W., Canestrelli, D., Martínez-Solano, I., 2017. Hybrid zone formation and contrasting outcomes of secondary contact over transsects in common toads. Mol. Ecol. 26, 5663–5675. https://doi.org/10.1111/mec.14273.

Arntzen, J.W., López-Delgado, J., van Riemsdijk, I., Wielstra, B., 2020. A genomic footprint of a moving hybrid zone in marbled news. J. Zool. Syst. Evol. Res. 59, 459–465. https://doi.org/10.1111/jzs.12439.

Avise, J.C., 2000. Phylogeography: The History and Formation of Species. Harvard University Press.

Bailon, S., 2003. Le Pélobate cultripède: Pelobates cultripes (Cuvier, 1829). In: Pascal, M., Lorvelec, O., Vigne, J.-D., Keith, P., Clerguea, P. (Eds.), Évolution Holocène de La Faune de Vertébrés de France: Invasions et Disparitions. Rapport INRA/CNRS/ MNHN, pp. 147–148.

Ballard, J.W.O., Whitlock, M.C., 2004. The incomplete natural history of mitochondria. Mol. Ecol. 13, 729–744. https://doi.org/10.1046/j.1365-294X.2003.02063.x.

Beheregaray, L.B., 2008. Twenty years of phylogeography: the state of the field and the challenges for the Southern Hemisphere. Mol. Ecol. 17, 3754–3774. https://doi.org/10.1111/j.1365-294X.2008.03857.x.

Booth-Rea, G., Ranero, R., Grevemeyer, I., 2018. The Alboran volcanic-arc modulated the Messianian faunal exchange and salinity crisis. Sci. Rep. 8, 13015. https://doi.org/10.1038/s41598-018-31307-7.

Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C.H., Xie, D., Suchard, M.A., Rambaut, A., Drummond, A.J., 2014. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 10, 1–7. https://doi.org/10.1371/journal.pcbi.1003537.

Busack, S.D., Jericho, B.G., Maxson, L.R., Uzzell, P., 1988. Evolutionary relationships of salamanders in the genus Triturus: the view from immunology 44, 307–316.

Cuttelod, A., García, N., Malak, D.A., Temple, J.H., Katariya, V., 2009. The Mediterranean: a biodiversity hotspot under threat, Wildlife in a Changing World – an analysis of the 2008 IUCN Red List of Threatened Species. doi: 10.1111/j.1439-0485.2010.00364.x.

Danecek, P., Bonfield, J.K., Liddle, J., Marshall, J., Ohan, V., Pollard, M.O., Whitwham, A., Keane, T., McCarthy, S.A., Davies, R.M., Li, H., 2021. Twelve years of SAMtools and BCFtools. GigaScience 10, 1–4. https://doi.org/10.1093/gigascience/giab008.

Davey, J.W., Hohenlohe, P.A., Etter, P.D., Boone, J.Q., Catchen, J.M., Blaxter, M.L., 2011. Genome-wide genetic marker discovery and genotyping using next-generation sequencing. Nat. Rev. Genet. 12, 499–510. https://doi.org/10.1038/nrg3012.

Díaz-Rodríguez, J., Gonçalves, H., Sequeira, F., Sousa-Neves, T., Tejedo, M., Ferrand, N., Martínez-Solano, I., 2015. Molecular evidence for cryptic candidate species in Iberian Pelodytes (Anura, Pelodytidae). Mol. Phylogenet. Evol. 83, 224–241. https://doi.org/10.1016/j.ympev.2014.12.007.

Díaz-Rodríguez, J., Gehara, M., Márquez, R., Vences, M., Gonçalves, H., Sequeira, F., Martínez-Solano, I., Tejedo, M., 2017. Integration of molecular, bioacoustical and morphological data reveals two new cryptic species of Pelodytes (Anura, Pelodytidae) from the Iberian Peninsula. Zootaxa 4243, 1–41. 10.11646/zootaxa.4243.1.1.

Dorda, J., Esteban, M., 1986. Evolution and subspeciation of Iberian Triturus marmoratus. Stud. Herpetol. 159–164.

Dufresnes, C., 2019. Amphibians of Europe, North Africa and the Middle East: a photographic guide. Bloomsbury, London.

Dufresnes, C., Di Santo, L., Leuenberger, J., Schuerch, J., Mazepa, G., Grandjean, N., Canestrelli, D., Perrin, N., Dubey, S., 2017. Cryptic invasion of Italian pool frogs (Pelophylax bergeri) across Western Europe unraveled by multilocus phylogeography. Biol. Invasions 19, 1407–1420. https://doi.org/10.1007/s10530-016-1359-z.

Dufresnes, C., Berroneau, M., Dubey, S., Litvinchuk, S.N., Perrin, N., 2020a. The effect of phylogeographic history on species boundaries: a comparative framework in *Hyla* tree frogs. Sci. Rep. 10, 5502. https://doi.org/10.1038/s41598-020-62382-4.

Dufresnes, C., Hernandez, A., 2021. Phylogeographic advances in midwife toads (*Alytes*) support the existence of a novel taxon endemic to the Central Pyrenees. J. Zool. Syst. Evol. Res. 59, 2170–2179. https://doi.org/10.1111/jzs.12564.

Dufresnes, C., Jablonski, D., 2022. A genomics revolution in amphibian taxonomy. Science 377, 1272. https://doi.org/10.1126/science.ade5002.

Dufresnes, C., Strachinis, I., Suriadna, N., Mykytynets, G., Cogălniceanu, D., Székely, P., Vukov, T., Arntzen, J.W., Wielstra, B., Lymberakis, P., Geffen, E., Gafny, S., Kumlutaş, Y., Ilgaz, Ç., Candan, K., Mizsei, E., Szabolcs, M., Kolenda, K., Smirnov, N., Géniez, P., Lukanov, S., Crochet, P.A., Dubey, S., Perrin, N., Litvinchuk, S.N., Denoël, M., 2019. Phylogeography of a cryptic speciation continuum in Eurasian

spadefoot toads (Pelobates). Mol. Ecol. 28, 3257–3270. https://doi.org/10.1111/mec.15133.

Dufresnes, C., Nicieza, A.G., Litvinchuk, S.N., Rodrigues, N., Jeffries, D.L., Vences, M., Perrin, N., Martínez-Solano, I., 2020b. Are glacial refugia hotspots of speciation and cytonuclear discordances? Answers from the genomic phylogeography of Spanish common frogs. Mol. Ecol. 29, 986–1000. https://doi.org/10.1111/mec.15368.

Dufresnes, C., Pribille, M., Alard, B., Gonçalves, H., Amat, F., Crochet, P.A., Dubey, S., Perrin, N., Fumagalli, L., Vences, M., Martínez-Solano, I., 2020c. Integrating hybrid zone analyses in species delimitation: lessons from two anuran radiations of the Western Mediterranean. Heredity 124, 423–438. https://doi.org/10.1038/s41437-020-0294-z.

Dufresnes, C., Poyarkov, N., Jablonski, D., 2023. Acknowledging more biodiversity without more species. Proc. Natl. Acad. Sci. 120, e2302424120 https://doi.org/10.1073/pnas.2302424120.

Ehl, S., Vences, M., Veith, M., 2019. Reconstructing evolution at the community level: a case study on Mediterranean amphibians. Mol. Phylogenet. Evol. 134, 211–225. https://doi.org/10.1016/j.ympev.2019.02.013.

Espregueira Themudo, G., Arntzen, J.W., 2007a. Newts under siege: Range expansion of Triturus pygmaeus isolates populations of its sister species. Divers. Distrib. 13, 580–586. https://doi.org/10.1111/j.1472-4642.2007.00373.x.

Espregueira Themudo, G., Arntzen, J.W., 2007b. Molecular identification of marbled news and a justification of species status for Triturus marmoratus and T. pygmaeus. Herpetol. J. 17, 24–30.

Espregueira Themudo, G., Nieman, A.M., Arntzen, J.W., 2012. Is dispersal guided by the environment? A comparison of interspecific gene flow estimates among differentiated regions of a newt hybrid zone. Mol. Ecol. 21, 5324–5335. https://doi.org/10.1011/mec.12026.

Figueiredo, P.M., Cabral, J., Rockwell, T.K., 2014. Recognition of Pleistocene marine terraces in the southwest of Portugal (Iberian Peninsula): evidences of regional Quaternary uplift. Ann. Geophys. 56 https://doi.org/10.4401/ag-6276.

Fijarczyk, A., Nadachowska, K., Hofman, S., Litvinchuk, S.N., Babik, W., Stuglik, M., Gollmann, G., Choleva, L., Cogălniceanu, D., Vukov, T., Dzukić, G., Szymura, J.M., 2011. Nuclear and mitochondrial phylogeography of the European fire-bellied toads Bombina bombina and Bombina variegata supports their independent histories. Mol. Ecol. 20, 3381–3398. https://doi.org/10.1111/j.1365-294X.2011.05175.x.

Frost, D.R., Grant, T., Faivovich, J., Bain, R.H., Haas, A., Haddad, C.F.B., De Sá, R.O., Channing, A., Wilkinson, M., Donnellan, S.C., Raxworthy, C.J., Campbell, J.A., Blotto, B.L., Moler, P., Drewes, R.C., Nussbaum, R.A., Lynch, J.D., Green, D.M., Wheeler, W.C., 2006. The amphibian tree of life. Bull. Am. Museum Nat. Hist. 297, 1–370.

García-París, M., Herrero, P., Martín, C., Dorda, J., Esteban, M., Arano, B., 1993. Morphological characterization, cytogenetic analysis, and geographical distribution of the Pygmy marbled newt Triturus marmoratus pygmaeus (Wolterstorff, 1905) (Caudata: Salamandridae). Bijdr. Tot Dierkd. 63, 3–14.

García-París, M., Arano, B., Herrero, P., 2001. Molecular characterization of the contact zone between Triturus pygmaeus and *T. marmoratus* (Caudata: Salamandridae) in Central Spain and their taxonomic assessment. Rev. Española Herpetol. 15, 115–126.

Gernandt, D.S., Aguirre Dugua, X., Vázquez-Lobo, A., Willyard, A., Moreno Letelier, A., Pérez de la Rosa, J.A., Pinero, D., Liston, A., 2018. Multi-locus phylogenetics, lineage sorting, and reticulation in *Pinus* subsection Australes. Am. J. Bot. 105, 711–725. https://doi.org/10.1002/ajb2.1052.

Gómez, A., Lunt, D.H., 2007. Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. In: Phylogeography of Southern European Refugia. Springer, Dordrecht, pp. 155–188.

Gonçalves, H., Martínez-Solano, I., Pereira, R.J., Carvalho, B., García-París, M., Ferrand, N., 2009. High levels of population subdivision in a morphologically conserved Mediterranean toad (*Alytes cisternax*) result from recent, multiple refugia: evidence from mtDNA, microsatellites and nuclear genealogies. Mol. Ecol. 18, 5143–5160. https://doi.org/10.1111/j.1365-294X.2009.04426.x.

Gonçalves, H., Maia-Carvalho, B., Sousa-Neves, T., García-París, M., Sequeira, F., Ferrand, N., Martínez-Solano, I., 2015. Multilocus phylogeography of the common midwife toad, *Alytes obstetricans* (Anura, Alytidae): contrasting patterns of lineage diversification and genetic structure in the Iberian refugium. Mol. Phylogenet. Evol. 93, 363–379. https://doi.org/10.1016/j.ympev.2015.08.009.

Gutiérrez-Rodríguez, J., Barbosa, A.M., Martínez-Solano, I., 2017a. Integrative inference of population history in the Ibero-Maghrebian endemic *Pleurodeles waltli* (Salamandridae). Mol. Phylogenet. Evol. 112, 122–137. https://doi.org/10.1016/j.ympev.2017.04.022.

Gutiérrez-Rodríguez, J., Barbosa, A.M., Martínez-Solano, I., 2017b. Present and past climatic effects on the current distribution and genetic diversity of the Iberian spadefoot toad (Pelobates cultripes): an integrative approach. J. Biogeogr. 44, 245–258. https://doi.org/10.1111/jbi.12791.

Hewitt, G., 2000. The genetic legacy of the quaternary ice ages. Nature 405, 907–913. https://doi.org/10.1038/35016000.

Hewitt, G., 2011a. Mediterranean peninsulas: the evolution of hotspots. Biodiversity Hotspots. 123–147. https://doi.org/10.1007/978-3-642-20992-5_2.

Hewitt, G., 2011b. Quaternary phylogeography: the roots of hybrid zones. Genetica 139, 617–638. https://doi.org/10.1007/s10709-011-9547-3.

Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. https://doi.org/10.1002/joc.1276.

Huang, Z.F., Chiba, H., Jin, J., Kizhakke, A.G., Wang, M., Kunte, K., Fan, X.L., 2019. A multilocus phylogenetic framework of the tribe Aeromachini (Lepidoptera: Hesperiidae: Hesperiinae), with implications for taxonomy and biogeography. Syst. Entomol. 44, 163–178. https://doi.org/10.1111/syen.12322.
Salamandra salamandra longirostris. Conserv. Genet. 19, 1411–1424. https://doi.org/10.1007/s10592-018-1107-10.

Arntzen, J.W., 2018. Morphological and molecular characters to describe a marbled newt hybrid zone in the Iberian peninsula. Contrib. to Zool. 87, 167–185. https://doi.org/10.1163/18759866-08703003.

Arntzen, J.W., 2023. Morphological and genetic diversification of Old-World marbled news, with the description of a new and 'not-at-all-cryptic' subspecies from the Iberian Peninsula (*Triturus*, *Salamandridae*). Contrib. to Zool. 1–26 https://doi.org/10.1163/18759866-bja10055.

Arntzen, J.W., Espregueira Themudo, G., 2008. Environmental parameters that determine species geographical range limits as a matter of time and space. J. Biogr. 35, 1177–1186. https://doi.org/10.1111/j.1365-2699.2007.01875.x.

Arntzen, J.W., Espregueira Themudo, G., Wielstra, B., 2007. The phylogeny of crested news (Triturus cristatus superspecies): nuclear and mitochondrial genetic characters suggest a hard polytypy, in line with the paleogeography of the centre of origin. Contrib. to Zool. 76, 261–278. https://doi.org/10.1163/18759866-07604005.

Arntzen, J.W., Wielstra, B., Wallis, G.P., 2014. The modality of nine *Triturus* newt hybrid zones assessed with nuclear, mitochondrial and morphological data. Biol. J. Linn. Soc. 113, 604–622. https://doi.org/10.1111/bij.12358.

Arntzen, J.W., de Vries, W., Canestrelli, D., Martínez-Solano, I., 2017. Hybrid zone formation and contrasting outcomes of secondary contact over transsects in common toads. Mol. Ecol. 26, 5663–5675. https://doi.org/10.1111/mec.14273.

Arntzen, J.W., López-Delgado, J., van Riemsdijk, I., Wielstra, B., 2020. A genomic footprint of a moving hybrid zone in marbled news. J. Zool. Syst. Evol. Res. 59, 459–465. https://doi.org/10.1111/jzs.12439.

Avise, J.C., 2000. Phylogeography: The History and Formation of Species. Harvard University Press.

Bailon, S., 2003. Le Pélobate cultripède: *Pelobates cultripes* (Cuvier, 1829). In: Pascal, M., Lorvelec, O., Vigne, J.-D., Keith, P., Clergeau, P. (Eds.), Évolution Holocène de La Faune de Vertébrés de France: Invasions et Disparitions. Rapport INRA/CNRS/ MNHN, pp. 147–148.

Ballard, J.W.O., Whitlock, M.C., 2004. The incomplete natural history of mitochondria. Mol. Ecol. 13, 729–744. https://doi.org/10.1046/j.1365-294X.2003.02063.x.

Beheregaray, L.B., 2008. Twenty years of phylogeography: the state of the field and the challenges for the Southern Hemisphere. Mol. Ecol. 17, 3754–3774. https://doi.org/10.1111/j.1365-294X.2008.03857.x.

Booth-Rea, G., Ranero, R., Grevemeyer, I., 2018. The Alboran volcanic-arc modulated the Messianian faunal exchange and salinity crisis. Sci. Rep. 8, 13015. https://doi.org/10.1038/s41598-018-31307-7.

Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C.H., Xie, D., Suchard, M.A., Rambaut, A., Drummond, A.J., 2014. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 10, 1–7. https://doi.org/10.1371/journal.pcbi.1003537.

Busack, S.D., Jericho, B.G., Maxson, L.R., Uzzell, P., 1988. Evolutionary relationships of salamanders in the genus *Triturus*: the view from immunology 44, 307–316.

Cuttelod, A., García, N., Malak, D.A., Temple, J.H., Katariya, V., 2009. The Mediterranean: a biodiversity hotspot under threat, Wildlife in a Changing World – an analysis of the 2008 IUCN Red List of Threatened Species. doi: 10.1111/j.1439-0485.2010.00364.x.

Danecek, P., Bonfield, J.K., Liddle, J., Marshall, J., Ohan, V., Pollard, M.O., Whitwham, A., Keane, T., McCarthy, S.A., Davies, R.M., Li, H., 2021. Twelve years of SAMtools and BCFtools. GigaScience 10, 1–4. https://doi.org/10.1093/gigascience/giab008.

Davey, J.W., Hohenlohe, P.A., Etter, P.D., Boone, J.Q., Catchen, J.M., Blaxter, M.L., 2011. Genome-wide genetic marker discovery and genotyping using next-generation sequencing. Nat. Rev. Genet. 12, 499–510. https://doi.org/10.1038/nrg3012.

Díaz-Rodríguez, J., Gonçalves, H., Sequeira, F., Sousa-Neves, T., Tejedo, M., Ferrand, N., Martínez-Solano, I., 2015. Molecular evidence for cryptic candidate species in Iberian *Pelodytes* (Anura, Pelodytidae). Mol. Phylogenet. Evol. 83, 224–241. https://doi.org/10.1016/j.ympev.2014.12.007.

Díaz-Rodríguez, J., Gehara, M., Márguez, R., Vences, M., Gonçalves, H., Sequeira, F., Martínez-Solano, I., Tejedo, M., 2017. Integration of molecular, bioacoustical and morphological data reveals two new cryptic species of *Pelodytes* (Anura, Pelodytidae) from the Iberian Peninsula. Zootaxa 4243, 1–41. 10.11646/zootaxa.4243.1.1.

Dorda, J., Esteban, M., 1986. Evolution and subspeciation of Iberian *Triturus marmoratus*. Stud. Herpetol. 159–164.

Dufresnes, C., 2019. Amphibians of Europe, North Africa and the Middle East: a photographic guide. Bloomsbury, London.

Dufresnes, C., Di Santo, L., Leuenberger, J., Schuerch, J., Mazepa, G., Grandjean, N., Canestrelli, D., Perrin, N., Dubey, S., 2017. Cryptic invasion of Italian pool frogs (*Pelophylax bergeri*) across Western Europe unraveled by multilocus phylogeography. Biol. Invasions 19, 1407–1420. https://doi.org/10.1007/s10530-016-1359-z.

Dufresnes, C., Berroneau, M., Dubey, S., Litvinchuk, S.N., Perrin, N., 2020a. The effect of phylogeographic history on species boundaries: a comparative framework in *Hyla* tree frogs. Sci. Rep. 10, 5502. https://doi.org/10.1038/s41598-020-62382-4.

Dufresnes, C., Hernandez, A., 2021. Phylogeographic advances in midwife toads (*Alytes*) support the existence of a novel taxon endemic to the Central Pyrenees. J. Zool. Syst. Evol. Res. 59, 2170–2179. https://doi.org/10.1111/jzs.12564.

Dufresnes, C., Jablonski, D., 2022. A genomics revolution in amphibian taxonomy. Science 377, 1272. https://doi.org/10.1126/science.ade5002.

Dufresnes, C., Strachinis, I., Suriadna, N., Mykytynets, G., Cogălniceanu, D., Székely, P., Vukov, T., Arntzen, J.W., Wielstra, B., Lymberakis, P., Geffen, E., Gafny, S., Kumlutaş, Y., Ilgaz, Ç., Candan, K., Mizsei, E., Szabolcs, M., Kolenda, K., Smirnov, N., Géniez, P., Lukanov, S., Crochet, P.A., Dubey, S., Perrin, N., Litvinchuk, S.N., Denoël, M., 2019. Phylogeography of a cryptic speciation continuum in Eurasian

spadefoot toads (*Pelobates*). Mol. Ecol. 28, 3257–3270. https://doi.org/10.1111/mec.15133.

Dufresnes, C., Nicieza, A.G., Litvinchuk, S.N., Rodrigues, N., Jeffries, D.L., Vences, M., Perrin, N., Martínez-Solano, I., 2020b. Are glacial refugia hotspots of speciation and cytonuclear discorcances? Answers from the genomic phylogeography of Spanish common frogs. Mol. Ecol. 29, 986–1000. https://doi.org/10.1111/mec.15368.

Dufresnes, C., Pribille, M., Alard, B., Gonçalves, H., Amat, F., Crochet, P.A., Dubey, S., Perrin, N., Fumagalli, L., Vences, M., Martínez-Solano, I., 2020c. Integrating hybrid zone analyses in species delimitation: lessons from two anuran radiations of the Western Mediterranean. Heredity 124, 423–438. https://doi.org/10.1038/s41437-020-0294-z.

Dufresnes, C., Poyarkov, N., Jablonski, D., 2023. Acknowledging more biodiversity without more species. Proc. Natl. Acad. Sci. 120, e2302424120 https://doi.org/10.1073/pnas.2302424120.

Ehl, S., Vences, M., Veith, M., 2019. Reconstructing evolution at the community level: a case study on Mediterranean amphibians. Mol. Phylogenet. Evol. 134, 211–225. https://doi.org/10.1016/j.ympev.2019.02.013.

Espregueira Themudo, G., Arntzen, J.W., 2007a. Newts under siege: Range expansion of *Triturus pygmaeus* isolates populations of its sister species. Divers. Distrib. 13, 580–586. https://doi.org/10.1111/j.1472-4642.2007.00373.x.

Espregueira Themudo, G., Arntzen, J.W., 2007b. Molecular identification of marbled news and a justification of species status for *Triturus marmoratus* and *T. pygmaeus*. Herpetol. J. 17, 24–30.

Espregueira Themudo, G., Nieman, A.M., Arntzen, J.W., 2012. Is dispersal guided by the environment? A comparison of interspecific gene flow estimates among differentiated regions of a newt hybrid zone. Mol. Ecol. 21, 5324–5335. https://doi.org/10.1111/mec.12026.

Figueiredo, P.M., Cabral, J., Rockwell, T.K., 2014. Recognition of Pleistocene marine terraces in the southwest of Portugal (Iberian Peninsula): evidences of regional Quaternary uplift. Ann. Geophys. 56 https://doi.org/10.4401/ag-6276.

Fijarczyk, A., Nadachowska, K., Hofman, S., Litvinchuk, S.N., Babik, W., Stuglik, M., Hollmann, G., Choleva, L., Cogălniceanu, D., Vukov, T., Dzukić, G., Szymura, J.M., 2011. Nuclear and mitochondrial phylogeography of the European fire-bellied toads *Bombus bombus* and *Bombus variegata* supports their independent histories. Mol. Ecol. 20, 3381–3398. https://doi.org/10.1111/j.1365-294X.2011.05175.x.

Frost, D.R., Grant, T., Faivovich, J., Bain, R.H., Haas, A., Haddad, C.F.B., De Sá, R.O., Channing, A., Wilkinson, M., Donnellan, S.C., Raxworthy, C.J., Campbell, J.A., Blotto, B.L., Moler, P., Drewes, R.C., Nussbaum, R.A., Lynch, J.D., Green, D.M., Wheeler, W.C., 2006. The amphibian tree of life. Bull. Am. Museum Nat. Hist. 297, 1–370.

García-París, M., Herrero, P., Martín, C., Dorda, J., Esteban, M., Arano, B., 1993. Morphological characterization, cytogenetic analysis, and geographical distribution of the Pygmy marbled newt *Triturus marmoratus* *pygmaeus* (Wolterstorff, 1905) (Caudata: Salamandridae). Bijdr. Tot Dierkd. 63, 3–14.

García-París, M., Arano, B., Herrero, P., 2001. Molecular characterization of the contact zone between *Triturus pygmaeus* and *T. marmoratus* (Caudata: Salamandridae) in Central Spain and their taxonomic assessment. Rev. Española Herpetol. 15, 115–126.

Gernandt, D.S., Aguirre Dugua, X., Vázquez-Lobo, A., Willyard, A., Moreno Letelier, A., Pérez de la Rosa, J.A., Pinero, D., Liston, A., 2018. Multi-locus phylogenetics, lineage sorting, and reticulation in *Pinus* subsection Australes. Am. J. Bot. 105, 711–725. https://doi.org/10.1002/ajb2.1052.

Gómez, A., Lunt, D.H., 2007. Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. In: Phylogeography of Southern European Refugia. Springer, Dordrecht, pp. 155–188.

Gonçalves, H., Martínez-Solano, I., Pereira, R.J., Carvalho, B., García-París, M., Ferrand, N., 2009. High levels of population subdivision in a morphologically conserved Mediterranean toad (*Alytes cisternatus*) result from recent, multiple refugia: evidence from mtDNA, microsatellites and nuclear genealogies. Mol. Ecol. 18, 5143–5160. https://doi.org/10.1111/j.1365-294X.2009.04426.x.

Gonçalves, H., Maia-Carvalho, B., Sousa-Neves, T., García-París, M., Sequeira, F., Ferrand, N., Martínez-Solano, I., 2015. Multilocus phylogeography of the common midwife toad, *Alytes obstetricans* (Anura, Alytidae): contrasting patterns of lineage diversification and genetic structure in the Iberian refugium. Mol. Phylogenet. Evol. 93, 363–379. https://doi.org/10.1016/j.ympev.2015.08.009.

Gutiérrez-Rodríguez, J., Barbosa, A.M., Martínez-Solano, I., 2017a. Integrative inference of population history in the Ibero-Maghrebian endemic *Pleurodeles waltli* (Salamandridae). Mol. Phylogenet. Evol. 112, 122–137. https://doi.org/10.1016/j.ympev.2017.04.022.

Gutiérrez-Rodríguez, J., Barbosa, A.M., Martínez-Solano, I., 2017b. Present and past climatic effects on the current distribution and genetic diversity of the Iberian spadefoot toad (*Pelobates cultripes*): an integrative approach. J. Biogeogr. 44, 245–258. https://doi.org/10.1111/jbi.12791.

Hewitt, G., 2000. The genetic legacy of the quaternary ice ages. Nature 405, 907–913. https://doi.org/10.1038/35016000.

Hewitt, G., 2011a. Mediterranean peninsulas: the evolution of hotspots. Biodiversity Hotspots. 123–147. https://doi.org/10.1007/978-3-642-20992-5_2.

Hewitt, G., 2011b. Quaternary phylogeography: the roots of hybrid zones. Genetica 139, 617–638. https://doi.org/10.1007/s10709-011-9547-3.

Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. https://doi.org/10.1002/joc.1276.

Huang, Z.F., Chiba, H., Jin, J., Kizhakke, A.G., Wang, M., Kunte, K., Fan, X.L., 2019. A multilocus phylogenetic framework of the tribe Aeromachini (Lepidoptera: Hesperiidae: Hesperiinae), with implications for taxonomy and biogeography. Syst. Entomol. 44, 163–178. https://doi.org/10.1111/syen.12322.
Hurst, G.D.D., Jiggins, F.M., 2005. Problems with mitochondrial DNA as a marker in population, phylogeographic and phylogenetic studies: the effects of inherited symbionts. Proc. R. Soc. B Biol. Sci. 272, 1525–1534. https://doi.org/10.1098/rspb.2005.3056.

Hutter, C.R., Cobb, K.A., Portik, D.M., Travers, S.L., Wood, P.L., Brown, R.M., 2022. FrogCap: a modular sequence capture probe set for phylogenomics and population genetics for all frogs, assessed across multiple phylogenetic scales. Mol. Ecol. Resour. 22, 1100–1119. https://doi.org/10.1111/1755-0998.13517.

Irwin, D.E., Rubtsov, A.S., Panov, E.N., 2009. Mitochondrial introgression and replacement between yellowhammers (Emberiza citrinella) and pine buntings (Emberiza leucocephalos) (Aves: Passeriformes). Biol. J. Linn. Soc. 98, 422–438. https://doi.org/10.1111/j.1095-8312.2009.01282.x.

Jiang, H., Lei, R., Ding, S.W., Zhu, S., 2014. Skewer: a fast and accurate adapter trimmer for next-generation sequencing paired-end reads. BMC Bioinf. 15, 1–12. https://doi.org/10.1186/1471-2105-15-182.

Jones, M.R., Good, J.M., 2016. Targeted capture in evolutionary and ecological genomics. Mol. Ecol. 25, 185–202. https://doi.org/10.1111/mec.13304.

Kalaentzis, K., Arntzen, J.W., Avci, A., van den Berg, V., Beukema, W., France, J., Olgun, K., van Riemsdijk, I., Üzüm, N., de Visser, M.C., Wielstra, B., 2023. Hybrid zone analysis confirms cryptic species of banded newt and does not support competitive displacement since secondary contact. Ecol. Evol. 13, e10442. https://doi.org/10.1002/ece3.10442.

Kopelman, N.M., Mayzel, J., Jakobsson, M., Rosenberg, N.A., Mayrose, I., 2015. Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 15, 1179–1191. https://doi.org/10.1111/1755-0998.12387.

Lewis, P.O., 2001. A likelihood approach to estimating phylogeny from discrete morphological character data. Syst. Biol. 50, 913–925. https://doi.org/10.1080/163501501753462876.

Li, H., 2011. A statistical framework for SNP calling, mutation discovery, association mapping and population genetic parameter estimation from sequencing data. Bioinformatics 27, 2987–2993. https://doi.org/10.1093/bioinformatics/btr509.

Li, H., 2013. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM 00, pp. 1–3.

Lischer, H.E.L., Excoffier, L., 2012. PGDSpider: An automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics 28, 298–299. https://doi.org/10.1093/bioinformatics/btr642.

Lonergan, L., White, N., 1997. Origin of the Betic-Rif mountain belt. Tectonics 16, 504–522. https://doi.org/10.1029/96TC03937.

López-Delgado, J., van Riemsdijk, I., Arntzen, J.W., 2020. Tracing species replacement in Iberian marbled newts. Ecol. Evol. 11, 402–414. https://doi.org/10.1002/ece3.7060.

Macaluso, L., Bertini, A., Carnevale, G., Eronen, J.T., Martinetto, E., Saarinen, J., Villa, A., Capasso, F., Delfino, M., 2023. A combined palaeomodelling approach reveals the role as selective refugia of the Mediterranean peninsulas. Palaeogeogr. Palaeoclimatol. Palaeoecol. 625, 111699. https://doi.org/10.1016/j.paleo.2023.111699.

Maddison, W.P., Maddison, D.R., 2023. Mesquite: A Modular System for Evolutionary Analysis [WWW Document]. URL http://www.mesquiteproject.org (accessed 6.12.23).

Maia-Carvalho, B., Vale, C.G., Sequeira, F., Ferrand, N., Martínez-Solano, Í., Gonçalves, H., 2018. The roles of allopatric fragmentation and niche divergence in intraspecific lineage diversification in the common midwife toad (Alytes obstetricans). J. Biogeogr. 45, 2146–2158. https://doi.org/10.1111/jbi.13405.

Martin, J., 2014. Miocene Atlantic-Mediterranean Seaways in the Betic. Rev. La Soc. Geol. Espana 27, 175–186.

Martinez-Solano, I., 2004. Phylogeography of Iberian *Discoglossus* (Anura: Discoglossidae). J. Zool. Syst. Evol. Res. 42, 298–305. https://doi.org/10.1111/j.1439-0469.2004.00257.x.

Martinez-Solano, I., Teixeira, J., Buckley, D., García-París, M., 2006. Mitochondrial DNA phylogeography of *Lissotriton boscai* (Caudata, Salamandridae): evidence for old, multiple refugia in an Iberian endemic. Mol. Ecol. 15, 3375–3388. https://doi.org/10.1111/j.1365-294X.2006.03013.x.

McCartney-Melstad, E., Mount, G.G., Shaffer, H.B., 2016. Exon capture optimization in amphibians with large genomes. Mol. Ecol. Resour. 16, 1084–1094. https://doi.org/10.1111/1755-0998.12538.

McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A., Garimella, K., Altshuler, D., Gabriel, S., Daly, M., DePristo, M.A., 2010. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303. https://doi.org/10.1101/gr.107524.110.20.

Médail, F., Myers, N., 2004. Mediterranean Basin. In: Mittermeier, R.A., Robles Gil, P., Hoffman, M., Pilgrim, J., Brooks, T., Mittermeier, C.G., Lamoreux, J., Da Fonseca, G. A.B. (Eds.), Hotspots Revisited: Earth's Biologically Richest and Most Endangered Terrestrial Ecoregions. Conservation International, Washington, pp. 144–147. https://doi.org/10.1016/j.cit.2004.09.001.

Médail, F., Quezel, P., 1999. Biodiversity hotspots in the Mediterranean Basin: setting global conservation priorities. Conserv. Biol. 13, 1510–1513.

Milá, B., Carranza, S., Guillaume, O., Clobert, J., 2010. Marked genetic structuring and extreme dispersal limitation in the Pyrenean brook newt *Calotriton asper* (Amphibia: Salamandridae) revealed by genome-wide AFLP but not mtDNA. Mol. Ecol. 19, 108–120. https://doi.org/10.1111/j.1365-294X.2009.00441.x.

Miller, M.A., Pfeiffer, W., Schwartz, T., 2010. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. 2010 Gatew. Comput. Environ. Work. GCE 2010. doi: 10.1109/GCE.2010.5676129.

Myers, N., Mittermeler, R.A., Mittermeler, C.G., Da Fonseca, G.A.B., Kent, J., 2000. Biodiversity hotspots for conservation priorities. Nature 403, 853–858. https://doi.org/10.1038/35002501.

Newman, C.E., Austin, C.C., 2016. Sequence capture and next-generation sequencing of ultraconserved elements in a large-genome salamander. Mol. Ecol. 25, 6162–6174. https://doi.org/10.1111/mec.13909.

Oosterbroek, P., Arntzen, J.W., 1992. Area-cladograms of Circum-Mediterranean taxa in relation to Mediterranean palaeogeography. J. Biogeogr. 19, 3. https://doi.org/10.2307/2845616.

Pabijan, M., Zielniški, P., Dudek, K., Chloupek, M., Sotiropoulos, K., Liana, M., Babik, W., 2015. The dissection of a Pleistocene refugium: phylogeography of the smooth newt, *Lissotriton vulgaris*, in the Balkans. J. Biogeogr. 42, 671–683. https://doi.org/10.1111/jbi.12449.

Peñalver-Alcázar, M., Jiménez-Valverde, A., Aragón, P., 2021. Niche differentiation between deeply divergent phylogenetic lineages of an endemic newt: implications for species distribution models. Zoology 144. https://doi.org/10.1016/j.zool.2020.125852.

Rambaut, A., Drummond, A.J., Xie, D., Baele, G., Suchard, M.A., 2018. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67, 901–904. https://doi.org/10.1093/sysbio/syy032.

Recuero, E., Canestrelli, D., Vörös, J., Szabó, K., Poyarkov, N.A., Arntzen, J.W., Crnobrnja-Isailović, J., Kidov, A.A., Cogálniceanu, D., Caputo, F.P., Nascetti, G., Martínez-Solano, I., 2012. Multilocus species tree analyses resolve the radiation of the widespread *Bufo bufo* species group (Anura, Bufonidae). Mol. Phylogenet. Evol. 62, 71–86. https://doi.org/10.1016/j.ympev.2011.09.008.

Recuero, E., García-París, M., 2011. Evolutionary history of *Lissotriton helveticus*: multilocus assessment of ancestral vs. recent colonization of the Iberian Peninsula. Mol. Phylogenet. Evol. 60, 170–182. https://doi.org/10.1016/j.ympev.2011.04.006.

Rowe, G., Harris, D.J., Beebee, T.J.C., 2006. Lusitania revisited: a phylogeographic analysis of the natterjack toad *Bufo calamita* across its entire biogeographical range. Mol. Phylogenet. Evol. 39, 335–346. https://doi.org/10.1016/j.ympev.2005.08.021.

Sambrook, J., Russell, D.W., 2001. Molecular Cloning: A Laboratory Manual, 3rd ed. Coldspring-Harbour Laboratory Press, London, United Kingdom.

Sánchez-Montes, G., Recuero, E., Barbosa, A.M., Martínez-Solano, I., 2019. Complementing the Pleistocene biogeography of European amphibians: testimony from a southern Atlantic species. J. Biogeogr. 46, 568–583. https://doi.org/10.1111/jbi.13515.

Sequeira, F., Bessa-Silva, A., Tarroso, P., Sousa-Neves, T., Vallinoto, M., Gonçalves, H., Martínez-Solano, I., 2020. Discordant patterns of introgression across a narrow hybrid zone between two cryptic lineages of an Iberian endemic newt. J. Evol. Biol. 33, 202–216. https://doi.org/10.1111/jeb.13562.

Sequeira, F., Arntzen, J.W., van Gulik, D., Hajema, S., Diaz, R.L., Wagt, M., van Riemsdijk, I., 2022. Genetic traces of hybrid zone movement across a fragmented habitat. J. Evol. Biol. 35, 400–412. https://doi.org/10.1111/jeb.13982.

Stamatakis, A., 2014. RAXML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313. https://doi.org/10.1093/bioinformatics/btu033.

Steinfartz, S., Vicario, S., Arntzen, J.W., Caccone, A., 2007. A Bayesian approach on molecules and behavior: reconsidering phylogenetic and evolutionary patterns of the Salamandridae with emphasis on *Triturus* newts. J. Exp. Zool. Part B Mol. Dev. Evol. 308, 139–162. https://doi.org/10.1002/jez.b.21119.

Sun, C., Shepard, D.B., Chong, R.A., Arriaza, J.L., Hall, K., Castoe, T.A., Feschotte, C., Pollock, D.D., Mueller, R.L., 2012. LTR retrotransposons contribute to genomic gigantism in plethodontid salamanders. Genome Biol. Evol. 4, 168–183. https://doi.org/10.1093/gbe/evr139.

Teixeira, J., Martínez-Solano, I., Buckley, D., Tarroso, P., García-París, M., Ferrand, N., 2015. Genealogy of the nuclear *β*-fibrinogen intron 7 in *Lissotriton boscai* (Caudata, Salamandridae): concordance with mtDNA and implications for phylogeography and speciation. Contrib. to Zool. 84, 193–215. https://doi.org/10.1163/18759866-08403002.

van Riemsdijk, I., Arntzen, J.W., Babik, W., Bogaerts, S., Franzen, M., Kalaentzis, K., Litvinchuk, S.N., Olgun, K., Wijnands, J.W.P.M., Wielstra, B., 2022. Next-generation phylogeography of the banded newts (*Ommatotriton*): a phylogenetic hypothesis for three ancient species with geographically restricted interspecific gene flow and deep intraspecific genetic structure. Mol. Phylogenet. Evol. 167, 107361. https://doi.org/10.1016/j.ympev.2021.107361.

Vences, M., Hauswaldt, S.J., Steinfartz, S., Rupp, O., Goesmann, A., Künzel, S., Orozco-terWengel, P., Vieites, D.R., Nieto-Roman, S., Haas, S., Laugsch, C., Gehara, M., Bruchmann, S., Pabijan, M., Ludewig, A.-K., Rudert, D., Angelini, C., Borkin, L.J., Crochet, P.-A., Crottini, A., Dubois, A., Ficetola, G.F., Galán, P., Geniez, P., Hachtel, M., Jovanovic, O., Litvinchuk, S.N., Lymberakis, P., Ohler, A., Smirnov, N. A., 2013. Radically different phylogeographies and patterns of genetic variation in two European brown frogs, genus *Rana*. Mol. Phylogenet. Evol. 68, 657–670. https://doi.org/10.1016/j.ympev.2013.04.014.

Wickham, H., 2007. Reshaping Data with the reshape Package. J. Stat. Softw. 21, 1–20. https://doi.org/10.5325/soundings.101.1.0006.

Wickham, H., 2016. ggplot2: Elegant Graphics for Data Analysis, 2nd ed. Springer, Houston, Texas, USA. doi: 10.1007/978-3-319-24277-4.

Wielstra, B., Arntzen, J.W., 2016. Description of a new species of crested newt, previously subsumed in *Triturus ivanbureschi* (Amphibia: Caudata Salamandridae). Zootaxa 4109, 73–80. https://10.11646/zootaxa.4109.1.6.

Wielstra, B., Crnobrnja-Isailović, J., Litvinchuk, S.N., Reijnen, B.T., Skidmore, A.K., Sotiropoulos, K., Toxopeus, A.G., Tzankov, N., Vukov, T., Arntzen, J.W., 2013. Tracing glacial refugia of *Triturus* newts based on mitochondrial DNA phylogeography and species distribution modeling. Front. Zool. 10, 1–14. https://doi.org/10.1186/1742-9994-10-13.
Hurst, G.D.D., Jiggins, F.M., 2005. Problems with mitochondrial DNA as a marker in population, phylogeographic and phylogenetic studies: the effects of inherited symbionts. Proc. R. Soc. B Biol. Sci. 272, 1525-1534. https://doi.org/10.1098/rspb.2005.3056.

Hutter, C.R., Cobb, K.A., Portik, D.M., Travers, S.L., Wood, P.L., Brown, R.M., 2022. FrogCap: a modular sequence capture probe set for phylogenomics and population genetics for all frogs, assessed across multiple phylogenetic scales. Mol. Ecol. Resour. 22, 1100-1119. https://doi.org/10.1111/1755-0998.13517.

Irwin, D.E., Rubtsov, A.S., Panov, E.N., 2009. Mitochondrial introgression and replacement between yellowhammers (Emberiza citrinella) and pine buntings (Emberiza leucocephalos) (Aves: Passeriformes). Biol. J. Linn. Soc. 98, 422-438. https://doi.org/10.1111/j.1095-8312.2009.01282.x.

Jiang, H., Lei, R., Ding, S.W., Zhu, S., 2014. Skewer: a fast and accurate adapter trimmer for next-generation sequencing paired-end reads. BMC Bioinf. 15, 1-12. https://doi.org/10.1186/1471-2105-15-182.

Jones, M.R., Good, J.M., 2016. Targeted capture in evolutionary and ecological genomics. Mol. Ecol. 25, 185-202. https://doi.org/10.1111/mec.13304.

Kalaentzis, K., Arntzen, J.W., Avci, A., van den Berg, V., Beukema, W., France, J., Olgun, K., van Riemsdijk, I., Üzüm, N., de Visser, M.C., Wielstra, B., 2023. Hybrid zone analysis confirms cryptic species of banded newt and does not support competitive displacement since secondary contact. Ecol. Evol. 13, e10442. https://doi.org/10.1002/eec.31042.

Kopelman, N.M., Mayzel, J., Jakobsson, M., Rosenberg, N.A., Mayrose, I., 2015. Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 15, 1179-1191. https://doi.org/10.1111/1755-0998.12387.

Lewis, P.O., 2001. A likelihood approach to estimating phylogeny from discrete morphological character data. Syst. Biol. 50, 913-925. https://doi.org/10.1080/163501501753462876.

Li, H., 2011. A statistical framework for SNP calling, mutation discovery, association mapping and population genetic parameter estimation from sequencing data. Bioinformatics 27, 2987-2993. https://doi.org/10.1093/bioinformatics/btr509.

Li, H., 2013. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM 00, pp. 1-3.

Lischer, H.E.L., Excoffier, L., 2012. PGDSpider: An automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics 28, 298-299. https://doi.org/10.1093/bioinformatics/btr642.

Lonergan, L., White, N., 1997. Origin of the Betic-Rif mountain belt. Tectonics 16, 504-522. https://doi.org/10.1029/96TC03937.

López-Delgado, J., van Riemsdijk, I., Arntzen, J.W., 2020. Tracing species replacement in Iberian marbled newts. Ecol. Evol. 11, 402-414. https://doi.org/10.1002/eec.7060.

Macaluso, L., Bertini, A., Carnevale, G., Eronen, J.T., Martinetto, E., Saarinen, J., Villa, A., Capasso, F., Delfino, M., 2023. A combined palaeomodelling approach reveals the role as selective refugia of the Mediterranean peninsulas. Palaeogeogr. Palaeoclimatol. Palaeoecol. 625, 111699. https://doi.org/10.1016/j.paleo.2023.111699.

Maddison, W.P., Maddison, D.R., 2023. Mesquite: A Modular System for Evolutionary Analysis [WWW Document]. URL http://www.mesquiteproject.org (accessed 6.12.23).

Maia-Carvalho, B., Vale, C.G., Sequeira, F., Ferrand, N., Martínez-Solano, Í., Gonçalves, H., 2018. The roles of allopatric fragmentation and niche divergence in intraspecific lineage diversification in the common midwife toad (Alytes obstetricans). J. Biogeogr. 45, 2146-2158. https://doi.org/10.1111/jbi.13405.

Martin, J., 2014. Miocene Atlantic-Mediterranean Seaways in the Betic. Rev. La Soc. Geol. Espana 27, 175-186.

Martinez-Solano, I., 2004. Phylogeography of Iberian *Discoglossus* (Anura: Discoglossidae). J. Zool. Syst. Evol. Res. 42, 298-305. https://doi.org/10.1111/j.1439-0469.2004.00257.x.

Martinez-Solano, I., Teixeira, J., Buckley, D., García-París, M., 2006. Mitochondrial DNA phylogeography of *Lissotriton boscai* (Caudata, Salamandridae): evidence for old, multiple refugia in an Iberian endemic. Mol. Ecol. 15, 3375-3388. https://doi.org/10.1111/j.1365-294X.2006.03013.x.

McCartney-Melstad, E., Mount, G.G., Shaffer, H.B., 2016. Exon capture optimization in amphibians with large genomes. Mol. Ecol. Resour. 16, 1084-1094. https://doi.org/10.1111/1755-0998.12538.

McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A., Garimella, K., Altshuler, D., Gabriel, S., Daly, M., DePristo, M.A., 2010. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297-1303. https://doi.org/10.1101/gr.107524.110.20.

Médail, F., Myers, N., 2004. Mediterranean Basin. In: Mittermeier, R.A., Robles Gil, P., Hoffman, M., Pilgrim, J., Brooks, T., Mittermeier, C.G., Lamoreux, J., Da Fonseca, G. A.B. (Eds.), Hotspots Revisited: Earth's Biologically Richest and Most Endangered Terrestrial Ecoregions. Conservation International, Washington, pp. 144-147.

Médail, F., Quezel, P., 1999. Biodiversity hotspots in the Mediterranean Basin: setting global conservation priorities. Conserv. Biol. 13, 1510-1513.

Milá, B., Carranza, S., Guillaume, O., Clobert, J., 2010. Marked genetic structuring and extreme dispersal limitation in the Pyrenean brook newt *Calotriton asper* (Amphibia: Salamandridae) revealed by genome-wide AFLP but not mtDNA. Mol. Ecol. 19, 108-120. https://doi.org/10.1111/j.1365-294X.2009.04441.x.

Miller, M.A., Pfeiffer, W., Schwartz, T., 2010. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. 2010 Gatew. Comput. Environ. Work. GCE 2010. doi: 10.1109/GCE.2010.5676129.

Myers, N., Mittermeler, R.A., Mittermeler, C.G., Da Fonseca, G.A.B., Kent, J., 2000. Biodiversity hotspots for conservation priorities. Nature 403, 853-858. https://doi.org/10.1038/35002501.

Newman, C.E., Austin, C.C., 2016. Sequence capture and next-generation sequencing of ultraconserved elements in a large-genome salamander. Mol. Ecol. 25, 6162-6174. https://doi.org/10.1111/mec.13909.

Oosterbroek, P., Arntzen, J.W., 1992. Area-cladograms of Circum-Mediterranean taxa in relation to Mediterranean palaeogeography. J. Biogeogr. 19, 3. https://doi.org/10.2307/2845616.

Pabijan, M., Zielniški, P., Dudek, K., Chloupek, M., Sotiropoulos, K., Liana, M., Babik, W., 2015. The dissection of a Pleistocene refugium: phylogeography of the smooth newt, *Lissotriton vulgaris*, in the Balkans. J. Biogeogr. 42, 671-683. https://doi.org/10.1111/jbi.12449.

Peñalver-Alcázar, M., Jiménez-Valverde, A., Aragón, P., 2021. Niche differentiation between deeply divergent phylogenetic lineages of an endemic newt: implications for species distribution models. Zoology 144. https://doi.org/10.1016/j.zool.2020.125852.

Rambaut, A., Drummond, A.J., Xie, D., Baele, G., Suchard, M.A., 2018. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67, 901-904. https://doi.org/10.1093/sysbio/syy032.

Recuero, E., Canestrelli, D., Vörös, J., Szabó, K., Poyarkov, N.A., Arntzen, J.W., Crnobrnja-Isailović, J., Kidov, A.A., Cogálniceanu, D., Caputo, F.P., Nascetti, G., Martínez-Solano, I., 2012. Multilocus species tree analyses resolve the radiation of the widespread *Bufo bufo* species group (Anura, Bufonidae). Mol. Phylogenet. Evol. 62, 71-86. https://doi.org/10.1016/j.ympev.2011.09.008.

Recuero, E., García-París, M., 2011. Evolutionary history of *Lissotriton helveticus*: multilocus assessment of ancestral vs. recent colonization of the Iberian Peninsula. Mol. Phylogenet. Evol. 60, 170-182. https://doi.org/10.1016/j.ympev.2011.04.006.

Rowe, G., Harris, D.J., Beebee, T.J.C., 2006. Lusitania revisited: a phylogeographic analysis of the natterjack toad *Bufo calamita* across its entire biogeographical range. Mol. Phylogenet. Evol. 39, 335-346. https://doi.org/10.1016/j.ympev.2005.08.021.

Sambrook, J., Russell, D.W., 2001. Molecular Cloning: A Laboratory Manual, 3rd ed. Coldspring-Harbour Laboratory Press, London, United Kingdom.

Sánchez-Montes, G., Recuero, E., Barbosa, A.M., Martínez-Solano, I., 2019. Complementing the Pleistocene biogeography of European amphibians: testimony from a southern Atlantic species. J. Biogeogr. 46, 568-583. https://doi.org/10.1111/jbi.13515.

Sequeira, F., Bessa-Silva, A., Tarroso, P., Sousa-Neves, T., Vallinoto, M., Gonçalves, H., Martínez-Solano, I., 2020. Discordant patterns of introgression across a narrow hybrid zone between two cryptic lineages of an Iberian endemic newt. J. Evol. Biol. 33, 202-216. https://doi.org/10.1111/jeb.13562.

Sequeira, F., Arntzen, J.W., van Gulik, D., Hajema, S., Diaz, R.L., Wagt, M., van Riemsdijk, I., 2022. Genetic traces of hybrid zone movement across a fragmented habitat. J. Evol. Biol. 35, 400-412. https://doi.org/10.1111/jeb.13982.

Stamatakis, A., 2014. RAXML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312-1313. https://doi.org/10.1093/bioinformatics/btu033.

Steinfartz, S., Vicario, S., Arntzen, J.W., Caccone, A., 2007. A Bayesian approach on molecules and behavior: reconsidering phylogenetic and evolutionary patterns of the Salamandridae with emphasis on *Triturus* newts. J. Exp. Zool. Part B Mol. Dev. Evol. 308, 139-162. https://doi.org/10.1002/jez.b.21119.

Sun, C., Shepard, D.B., Chong, R.A., Arriaza, J.L., Hall, K., Castoe, T.A., Feschotte, C., Pollock, D.D., Mueller, R.L., 2012. LTR retrotransposons contribute to genomic gigantism in plethodontid salamanders. Genome Biol. Evol. 4, 168-183. https://doi.org/10.1093/gbe/evr139.

Teixeira, J., Martínez-Solano, I., Buckley, D., Tarroso, P., García-París, M., Ferrand, N., 2015. Genealogy of the nuclear *β*-fibrinogen intron 7 in *Lissotriton boscai* (Caudata, Salamandridae): concordance with mtDNA and implications for phylogeography and speciation. Contrib. to Zool. 84, 193-215. https://doi.org/10.1163/18759866-08403002.

van Riemsdijk, I., Arntzen, J.W., Babik, W., Bogaerts, S., Franzen, M., Kalaentzis, K., Litvinchuk, S.N., Olgun, K., Wijnands, J.W.P.M., Wielstra, B., 2022. Next-generation phylogeography of the banded newts (*Ommatotriton*): a phylogenetic hypothesis for three ancient species with geographically restricted interspecific gene flow and deep intraspecific genetic structure. Mol. Phylogenet. Evol. 167, 107361. https://doi.org/10.1016/j.ympev.2021.107361.

Vences, M., Hauswaldt, S.J., Steinfartz, S., Rupp, O., Goesmann, A., Künzel, S., Orozco-terWengel, P., Vieites, D.R., Nieto-Roman, S., Haas, S., Laugsch, C., Gehara, M., Bruchmann, S., Pabijan, M., Ludewig, A.-K., Rudert, D., Angelini, C., Borkin, L.J., Crochet, P.-A., Crottini, A., Dubois, A., Ficetola, G.F., Galán, P., Geniez, P., Hachtel, M., Jovanovic, O., Litvinchuk, S.N., Lymberakis, P., Ohler, A., Smirnov, N. A., 2013. Radically different phylogeographies and patterns of genetic variation in two European brown frogs, genus *Rana*. Mol. Phylogenet. Evol. 68, 657-670. https://doi.org/10.1016/j.ympev.2013.04.014.

Wickham, H., 2007. Reshaping Data with the reshape Package. J. Stat. Softw. 21, 1-20. https://doi.org/10.5325/soundings.101.1.0006.

Wickham, H., 2016. ggplot2: Elegant Graphics for Data Analysis, 2nd ed. Springer, Houston, Texas, USA. doi: 10.1007/978-3-319-24277-4.

Wielstra, B., Arntzen, J.W., 2016. Description of a new species of crested newt, previously subsumed in *Triturus ivanbureschi* (Amphibia: Caudata Salamandridae). Zootaxa 4109, 73-80. https://10.11646/zootaxa.4109.1.6.

Wielstra, B., Crnobrnja-Isailović, J., Litvinchuk, S.N., Reijnen, B.T., Skidmore, A.K., Sotiropoulos, K., Toxopeus, A.G., Tzankov, N., Vukov, T., Arntzen, J.W., 2013. Tracing glacial refugia of *Triturus* newts based on mitochondrial DNA phylogeography and species distribution modeling. Front. Zool. 10, 1-14. https://doi.org/10.1186/1742-9994-10-13.
Wielstra, B., Burke, T., Butlin, R.K., Avcı, A., Üzüm, N., Bozkurt, E., Olgun, K., Arntzen, J. W., 2017. A genomic footprint of hybrid zone movement in crested newts. Evol. Lett. 1, 93–101. https://doi.org/10.1002/evl3.9.

Wielstra, B., Canestrelli, D., Cvijanović, M., Denoël, M., Fijarczyk, A., Jablonski, D., Liana, M., Naumov, B., Olgun, K., Pabijan, M., Pezzarossa, A., Popgeorgiev, G., Salvi, D., Si, Y., Sillero, N., Sotiropoulos, K., Zieliński, P., Babik, W., 2018. The distributions of the six species constituting the smooth newt species complex (*Lissotriton vulgaris* sensu lato and *L. montandoni*) - An addition to the New Atlas of Amphibians and Reptiles of Europe. Amphib. Reptil. 39, 252–259. https://doi.org/10.1163/15685381-17000128.

Wielstra, B., McCartney-Melstad, E., Arntzen, J.W., Butlin, R.K., Shaffer, H.B., 2019. Phylogenomics of the adaptive radiation of *Triturus* newts supports gradual ecological niche expansion towards an incrementally aquatic lifestyle. Mol. Phylogenet. Evol. 133, 120–127. https://doi.org/10.1016/j.ympev.2018.12.032.

Zeisset, I., Beebee, T.J.C., 2008. Amphibian phylogeography: a model for understanding historical aspects of species distributions. Heredity 101, 109–119. https://doi.org/10.1038/hdy.2008.30.

Zheng, X., Levine, D., Shen, J., Gogarten, S.M., Laurie, C., Weir, B.S., 2012. A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics 28, 3326–3328. https://doi.org/10.1093/bioinformatics/bts606.
Wielstra, B., Burke, T., Butlin, R.K., Avcı, A., Üzüm, N., Bozkurt, E., Olgun, K., Arntzen, J. W., 2017. A genomic footprint of hybrid zone movement in crested newts. Evol. Lett. 1, 93–101. https://doi.org/10.1002/evl3.9.

Wielstra, B., Canestrelli, D., Cvijanović, M., Denoël, M., Fijarczyk, A., Jablonski, D., Liana, M., Naumov, B., Olgun, K., Pabijan, M., Pezzarossa, A., Popgeorgiev, G., Salvi, D., Si, Y., Sillero, N., Sotiropoulos, K., Zieliński, P., Babik, W., 2018. The distributions of the six species constituting the smooth newt species complex (*Lissotriton vulgaris* sensu lato and *L. montandoni*) - An addition to the New Atlas of Amphibians and Reptiles of Europe. Amphib. Reptil. 39, 252–259. https://doi.org/10.1163/15685381-17000128.

Wielstra, B., McCartney-Melstad, E., Arntzen, J.W., Butlin, R.K., Shaffer, H.B., 2019. Phylogenomics of the adaptive radiation of *Triturus* newts supports gradual ecological niche expansion towards an incrementally aquatic lifestyle. Mol. Phylogenet. Evol. 133, 120–127. https://doi.org/10.1016/j.ympev.2018.12.032.

Zeisset, I., Beebee, T.J.C., 2008. Amphibian phylogeography: a model for understanding historical aspects of species distributions. Heredity 101, 109–119. https://doi.org/10.1038/hdy.2008.30.

Zheng, X., Levine, D., Shen, J., Gogarten, S.M., Laurie, C., Weir, B.S., 2012. A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics 28, 3326–3328. https://doi.org/10.1093/bioinformatics/bts606.