Tracing species replacement in Iberian marbled newts

Abstract Secondary contact between closely related species can lead to the formation of hybrid zones, allowing for interspecific gene flow. Hybrid zone movement can take place if one of the species possesses a competitive advantage over the other, ultimately resulting in species replacement. Such hybrid zone displacement is predicted to leave a genomic footprint across the landscape in the form of asymmetric gene flow (or introgression) of selectively neutral alleles from the displaced to the advancing species. Hybrid zone movement has been suggested for marbled newts in the Iberian Peninsula, supported by asymmetric gene flow and a distribution relict (i.e., an enclave) of Triturus marmoratus in the range of T. pygmaeus. We developed a panel of nuclear and mitochondrial SNP markers to test for the presence of a T. marmoratus genomic footprint in the Lisbon peninsula, south of the enclave. We found no additional populations of T. marmoratus. Analysis with the software Structure showed no genetic traces of T. marmoratus in T. pygmaeus. A principal component analysis showed some variation within the local T. pygmaeus, but it is unclear if this represents introgression from T. marmoratus. The results may be explained by (a) species replacement without introgressive hybridization and (b) displacement with hybridization followed by the near‐complete erosion of the footprint by purifying selection. We predict that testing for a genomic footprint north of the reported enclave would confirm that species replacement in these marbled newts occurred with hybridization.

Enclaves form when the population of one species is surrounded by populations of a mutually exclusive, competing species, therewith becoming genetically isolated from the remainder of the receding species' range (Arntzen, 1978;Littlejohn & Roberts, 1975). Enclaves can therefore illustrate historical species replacement, particularly in ground-dwelling organisms with low dispersal capabilities, such as amphibians. The formation of enclaves may either take place with or without gene flow between the species involved.
The Eurasian newt genus Triturus has been previously employed to investigate moving hybrid zones and enclave formation, for nine species at different spatial, environmental and phylogenetic settings (Arntzen & Wallis, 1991;Arntzen et al., 2014;Wielstra et al., 2017).
The genus includes the northern marbled newt Triturus marmoratus (Latreille, 1880), which inhabits central and southern France as well as the northern part of the Iberian Peninsula, and the pygmy marbled newt T. pygmaeus (Wolterstorff, 1905), which occupies the southwestern parts of the Iberian Peninsula. These are sister species that engage in a hybrid zone spanning ca. 600 km across the west of the Iberian Peninsula (Figure 1a). We previously documented an enclave of T. marmoratus in the northwest of the Lisbon Peninsula, near the town Caldas da Rainha, composed of pure and introgressed populations ( Figure 1b). Another remarkable spatial signature is the presence of T. pygmaeus in and along the coastal dunes, approximately from Caldas da Rainha to the Aveiro Lagoon ( Figure 1). These observations suggest the northward competitive advance of T. pygmaeus over a stretch of more than 200 km, after which its range expansion was halted by the Rio Vouga estuary (Arntzen et al., in press).
The documented distribution pattern suggests that species replacement in Iberian marbled newts at the Portuguese Atlantic coast might take place through a moving contact zone, possibly leaving behind a tail of asymmetric introgression. We here test for the presence of genomic footprints of T. marmoratus in T. pygmaeus of the Lisbon Peninsula, south of the documented enclave, where we suspect that species turnover has taken place. For this study, we developed a panel of single-nucleotide polymorphism markers (SNPs; Garvin et al., 2010) that were species diagnostic in a large reference sample, available from an earlier study (Arntzen, 2018). In addition, F I G U R E 1 The Iberian Peninsula with records of Triturus marmoratus (solid round symbols) and T. pygmaeus (open round symbols). (a) Localities sampled for the evaluation of SNP marker diagnosticity. The boxed area includes the Lisbon Peninsula and the Caldas da Rainha area, with an inset of Europe for reference. The Aveiro Lagoon is shown by a four point star. (b) Localities sampled in the Lisbon Peninsula, with symbols as above. The continuous distribution with T. marmoratus in dark gray and T. pygmaeus in light gray is from Arntzen et al. (2009) and Arntzen (2018). Note the existence of a T. marmoratus enclave around Caldas da Rainha. The locality Juncal with both species and hybrids is shown by a cross. (c) Bar plot with individual admixture proportion (Structure Q scores) for T. marmoratus and T. pygmaeus reference populations and Lisbon Peninsula populations, with a detailed view of the Caldas da Rainha enclave (populations 1 and 2) we modeled the mutual species' distributions from climatic conditions of the present day, the Mid-Holocene and the Last Glacial Maximum, to investigate whether climate change would support the scenario of species replacement.

