An amphibian species pushed out of Britain by a moving hybrid zone

Abstract Classical theory states that hybrid zones will be stable in troughs of low population density where dispersal is hampered. Yet, evidence for moving hybrid zones is mounting. One possible reason that moving zones have been underappreciated is that they may drive themselves into oblivion and with just the superseding species remaining, morphological and genetic signals of past species replacement may be difficult to appreciate. Using genetic data (32 diagnostic single nucleotide polymorphisms) from a clinal hybrid zone of the common toad (Bufo bufo) and the spined toad (Bufo spinosus) in France for comparison, alleles of the latter species were documented in common toads in the south of Great Britain, at frequencies in excess of 10%. Because long distance dispersal across the Channel is unlikely, the conclusion reached was that the continental toad hybrid zone which previously extended into Britain, moved southwards and extirpated B. spinosus. Species distribution models for the mid‐Holocene and the present support that climate has locally changed in favour of B. bufo. The system bears resemblance with the demise of Homo neanderthalensis and the rise of Homo sapiens and provides an example that some paleoanthropologists demanded in support of a hominin “leaky replacement” scenario. The toad example is informative just because surviving pure B. spinosus and an extant slowly moving interspecific hybrid zone are available for comparison.


| INTRODUC TI ON
Hybrid zone studies have a prominent place in evolutionary research because they provide opportunities to study the mechanisms that keep species apart in the face of gene flow (Abbott et al., 2013;Ravinet et al., 2017). They can be seen as natural arenas where combinations of mutations from different species are functionally tested (Zieliński et al., 2019). Most recognized hybrid zones have come into existence upon the secondary contact of groups of organisms that are morphologically and genetically differentiated. Hybrid zones that are independent of the environment and in which the admixed offspring has lower fitness than the parents, are termed "tension zones". These are maintained by an equilibrium between dispersal into the zone and selection against hybrids (Barton, 1979b;Barton & Gale, 1993;Key, 1968). Classical hybrid zone theory states that tension zones will be stable, or at least persistent, in troughs of low population density where dispersal is reduced. These troughs would be frequent and deep enough to resist large-scale movement due to species fitness differences or asymmetry of hybridization (Barton, 1979b; Barton & Hewitt, 1985). The perspective of hybrid zone stability has, however, eroded under a flood of empirical data from three lines of evidence, namely (a) historical references and observations over time, (b) isolated populations surrounded by the superseding species (enclaves) and (c) alleles left behind by the receding species (genetic footprints) (Arntzen & Wallis, 1991;Buggs, 2007;Wielstra, 2019). One possible reason that moving hybrid zones have been underappreciated is that they may drive themselves into oblivion and with just the superseding species remaining, morphological and genetic signals of past species replacement may be difficult to appreciate. As Levin (2006) stated, finding the "ghost of hybridization past" can be elusive.
This study follows the discovery and description of an amphibian hybrid zone across France that formed following species' northerly range expansions, starting with the Holocene glacial retreat (Arntzen, McAtear, Butôt, & Martínez-Solano, 2018;Recuero et al., 2012;Trujillo, Gutiérrez-Rodríguez, Arntzen, & Martínez-Solano, 2017). The species involved are the common toad Bufo bufo to the northeast and the spined toad Bufo spinosus to the southwest of their mutual species border (Figure 1). Evidence is mounting for this hybrid zone to be moving along the Channel coast (Arntzen et al., 2016;. To test this hypothesis this study set out to investigate if-when the Channel had yet to be formed-the Bufo species contact might have extended into Great Britain.

| MATERIAL S AND ME THODS
Tissue material was obtained as toe tips from adults and the tips of tailfins from larvae for 145 individuals at 15 localities in Great Britain (Table 1) and as extracted DNAs from 358 individuals from Dorset, United Kingdom (Coles, Reading, & Jehle, 2019) and 318 individuals at 21 localities in the southwest of Norway (Roth & Jehle, 2016).
New material was collected in accordance with the UK's Wildlife and Countryside Act, 1981. Fluorescence-based genotyping (Semagn, Babu, Hearne, & Olsen, 2014) was used in the Kompetitive Allele-Specific PCR (KASP) genotyping system at the SNP (single nucleotide polymorphism) genotyping facility of the Institute of Biology, Leiden University. Primer design, PCR setup and data visualization followed Arntzen et al. (2016) and .
SNP's were selected on the basis of species diagnosticity. The data for Britain and Norway were analyzed jointly with published data that represent a Bufo spinosus to Bufo bufo hybrid zone in the northwest of France, along with reference populations from Spain, southern France and the Channel island Jersey (B. spinosus) and northern France and the Netherlands (B. bufo) (Arntzen et al., 2016;Arntzen, de Vries, Canestrelli, & Martínez-Solano, 2017;Arntzen, Wilkinson, Butôt, & Martínez-Solano, 2014;. The full panel of SNP's was studied with the exclusion of c10orf2 and klhl that proved not diagnostic, banp that showed a widely displaced cline and egflam that was strongly out of Hardy-Weinberg equilibrium . Data for individuals with more than nine SNPs missing were discarded. For the remainder, missing data amounted to 2.8% (N = 974) for 31 nuclear markers and 0.1% (N = 1) for the one mitochondrial marker.
A two-species distribution model was constructed by contrasting presence data for both species (Arntzen, Canestrelli, & Martínez-Solano, 2019). It describes the contiguous ranges of Bufo spinosus and B. bufo across France and into Italy from environmental parameters and showed a good fit to the underlying data of 404 genetically investigated toad populations (AUC = 0.97 ± 0.007). The model was re-estimated to only include climate variables (bio01-bio19) and then applied to the mid-Holocene climate reconstructions of WorldClim v. 1.4 (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005; the nine model codes are BC, CC, CN, HE, HG, IP, ME, MG and MR).
The parameters bio14 and bio15 were excluded from the analysis because of their regular high discrepancy between data sets (Varela, Lima-Ribeiro, & Terribile, 2015). In the absence of firm guidance of which climate reconstruction would be most appropriate to apply F I G U R E 1 Distribution of the spined toad, Bufo spinosus (red) and the common toad, Bufo bufo (blue) after Agasyan et al. (2009) and Arntzen et al. (2018). Study localities are shown by x-symbols. Inset: study localities in England. Note that Bufo toads are absent in Ireland and that their presence on the Isle of Man remains to be confirmed (NBNatlas, https ://speci es.nbnat las.org/speci es/NHMSY S0000 080159) The genetic data were analyzed with genepop version 4.2 (Rousset, 2008) and HIest (Fitzpatrick, 2012). Statistical data analysis was with spss version 20 (IBM SPSS, 2016). Habitat models were visualized with ilwis 3.6 (ILWIS, 2009). higher on the Isle of Wight (B01-06, 12.4% < F s < 15.8%) than in the south of mainland Britain (B07-10, 7.9% < F s < 11.2%) than to the north (B11-15 and N01-21, 0.4% < F s < 5.8%) ( differentiation between these three population groups is significant (Fisher's method, p < .001 for the three pairwise comparisons).

| RE SULTS
Allele frequencies for British and hybrid zone populations were to varying degrees correlated. The strongest signals were observed for the combination of southern British (B01-10) and northern hybrid zone populations (France, 12-16), with nine out of 15 pairwise comparisons statistically significant. Extremes in the marker spectrum are the loci pigg and pomc, with low respectively high frequencies of spinosus-alleles in a B. bufo genetic background in either country. The population from Norway has a relatively high frequency of spinosus-alleles (2.0% < F s < 4.8%), which is almost entirely due to the contribution of the loci aimp2 (F s = 70.5%) and med8 (F s = 25.1%).
The selected two-species distribution model for climatic data is represented by the logistic equation and not for the four others (0.55 < AUC <0.69).

| Genetic and climatic support for a past enclave
Hybrid zone movement may result in substantial unidirectional introgression of selectively neutral material from the local to the advancing species, leaving a genetic footprint (Scribner & Avise, 1993;Wielstra et al., 2017). In France, a weak but statistically significant Strengths of introgression may differ between genomic regions, depending on selection promoting or counteracting genetic exchange (Baack & Rieseberg, 2007;Roux, Tsagkogeorga, Bierne, & Galtier, 2013).
Regions that introgress more easily than the genomic background are likely to contain adaptive variants, whereas the flow of the genes involved in reproductive isolation and genes linked to them is impeded and yet others, perhaps the majority, may be selectively neutral (Barton, 1979a;Piálek & Barton, 1997). The populations in Norway are almost completely devoid of spinosus-alleles, except at the loci aimp2 and med8 ( Table 2). Considering that genetic variation usually decreases towards a species' range edge, the standing genetic variation observed in Norway is high (Roth & Jehle, 2016). A pattern where alleles reappear at localities away from the admixture zone has been observed in European house mice (Teeter et al., 2009) and is suggestive of differential negative selection, stronger within the area of admixture than away from it.
The rare alleles may, however, also reflect a different (i.e., northeastern) post-glacial colonization route for Scandinavia from a source carrying these variants, versus a more southern, homozygous origin for France and Great Britain. To resolve if the B. spinosus range extended into Scandinavia requires a wider phylogeographic study.
The frequencies of introgressed alleles in Britain and France are correlated across loci (Table 2). This result could be explained by the existence of a single, continuous hybrid zone for a long time, so that gene flow within species generated associations of allele frequencies between Britain and France. However, both sections of the hybrid zone became disconnected at c. 7.5 kya (see below) and ever since, British and French population groups followed independent evolutionary trajectories. Given the potential for differentiation, the observed correlations indicate that the introgression processes were similar in both countries, suggesting that a large proportion of the

| Historical biogeography
The Note: The correlation coefficient applied is that of Spearman (r s ), with values and significances shown as bold (*p < .05, **p < .01, ***p < .001 and no mark p > .05, not significant).
Abbreviations: A, average F s ; R, F s range; #, markers aimp2 and med8 excluded.
contact originated in the Holocene, through range expansion from glacial refugia in Spain and the south of France (B. spinosus) and the northern Balkans or beyond (B. bufo) . A relatively short period of climate warming during the Older Dryas (the Bøling-Allerød period beginning at 14.7 kya) provided a first opportunity for northward expansion in many vertebrate species (Montgomery, Provan, McCabe, & Yalden, 2014). Bufo toads may not have been amongst these early range expanders, or not have been dispersing fast enough to reach Ireland before it became disconnected from the European continent at c. 15 kya . Range expansions may have temporarily stagnated during the Younger Dryas (12.9-11.7 kya), the last cold event of the Pleistocene preceding the final warming at the beginning of the Holocene . However, both toad species did make it into Britain, traversing the now inundated "Doggerland" landmass. Under a warming climate and a rising sea level the former "Channel river" (Antoine et al., 2003) enlarged to become the Channel. Consequently, the later a species' expanding range hit upon the Channel coast, the more northerly the overland route to Britain was located. The extant but eroded enclave indicates that B. spinosus was the first species to arrive. A spatial argument underpinning this notion is that the distance from the nearest glacial refugium is shorter for B. spinosus (c. 900 km) than for B. bufo (c. 1,500 km) . Currently no data are available that compare locomotion speed, but given that they are morphologically similar species, differences may be small. On its way north B. spinosus reached Normandy and the Cotentin peninsula 8 kya the latest, because the species made it into Jersey, which was the last of the Channel Islands to be separated from the European continent (Johnston, 1981).
Entering Britain implies that B. spinosus expanded its range through the Low Countries, reaching as far north as the then adjoined Thames, Rhine and Meuse estuaries (Coles, 2000). For either species, the latest opportunity to make it into Britain was at 7.5 kya (Yalden, 1999), which also provides a minimum age for the B. spinosus -B. bufo species interaction. Given its current position the hybrid zone may have then travelled with an estimated speed of 40 m per year. Hybrid zone movement may have been slower if more time was available, or it may have been higher prior to reaching terrain that B. bufo finds more difficult to cross, such as the Collines de Normandie, a range of hills that appears to act as a barrier to further southward displacement (Arntzen et al., 2016).
Other anuran species native to the British Isles are the common frog, Rana temporaria, and the natterjack toad, Epidalia calamita, which may have had a glacial refugium in or near the southwest of Ireland, on what was then dry land (Forbes, 1846;Rowe, Harris, & Beebee, 2006;Teacher, Garner, & Nichols, 2009). A similar scenario for Bufo toads is unsupported because either species is absent from Ireland.

| Genetic footprints, and parallels with the demise of Homo neanderthalensis
The frequency of spinosus-alleles in Britain decreases from south to north, with most of the decrease (from c. 15% to 5%) taking place over a distance of 200 km, from the South of England to the Midlands.
Yet, data are currently insufficient to determine the extent of the past B. spinosus enclave in the south of Britain. This not just reflects limited spatial and genetic sampling, but also population genetic processes including drift, selection and neutral diffusion that can reduce, amplify or distort spatial genetic signals, as well as habitat dependent demographic parameters such as reproductive success and the propensity for dispersal (see Excoffier, Foll, & Petit, 2009;Harrison & Larson, 2014 for reviews (Currat & Excoffier, 2004).
An erratic, large-scale pattern appears to describe the spatial-genetic interactions of the lineages of anatomically modern humans and for example, Neanderthals, with frequent interbreeding when hominin species' ranges overlapped (Petr, Pääbo, Kelso, & Vernot, 2019;Slon et al., 2018;Vernot & Akey, 2014). The inferred species transition that was labelled as "leaky replacement" (Pääbo, 2015; see also Gibbons, 2011), parallels that of a moving hybrid zone leaving a genetic footprint. One difference may be that leaky replacement inside well-developed hybrid zones early generation hybrids may be rare, for example due to assortative mating, or to habitat segregation as in mosaic hybrid zones in which taxa map onto patches of interdigitated habitats. Third, the author seems to refer to stable hybrid zones, therewith ignoring that hybrid zones may be moving and ephemeral (Buggs, 2007;Wielstra, 2019). Fourth, uncovering traces of introgression, as in case of the H. sapiens -H. neanderthalensis interaction, requires extensive data that may be available in hominin paleogenomics (Fu et al., 2016;Racimo et al., 2015;Sankararaman, Mallick, Patterson, & Reich, 2016) but not in most non-model organisms. Finally, to brush the leaky replacement scenario aside as "human exceptionalism" (Varki, 2016- (Carrión et al., 2018;Jennings, Finlayson, Fa, & Finlayson, 2011). Norway (RJ, SR).

AUTH O R CO NTR I B UTI O N S
The project was designed and carried out by JWA, who also prepared the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
Newly obtained SNP-data are provided in Supplementary Material S1