| Fieldwork, sampling, and DNA preparation
Fieldwork was carried out in 2018 and 2019 in the Lisbon Peninsula, where we searched for standing water bodies containing marbled newts. Twenty-five localities with newt populations were found in an area spanning ca. 2,000 km 2 , ranging from the Tejo river in the east, the Tejo estuary in the south, and the Atlantic Ocean in the west.
In the north, we sampled up to the town Caldas da Rainha where our survey adjoined the area previously investigated (Espregueira Themudo & Arntzen, 2007). We specifically included the Serra de Sintra mountains in the southwest of the Lisbon peninsula. This was because we suspected the (past) presence of T. marmoratus in these mountains, on account of (a) the distribution of T. marmoratus at higher altitudes than T. pygmaeus elsewhere along the species' mutual contact (Arntzen & Espregueira Themudo, 2008), and (b) a pilot study that revealed weak evidence for the presence of T. marmoratus genetic material at Sintra (Arntzen et al., in press).
Adult and larval newts were captured by dip netting or with funnel traps. To reduce sampling bias, for example, toward siblings from particular breeding pairs (Goldberg & Waits, 2010), we made an effort to include all accessible sections of the water bodies. Tail tip tissue samples were collected and stored in 96% ethanol. We also studied material from seven localities obtained earlier (Espregueira Themudo & Arntzen, 2007). DNA extraction of tissue samples followed the Kingfisher™ (Thermo Scientific) and DNeasy extraction kit (Qiagen) standard protocols.  (Bolger et al., 2014) under default settings, except for the parameter "minlen" that was set to exclude read lengths below 60 base pairs. Accordingly, the number of paired reads dropped from 60.34 10 6 to 46.21 10 6 for T. pygmaeus and from 56.04 10 6 to 42.67 10 6 for T. marmoratus. Trinity v2.5.1 (Grabherr et al., 2011;Haas et al., 2013) was employed for de novo transcriptome assembly under default settings, resulting in 19.6 10 4 contigs for T. pygmaeus and 18.8 10 4 contigs for T. marmoratus. The transcriptome assembly was blasted against itself in order to identify and remove all hits occurring more than once, therewith discarding all contigs possibly representing paralogs (Altschul et al., 1990).

| SNP marker design
Paralog filtering resulted in 12.3 10 4 contigs for T. pygmaeus and 11.2 10 4 contigs for T. marmoratus. An outline of the procedure is shown in Figure 2.
Single-nucleotide polymorphism marker design followed the molecular inversion probes pipeline that encompasses advantages for targeted resequencing, including high specificity, a high level of multiplexing and no ascertainment bias (Hardenbol et al., 2003;Niedzicka et al., 2016). Given that the T. pygmaeus transcriptome had a higher number of contigs, the following exon selection step was carried out using this dataset. Gene models were constructed based on the Xenopus tropicalis (Gray, 1864) frog genome, available from Biomart ENSEMBL (genome version JGI4.2). To further remove potential paralogs, T. pygmaeus contigs were mapped to the gene models to select exons that map to single X. tropicalis genes; these unique hits, along with the exon boundary coordinates, were used to extract exon sequences from T. pygmaeus contigs.  Figure 1. Localities with an asterisk have previously been studied by Espregueira Themudo and Arntzen (2007).

TA B L E 2
Localities for the 32 study populations of marbled newts on the Lisbon peninsula, with sample sizes (N) occurred over two times, twice or once within the X. tropicalis models (Altschul et al., 1990). The filtering script prioritized exons found more than twice in the gene model ("priority 1"), followed by exons present twice ("priority 2") and, finally, exons occurring only once ("priority 3"). Priority 1 exons were chosen as follows: whenever exons were present three times in the gene model, the exon located in the center was selected and if found over three times, an exon away from the extremes of the model was randomly selected. Priority 2 exons were designated by selecting the longest exon. Priority 3 exons were discriminated against as they provide no spatial information given their single presence within a gene model, with the risk that they would contain noncoding sequence.
The model strictly selected exons based on length, selecting exons with a minimum of 100 base pairs to ensure sufficient length for primer attachment and a maximum of 400 base pairs, considering average exon length (Sakharkar et al., 2004). This noncoding sequence removal step resulted in 3,372 T. pygmaeus exons of which 42 exons were priority 1, 228 were priority 2 and 3,102 were priority 3.
Triturus pygmaeus exons were blasted against the T. marmoratus transcriptome contig database, with potential paralogs being manually removed from the BLAST output (Altschul et al., 1990). after performing a sequence alignment with Muscle v3.8.31 (Edgar, 2004).

| SNP detection and validation
Single-nucleotide polymorphism genotyping took place at the For SNP validation, we used 120 T. pygmaeus and 118 T. marmoratus from 43 reference populations across the Iberian Peninsula, located outside the documented hybrid zone of these species (Arntzen, 2018; Figure 1,  for the nuclear genes β-Fibrinogen intron 7 (BF), Calreticulin intron C (CC), and the Platelet-derived growth factor receptor α intron 11 (PDG). Single-nucleotide polymorphisms were considered species diagnostic for 54 nuclear markers with a low Cohen's kappa (1-ĸ) > 0.9.

| Population genetics
Hardy-Weinberg equilibrium and genotypic disequilibrium among pairs of nuclear loci were tested with the R package GENEPOP v1.0.5, under the Benjamini-Hochberg correction (Benjamini & Hochberg, 1995;Rousset, 2008). Gene flow between genetically distinct populations produces admixture linkage disequilibrium among those loci that have different allele frequencies in the founding populations (Pfaff et al., 2001). Admixture linkage disequilibrium was estimated following Barton and Gale (1993) and was based on The SNP data were analyzed with Structure software under the "admixture model," given that neighboring populations can interbreed (Falush et al., 2003;Pritchard et al., 2000). The program was run for 100,000 generations following 100,000 generations of burn-in, with 10 replicates. The number of potentially differentiated gene pools was analyzed over the 1 < K < 10 range. The best K was selected using the "Evanno-method" with StructureHarvester (Earl & vonHoldt, 2012

| Environmental modeling
A two-species distribution model was constructed by the comparison of presence-only data for both species under reference to contemporary climate conditions. The biological data employed were 108 T. marmoratus and T. pygmaeus records that were supported by molecular species identification (Arntzen, 2018;present paper (Clarke & Gorley, 2006).
Variables were retained using the criterion of partial independence at r s < .7 ( Figure 3). Accordingly, variables selected for analysis were

| RE SULTS
Single-nucleotide polymorphisms were considered species diagnostic for 55 out of 64 nuclear and mitochondrial markers (  (Table S3).
In the Structure analysis, the number of genetically differentiated groups supported by the data was K = 2. Minor support was found for K = 3, and no support was obtained for Structure K-values higher than that (Figure 4). Variables selected for analysis using the criterion of partial independence at r s < .7 were bio01, bio02, bio03, bio05, bio06, bio12, and bio17 The selected two-species distribution model is represented The spatial interpretation of the model is shown in Figure 6a.
Temporal extrapolations of the two-species distribution models (or F I G U R E 4 Estimation of the best number of K from an assumed range of 2-9 based on the Evanno method. ∆K was calculated as mean(|L″(K)|/sd(L(K))

F I G U R E 5 Principal component analysis of genetic variation at 54 biallelic nuclear loci in marbled newts from the Iberian Peninsula. (a)
Results for Triturus marmoratus (gray shading) and T. pygmaeus (colors) for reference and Lisbon peninsula populations are summarized by ellipses that represent the mean ± one standard deviation. The amount of the total variation explained along the first and second PCA-axis is 83.8% and 3.1%, respectively. (b) Contour plot of PCA2 scores for T. pygmaeus at the Lisbon peninsula. Results are shown for reference and Lisbon peninsula populations, with a per locality radius of ca. 100 km; see PCA2 color legend to the left. Note that Lisbon peninsula populations from the Serra de Sintra and adjacent to the T. marmoratus enclave have high scores (PCA2 > 1.5) and that centrally located populations have intermediate scores. Most reference populations have low scores (PCA2 < 1.5) "hindcasts") for the Mid-Holocene are roughly similar to the one for the present (Figure 6a,b) whereas one out of three models for the Last Glacial Maximum support a more southerly located mutual species border (Figure 6a,c).

| D ISCUSS I ON
We employed a panel of presumably unlinked neutral markers to test for genetic traces of species replacement with hybridization in two species of marbled newts in Portugal. A northward hybrid zone movement has been proposed for this system, with a documented  Portugal are different from those more to the northwest. It must be noted that our pilot study also showed substantial variation within T. pygmaeus (Arntzen et al., in press).
We identified the parameter "amount of precipitation in the driest quarter of the year" (bio17) as most closely and significantly associated with the two-species distributions ( Figure 6). Newts of  (Kupfer, 1998;Kupfer & Kneitz, 2000;Trochet et al., 2014), the inferred 200 km range extension could be achieved in a couple of centuries, even when competition with a related species slows down the advance (Arntzen & Wallis, 1991;Wielstra et al., 2017). The reconstructed paleo-climatic data yield little support for a wider, more southerly distribution of T. marmoratus in the Mid-Holocene ( Figure 6b) and for the Last Glacial Maximum the support is limited to one climate reconstruction out of the three that are available (Figure 6c). The advancement of  (Arntzen et al., 2009;Coughlan & Matute, 2020;Zuiderwijk, 1990). Hybridization between T. marmoratus and T. pygmaeus is not frequent, as is illustrated by the narrow characteristics of the hybrid zone for many characters (Arntzen, 2018), as well as by the local population Juncal (Figure 1), where species composition is bimodal (two hybrids were found along with twelve T. marmoratus and two T. pygmaeus; Espregueira Themudo & Arntzen, 2007).
Variants of this scenario are the retraction of T. marmoratus followed by expansion of parapatric T. pygmaeus into abandoned ponds, or that hybridization and introgression are too infrequent to be detected with the current biological and genetic sampling schemes.
Transcriptomic-derived markers may not be selectively neutral and possibly fail to detect asymmetric introgression; therefore, we performed tests to identify any markers under selection and concluded our markers displayed neutral behavior. An argument against this explanation is that T. marmoratus and T. pygmaeus do, at the present day, hybridize along the species contact zone (Arntzen, 2018;Arntzen et al., in press).
Secondly, species replacement may have occurred with hybridization. The genetic differentiation of T. pygmaeus in the Lisbon peninsula could be attributed to hybridization and purifying selection, with introgressed alleles removed at some loci and not others. Strong selection against hybrids (purifying selection), as suggested by the low numbers of admixed individuals reported so far, could have eroded the footprint. In support of this scenario is that simulation models suggest that narrow (bimodal) hybrid zones-as observed for T. marmoratus and T. pygmaeus-are consistent with strong selection against hybrids (Currat et al., 2008;Singhal & Moritz, 2011).
If the force of selection is sufficient, only sporadical take-up of advantageous alleles by the invading species may take place, as has been previously shown for toads (Bufo) in the UK (Arntzen, 2019).
Such a scenario would be comparable to the phenomenon of "leaky replacement" of Neanderthals by other hominins (Hedrick, 2013;Pardo-Diaz et al., 2012;Racimo et al., 2015). warming (Thomas et al., 2006;van de Vliet et al., 2014). This would parallel observations on another Triturus moving hybrid zone that was shown to have halted due to the loss of breeding sites (Visser et al., 2017). Understanding shifts in species distributions, particularly when driven by climate change and anthropogenic activities, are especially relevant when deciphering the dynamics of species replacement (Taylor et al., 2015).