Genetic traces of hybrid zone movement across a fragmented habitat

Abstract

Theoretical and empirical studies suggest that the structure and position of hybrid zones can change over time. Evidence for moving hybrid zones has been directly inferred by repeated sampling over time, or indirectly through the detection of genetic footprints left by the receding species and the resulting asymmetric patterns of introgression across markers. We here investigate a hybrid zone formed by two subspecies of the Iberian golden‐striped salamander, Chioglossa lusitanica, using a panel of 35 nuclear loci (31 SNPs and 4 allozymes) and one mitochondrial locus in a transect in central Portugal. We found concordant and coincident clines for most of the nuclear loci (n = 22, 63%), defining a narrow hybrid zone of ca. 6 km wide, with the centre positioned ca. 15 km south of the Mondego River. Asymmetric introgression was observed at another 14 loci. Their clines are displaced towards the north, with positions located either close to the Mondego River (n = 6) or further northwards (n = 8). We interpret these profiles as genetic traces of the southward displacement of C. lusitanica lusitanica by C. l. longipes over the wider Mondego River valley. We noted the absence of significant linkage disequilibrium, and we inferred low levels of effective selection per locus against hybrids, suggesting that introgression in the area of species replacement occurred under a neutral diffusion process. A species distribution model suggests that the C. lusitanica hybrid zone coincides with a narrow corridor of fragmented habitat. From the position of the displaced clines, we infer that patches of locally suitable habitat trapped some genetic variants that became disassociated from the southward moving hybrid zone. This study highlights the influence of habitat availability on hybrid zone movement.

Abstract

Admixture proportions are tri‐modal in a hybrid zone of salamander subspecies, in line with habitat fragmentation along the investigated transect.

Open in new tabDownload slide

INTRODUCTION

The study of hybridization is an important avenue of research to address fundamental questions in the field of evolutionary biology, in particular towards an improved understanding of the speciation process (Abbott et al., 2013). Ever since Hewitt's studies (Hewitt, 1988, 1996, 2000), it has become widely accepted that Pliocene and Pleistocene climatic changes profoundly affected the population structure of many extant groups of organisms, through multiple cycles of range contraction, fragmentation and subsequent range expansion. In many cases, populations evolved independently while occupying distinct glacial refugia, leading to the formation of secondary contact zones after post‐glacial range expansions. When hybridization results in reproductively successful offspring, parental populations may eventually merge. Conversely, when hybridization results in individuals with low fitness, selection may favour gene variants that strengthen the barrier to gene flow, therewith entering the route towards complete reproductive isolation (Barton & Hewitt, 1985; Pereira & Wake, 2009; Seehausen, 2006; Wu, 2001). However, a frequently observed outcome is the persistence of a more or less stable area of admixture in contact (hybrid) zones among taxa at varying levels of divergence (Abbott et al., 2013; Barton & Hewitt, 1985; Hewitt, 1988; Moore, 1977).

Hybrid zones in which admixed offspring are less fit than the parents, independent of the geographic position and local ecological gradients, are so called ‘tension zones’. These are maintained by a balance between selection against hybrid genotypes and dispersal (Barton, 1979; Barton & Gale, 1993; Key, 1968), and tend to stabilize in areas of low population density or at an ecological barrier to dispersal (Barton & Hewitt, 1985; Endler, 1977). Because selection affecting the structure of hybrid zones is independent of the environment, tension zones are also able to move across the landscape in response to gradients in population density, fitness differences and varying selection pressures between hybridizing taxa (Barton, 1979; Barton & Hewitt, 1985; Bazykin, 1969; Endler, 1977; Key, 1968). Analysis of allele frequency changes across the landscape from temporally spaced samples over a few decades, has provided unequivocal evidence of contemporary spatial shifting (i.e. moving) hybrid zones, particularly in response to recent altering of the landscape and climate change (Buggs, 2007; Dasmahapatra et al., 2002; Leaché et al., 2017; Ryan et al., 2018; Taylor et al., 2015). Because repeated sampling over long timescales is complicated or impossible in hybrid zones for which spatial shifts have elapsed over time frames of hundreds or thousands of years, most moving hybrid zones have been inferred based on the detection of a geographical asymmetry of allelic introgression (for a review see Wielstra, 2019).

Theoretical and empirical evidence supports the view that the genomes of interbreeding species are heterogeneous regarding their permeability to foreign alleles (Endler, 1977; Gompert et al., 2017; Wu, 2001). Hybridizing taxa are expected to freely exchange neutral or advantageous alleles, whereas alleles at genomic regions that contribute to reproductive isolation through fitness variation in hybrids are not expected to be easily exchanged across taxon boundaries. In the last case, allele frequency transitions (that represent the centre of the cline) are expected to share the same geographic position (Barton, 1979; Barton & Bengtsson, 1986; Barton & Hewitt, 1985). In the wake of a moving hybrid zone, a ‘tail of introgression’ or ‘genetic footprint’ corresponds to unlinked, selectively neutral alleles derived from the displaced species and incorporated by populations of the advancing species (Currat et al., 2008; Excoffier et al., 2009; Scribner & Avise, 1993). This tail of introgression is predicted to persist because it is only subject to drift (Barton & Hewitt, 1985; Currat et al., 2008). However, different cline spatial configurations may arise within the same hybrid zone, because the amount of introgression may differ for genomic regions and depend on the strength of barriers to genetic exchange across the species boundary as well as on species spatial dynamics (e.g. Gompert et al., 2017; Quilodrán et al., 2020). The study of dynamic hybrid zones provides the opportunity to deepen our understanding about the role of various evolutionary forces in generating patterns of introgression and help to understand how changes in the interactions between hybridizing species can impact their genetic cohesion (Abbott et al., 2013; Barton & Hewitt, 1985; Ryan et al., 2018; Taylor et al., 2015).

The golden‐striped salamander, Chioglossa lusitanica, is an endemic species restricted to the north‐western corner of the Iberian Peninsula. Two morphologically cryptic subspecies established a secondary contact zone in central Portugal, following population expansions from separate local glacial refugia (Alexandrino et al., 2000, 2007). This mutual range border is a promising system to test for asymmetric introgression following hybrid zone movement. A small panel of diagnostic genetic markers showed that allele frequencies changed in parallel, except for mtDNA and one nuclear locus (out of four studied) at which allele frequency transitions were situated more to the north, suggesting southward hybrid zone movement (Sequeira et al., 2005). We here expand upon this work by studying 31 gene coding markers combined with distribution modelling of C. lusitanica. We aim to test the hypothesis of genome‐wide traces of the receding southern subspecies (C. lusitanica lusitanica) in the advancing northern subspecies (C. l. longipes), and we explore whether ecological/environmental factors across the hybrid zone may influence spatial patterns of introgression.

MATERIALS AND METHODS

Tissue sampling, transect design and data availability

We collected tail‐tip tissue samples from adult and larval Chioglossa lusitanica from one Spanish and 14 Portuguese populations, following an approximately south to the north sampling scheme (Figure 1, Table 1). Because the C. lusitanica range is constricted in central Portugal and the distribution of the southern subspecies is relatively small, it was not possible to achieve a long and linear transect. Therefore, we follow Sequeira et al. (2005) and Arntzen (2006), who identified a narrow latitudinal corridor across the Mondego River just east of Coimbra city. Accordingly, distances for localities close to the Mondego River were measured as line‐of‐sight distances to that river, with negative values for southern (localities 4–7) and positive values for northern populations (localities 8–10). For populations further away (the southern localities 1–3 and the northern localities 11–15), the transect runs through this corridor (Figure 1). Populations studied by Alexandrino et al. (2000, 2002) from north of the Mondego (localities A1‐A7) and Sequeira et al. (2005) from south of the Mondego (localities S1–S4), with available mitochondrial DNA (Cytochrome b, cytb) and enzyme data at four loci (Adh1, PepC, PepD and Pgm1), were also considered.

SNP marker development from transcriptome data

Two transcriptomes from liver tissue of one adult male of C. l. lusitanica from Foz do Alge (locality 1, NCBI code: SAMN25275766) and one adult male of C. l. longipes from Valongo (locality 14, EU880308, Rancilhac et al., 2021) were sequenced commercially by ZF Screens, Leiden, on the Illumina HiSeq 2000 platform. The transcriptome reads were filtered with Trimmomatic v.0.36 (Bolger et al., 2014) and assembled with Trinity v.2.6.5 (Grabherr et al., 2011) using default parameters. Paralog gene copies were filtered using a BLAST run (Altschul et al., 1990) of the data for C. l. lusitanica onto itself. A pipeline was used to identify exon boundaries (Niedzicka et al., 2016) using the Xenopus (Silurana) tropicalis genome assembly JGI_4.2. This genome was also used to annotate markers; after testing for diagnostic differences between the two subspecies, just six markers used were successfully annotated, Table S9. A second blast of the selected exons against the transcriptome of C. l. longipes was used to exclude potentially undetected paralogs and to identify SNPs. A detailed description, including custom python and bash scripts, is provided in SNPs pipeline user manual in Appendix S1. Figure S1 shows the steps taken, along with the amount of data retained in each step of the pipeline.

DNA from the new material was extracted using the EasySpin^® Genomic DNA Tissue Kit (Citomed). All individuals were genotyped with the Kompetitive Allele‐Specific PCR (KASP) system, based on fluorescence‐based genotyping (Semagn et al., 2014) at the SNP genotyping facility of the Institute of Biology at Leiden University. SNP primer design, PCR setup and data visualization followed Arntzen et al. (2016). For 88 exon sequences, primers were designed with the Kraken software (LGC genomics, UK). After a first run, 18 SNPs were excluded as they failed for all individuals of one subspecies. From the retained 70 SNPs, 31 were selected for which the SNPs differentiated reference samples of C. l. longipes (n = 42 from localities 1 and 2) and C. l. lusitanica (n = 42 from localities 13 and 14) with a Cohen's kappa (κ, Cohen, 1960) score of (1−κ) ≥0.9, indicating a high diagnosticity of these markers. Data for individuals with more than nine SNPs missing were discarded. For the remainder, missing data amounted to 5.6% (Tables S1 and S2). A total of eight low frequency alleles at the enzyme coding loci listed above (Alexandrino et al., 2000, 2002; Sequeira et al., 2005) could not unequivocally be attributed to one subspecies or the other and were ignored, and frequencies of the remaining two alleles rescaled from zero to one.

Population genetic analyses

Pairwise linkage disequilibrium and Hardy–Weinberg proportions were tested with the R version of GenePop (Rousset, 2008). We used the Bonferroni correction to account for multiple comparisons (N = 35), because accounting for all tests may be overly conservative (e.g. Narum, 2006; Rice, 1989).

Equilibrium cline models

The R package HZAR (Derryberry et al., 2014) was used to estimate allele frequency clines. The shape of each cline is described as a sigmoidal curve with a centre (maximum slope) and two tails (one on each side) modelled as an exponential decay, following Szymura and Barton (1986, 1991). Five different models were tested, varying in the presence versus absence of left and right tails and their symmetry versus asymmetry. Minimum (p  _min) and maximum (p  _max) character frequencies were used to scale cline models in three different ways: with empirical estimates (using lower and higher observed frequencies in the data), with best‐fit values and with no scaling (i.e. p  _min = 0 and p  _max = 1). The combination of cline shapes and scaling thus produces a total of 15 different models for each cline, plus a null model (as in Prada & Hellberg, 2013). The Metropolis–Hasting algorithm in HZAR was run twice for each locus using the default value of 100,000 steps with a randomly selected seed and a burn‐in of 10%. Given our occasionally uneven sample sizes, it is important to note that HZAR cline estimates are made with sample size taken into account. Convergence and stability of the parameters of the models selected as ‘best’ were assessed by visual inspection of the MCMC trace files. After each run, we recorded the corrected Akaike information criterion (AICc) score for each model within each cline and chose the one with the lowest score to infer cline width, cline centre and other cline shape parameters. We compared the log L values of the best selected models across runs and therewith confirmed convergence towards the best model.

Admixture proportions (usually referred to as Structure Q‐scores) inferred by the software Structure (Pritchard et al., 2000) have often been used to fit a geographic consensus cline, from which both width and centre of a hybrid zone are estimated (Arntzen et al., 2016; Dufresnes et al., 2019; Sequeira et al., 2020; Wielstra, Burke, Butlin, Avcı, et al., 2017). However, following a recent study (Toyama et al., 2020), the ability of structure to accurately estimate admixture proportions is affected by several factors (e.g. sampling scheme, degree of differentiation between parental populations, age of the admixture event and number of markers) that may result in artificially steeper clines and narrower hybrid zones. Accordingly, we opted for estimating the centre of the hybrid zone by determining the geographic position (c′) where the cumulative admixture is highest. This cumulative curve was derived from the bell‐shaped curves that were built for each locus separately, by using the first derivative of the corresponding sigmoid curve that was estimated by HZAR.

To examine the coincidence of cline centres, we refitted the best‐fitting cline model for each locus, but fixed the value at position c′. The fits of the two nested models (fixed centre or simpler model versus unconstrained model or full model) were compared using a likelihood ratio test, under the assumption that twice the difference in log likelihood (G = −2∆LL) of two models asymptotically follows a χ  ² distribution, with the degrees of freedom equal to the difference in the number of parameters. Finally, we investigated whether clines were displaced or coincident by comparing the two log‐likelihood unit support limits interpreted as the 95% confidence interval (hereafter CI) of the cline centre for each locus with the CI of c′. Clines are considered displaced if their CI’s do not overlap with the CI of c′.

For every cline, ‘surfaces under and over the curve’ (SUOC) were estimated with a custom R script, measuring the area underneath the cline curve north of the hybrid zone centre, up to the point where the frequency of the marker reached a value of 0.95 (SUOC north), or the area above the curve south of the hybrid zone centre until it reached a value of 0.05 (SUOC south). Their asymmetry (⁠$a^$⁠) was then measured by calculating $a^$ = log (SUOC south/SUOC north) over the F  _south (frequency of alleles in the southern taxon) zero to unity range. If this was not feasible, the range was set to 0.1 < F  _south < 0.9.

Selection against hybrids

Clines that are positioned in the same location suggest the existence of a barrier to introgression. Following Barton and Gale (1993), the strength of the genetic barrier was estimated from admixture linkage disequilibrium (D′), selection against hybrids (s*) and cline width (w) under neutrality, with a published R script (van Riemsdijk et al., 2019). Input parameters used were a generation time of four years (Lima et al., 2001), an assumed recombination rate ^® of 0.5 for unlinked loci (Szymura & Barton, 1986), and periods of 10,000 or 40,000 years corresponding to possible times of secondary contact since the last and the before‐last glacial period (Sequeira et al., 2005). To calculate the width of a neutral cline and effective selection against hybrids, we first estimated dispersal distance (standard deviation of distance between parent and offspring per generation) using linkage disequilibrium (D′) calculated from our genomic data following Barton and Gale (1993). Dispersal distance (σ) was then estimated using the formula: σ = √rD′w  ². We then combined the estimates of dispersal and time since secondary contact to calculate the width of a cline under neutrality, as w = 2.51σ √t. Finally, using our estimate of dispersal distance, we calculated the effective selection against hybrids using the formula s = 8σ  ²/w  ² (Barton, 1979; Bazykin, 1969). In these calculations, population 14 was excluded because it contains heterozygote genotypes, whereas in population 13, there are none (Tables S2 and S4). The heterozygosity in population 14 is thus to be considered as potential remnants of hybrid zone movement, but is not indicative of the selection against hybrids we intend to calculate.

Species distribution model

We modelled the distribution of C. lusitanica over the Iberian Peninsula. Species records were those of the Portuguese and Spanish atlases, organized at the Universal Transverse Mercator (UTM) grid system, which altogether occupied 356 different grid cells (Sequeira & Alexandrino, 2008; Vences, 2002). Considered environmental data layers were altitude in m a.s.l.; mean annual number of frost days—minimum temperature ≤0.5 C; annual range in the normalized difference vegetation index (NDVI); NDVI March 1993 minus NDVI December 1992; mean annual number of days with precipitation ≥0.1 mm; mean relative air humidity in January at 07:00 h in %; mean relative air humidity in July at 07:00—hours (%); maximum precipitation in 24 h in mm; monthly vegetation index composite April 1992; idem, August 1992; idem, December 1992; permeability of the soil in three categories of increasing permeability; mean annual precipitation in mm; relative maximum precipitation; the running vegetation lifeforms classification; slope in degrees; mean annual solar radiation in kWh/m²/day; mean annual temperature in °C; mean temperature in July in °C; and annual temperature range in °C. For details and the origin of the environmental data, see table 1 in Arntzen and Espregueira Themudo (2008). The model for C. lusitanica was constructed with MaxEnt software (Phillips et al., 2006) under default settings.

RESULTS

Significant heterozygote deficits were found in populations at the core of the transect (localities 5–7 and 9) and not at the others (1–4, 8 and 10–14), and not repeatedly at the same locus (Table S4). No significant values of admixture linkage disequilibrium (D′) were found (Figure S2).

Model selection from a set of nested cline models generally converged well, as shown by near‐equivalent LogL values for 34 markers out of 36 (Table S10). Maximum‐likelihood estimates of the cline centre (c) were highly variable among loci and are positioned between −17 km and 30 km (Figure 2, Tables 2 and S5), while c′ is located at −15.0 km, that is well south of the Mondego River and at the foothills of the Lousã mountains (around localities 5 and 6; Figure 3). Likelihood ratio tests based on the logit–logistic model indicated that 14 out of 31 (45%) of SNP clines and three out of four allozymes clines (75%) could comfortably be constrained to c′, and can thus be considered coincident. However, 17 out of 31 (55%) of the SNP clines, one enzyme cline out of four and the mtDNA marker showed a significant drop in the likelihood when constrained like this (Table S6). Among these non‐coincident clines, several correspond to small shifts of less than one cline width (cox18, exon54 and galk, see Tables 2 and S5). The other 14 SNP loci showed nonoverlapping CIs when compared to the CI of c′ (−16.2 < CI < −13.4 km) (Figure 2; Tables S5 and S6). These 14 displaced clines were all shifted to the north, with different levels in the height of the tail (Figure 2 and Table S5). Eight of these loci, and also mtDNA, showed a central position north of the Mondego River (12.5 < c < 30.0 km, north of population 10), whereas the other subset of six loci had cline centre estimates that clustered around the Mondego River (populations 7–10), or half‐way between the Mondego River and the Lousã mountains (−10.5 < c < 2.8 km) (Table 2; Figures 3 and S4). The result of 14 clines displaced to the north, and none to the south is significantly different from equal expectations (binomial test, p < 0.0001). Moreover, the coinciding 22 nuclear loci clines are on average significantly skewed to the north (⁠$a^$ = 0.28 ± 0.66 and different from zero, one‐sample t‐test, d.f. = 21, t = 1.98, p < 0.05), whereas the displaced clines are not significantly skewed (⁠$a^$ = 0.06 ± 0.25, d.f. = 13, t = 0.93, p > 0.05) (Table S7).

The estimates of cline width (w) were highly variable among loci, ranging from 1.84 km to 76.1 km (Tables 2 and S5). A significant positive relationship was observed between position of the cline and cline width (Spearman's correlation coefficient r  _s = 0.75, N = 35, p < 0.0001). Overall, the non‐coincident loci were outliers, with widths that ranged from 16.8 to 76.1 km (average 29.3 km; 18.8 < CI < −50.6 km). An exception is the mtDNA cline width (w = 7.0 km; 3.7 < CI < 14.6 km) that is within the range of the average widths for coincident clines (w = 6.1 km; 2.8 < CI < 12.6 Km) (Tables 2, S5 and S6).

Admixture linkage disequilibrium (D′) in the zone centre was 0.005 (0.0002 < CI < 0.0104). Based on this estimate, and with an assumed 10 Kya since secondary contact, the lifetime dispersal (σ) was estimated at 0.45 km generation^−1/2 (0.10 < CI < 0.84 km generation^−1/2, Table S8). The expected cline width assuming no selection against hybrids was 56.4 km (12.5 < CI < 105.3 km), which is substantially wider than the width estimates for most of the clines. The average effective selection on a locus (s*) was 0.0063 (0.0003 < CI < 0.0145). The analysis with 40 Kya since secondary contact did not result in markedly different results (Table S8).

The distribution model projects the presence C. lusitanica across the north‐western corner of the Iberian Peninsula (Figure 3a). The single variable ‘precipitation’ contributed more to the model (63.7%) than the other variables taken together (36.3%). The full model has an AUC value of 0.964, and the AUC value for a model with just precipitation is 0.922. The full model almost completely coincides with the distribution records of the species for Portugal and extends somewhat over the records documented for Spain (results not shown). The MaxEnt response curve indicates that areas with a precipitation of >1200 mm/year are favourable, and areas with a precipitation <1000 mm/year are unfavourable to the species (Figure 3b). Precipitation levels in the contact zone region are shown in Figure 3c, with a colour code that highlights variation over the critical 1000–1200 mm/year range.

DISCUSSION

The analysis of the Chioglossa hybrid zone extended to encompass the larger panel of 35 diagnostic nuclear loci reveals a sharp genetic transition between the southern and northern subspecies, with coincident and concordant clines for most loci. Other nuclear clines are wider and significantly displaced to the north. The cline for mtDNA is also displaced north, yet narrow. Mitochondrial DNA usually forms narrower clines, due to a lower effective population size compared with nuclear markers (Excoffier et al., 2009; Polechová & Barton, 2011; Toews & Brelsford, 2012).

The sharp transition at the centre of the hybrid zone suggests that many genomic regions locally experience a barrier to gene flow. However, cline shape alone does not allow alternative mechanisms underlying the maintenance of sets of narrow clines to be distinguished (Kruuk et al., 1999). One potential cause of steep, narrow and coincident clines is environment‐dependent selection (Endler, 1977). The differential ecological–climatic modelling of the two Chioglossa subspecies suggested that they occupy habitats that represent different environmental conditions (Arntzen & Alexandrino, 2004). However, these differences are unlikely to play a prominent role in the hybrid zone, because they are subtle and neither coincident with the Mondego River, nor with the subspecies range borders as represented by the set of coincident clines adjacent to the Lousã mountains. Another potential cause for the maintenance of a set of narrow clines is strong selection against hybrids (Barton & Hewitt, 1985). However, average effective selection against Chioglossa hybrids is as low as in other hybrid zones where introgression operates under a neutral diffusion process (Baldassarre et al., 2014; Wielstra, Burke, Butlin, Avcı, et al., 2017). Strong selection is also incompatible with the observed absence of significant linkage disequilibrium. One crucial parameter particularly difficult to assess is lifetime dispersal. Although the genetic estimate of σ = 0.45 km gen^−1/2 (0.10 < CI < 0.84, Table S8) is similar to that observed for salamanders in general (Smith & Green, 2005), it is not supported by field data on Chioglossa. With documented travel distances of 700 m in adults (Arntzen, 1981, 1984) and the species’ propensity for dispersal by larval drift (Arntzen, 1995; Thiesmeier, 1994), actual dispersal is likely to exceed the genetic estimate. On the contrary, displacements of this species outside moist habitats are severely restricted (Arntzen, 1995), so that dispersal over wider areas will at best be indirect, probably similar to the situation in other stream‐associated salamanders (Grant et al., 2010; Miller et al., 2015) and riverine amphibian species (Vörös et al., 2016). A bimodal dispersal profile (high along streams and low across country) would necessitate analysis at the landscape level. Although the species distribution model from just the main parameter precipitation does not correspond well with documented occurrences that delimit the range border, at least over the southern part of the species’ range (Figure 3c), it highlights the possible mosaic of habitat suitability across the subspecies hybrid zone.

Steep and coincident clines could also arise from a strong physical barrier to dispersal (Barton & Gale, 1993). Collectively, these clines coincide with an area of low habitat suitability (Figure 3c). As noted earlier, the connection between the southern and northern subspecies appears to be limited to a narrow and seemingly interrupted corridor directly east of the city of Coimbra. The combined results from genetic and distribution modelling are in accordance with the presence of a tension zone, for which is predicted that clines cluster in population density troughs or areas of restricted dispersal, even when effective selection against hybrids is low (Barton, 1980; Bierne et al., 2011; Buggs, 2007; Endler, 1977; Gompert et al., 2017).

Our study revealed a set of markers with their cline centre displaced northward into the C. l. longipes geographic range (13 nuclear and mtDNA), whereas not a single marker shows a southward displacement. Asymmetric introgression could arise from pre‐zygotic effects, including sex‐biased dispersal and non‐random mate choice (García‐París et al., 2003; Lamb & Avise, 1986; Petit & Excoffier, 2009) and post‐zygotic effects, including sex‐biased offspring survival (Haldane, 1922; Laurie, 1997) and Dobzhansky–Muller incompatibilities (Dobzhansky, 1937; Turelli & Orr, 2000). Field data on Chioglossa showed that observed travel distances were approximately equal for both sexes (Arntzen, 1984), alluding to the absence of sex‐biased dispersal, but the occurrence of asymmetric assortative mating or post‐zygotic effects in Chioglossa remains unknown. Moreover, asymmetric introgression can be explained by other processes, such as adaptive introgression or hybrid zone movement (Barton & Hewitt, 1985; Buggs, 2007; Currat et al., 2008; Toews & Brelsford, 2012; Wielstra, 2019). The observed introgression may be caused by C. l. lusitanica alleles that perform well in the C. l. longipes genetic background and that introgressed ahead of neutral ones. However, the proportion of displaced markers is high (14 out of 36) and the pattern of introgression is asymmetric, suggesting that some kind of trend is operating rather than differential selection on individual loci (Barton & Hewitt, 1985; Bierne et al., 2011; Currat et al., 2008; Wielstra, 2019; Wielstra, Burke, Butlin, Avcı, et al., 2017). Accordingly, the observed biased introgression, with the presence of southern C. l. lusitanica‐specific alleles over a ca. 50 km long stretch north of the centre of the hybrid zone, is best interpreted as the signature of a moving hybrid zone, in which C. l. longipes incorporates genes from C. l. lusitanica when expanding its territory. Additionally, we found that most of the nuclear displaced clines are wide, suggesting that introgressed alleles behaved neutrally once left behind, as expected when hybrid zones move (Barton & Gale, 1993; Macholán et al., 2011; Wielstra, Burke, Butlin, Avcı, et al., 2017). The clines come in three sets, with positions at the southern (a), middle (b) and northern section of the corridor (c) (Figure 3d). Group (a) represents the centre of the hybrid zone, whereas sets (b) and (c) are interpreted as genetic remnants of the receding southern lineage, but to associate these one to one with environmental features such as altitude and precipitation would require more spatial detail (Figure 3c).

The two Chioglossa subspecies have been diverging within Iberia since the early Pleistocene (Alexandrino et al., 2002; see Gomez & Lunt, 2007 for a ‘refugia within refugia’ model). Because refugia were local, minor range fluctuations will have triggered repeated instances of spatial contact and genetic exchange, allowing for recombination in the subspecies’ differentiated gene pools, therewith eroding regions of the genome involved in hybrid incompatibilities and/or conferring adaptation to local habitat conditions (Baird, 1995; Bierne et al., 2011). Accordingly, in the absence of strong selection limiting gene exchange, it is possible that innate population size disparities between C. lusitanica subspecies underlie the movement of the hybrid zone, as was suggested for other hybrid zones from empirical and simulated data (Barton & Hewitt, 1985; Engler et al., 2013; Menon et al., 2020). Considering that fragmentation of the suitable habitat across the hybrid zone could have played an important role in enhancing barriers to dispersal, we propose that random genetic drift in more or less isolated populations delayed the advance or even trapped some genetic clines, dissociated from the southward moving unit (Barton, 1983; Plouviez et al., 2013; Souissi et al., 2018). In line with theoretical predictions, it is more likely that barriers to gene flow or allele frequency shifts correspond to a pre‐existing tension zone that has become trapped by a natural barrier to dispersal (Barton, 1979; Bierne et al., 2011). Indeed, at ca. 40%, the proportion of displaced clines found in Chioglossa is too high to be explained by local adaptation (Bierne et al., 2011). The bimodal dispersal ability of Chioglossa in combination with a fragmented habitat could lead to the long‐term persistence of this genetic mosaic.

Although hybrid zone movement has not previously been considered to be common (Buggs, 2007), the spatial and genetic signatures of moving hybrid zones are by now been well documented in a wide range of vertebrate taxa including fishes (Souissi et al., 2018), salamanders (Arntzen & Wallis, 1991; Wielstra, Burke, Butlin, Arntzen, et al., 2017), toads (Arntzen, 1978, 2019; Arntzen et al., 2017), birds (Carling & Zuckerberg, 2011) and mammals (Lado et al., 2018; Macholán et al., 2011). Two important reservations apply to our interpretations. First, given the limited opportunity for sampling we cannot rule out that clines in the south were missed, just as we would not have encountered some displaced clines if the northern leg of the transect had been short. Second, we are forced to rely on a single transect whereas several studies on multiple transects across the same hybrid zone revealed different patterns of introgression, governed by a variety of demographic and environmental parameters (Dufresnes et al., 2021; Harrison & Larson, 2016; Janoušek et al., 2012; Szymura & Barton, 1986, 1991; Teeter et al., 2010). Further analysis of additional transects (if these were to exist) and extending the sampling of the southern populations could help to test and refine our inferences. Nonetheless, the observation on unidirectional introgression for ca. 40% of the genetic markers provides a strong and consistent signal for a moving hybrid zone between the two Chioglossa subspecies.

The southward movement of the Chioglossa hybrid zone has previously been suggested (Sequeira et al., 2005). However, with the increased number of markers and the incorporation of a species distribution model, the current study shows the importance of considering the influence of habitat availability on the spatial signal of introgression, effectuated through a moving hybrid zone. We inferred that the fragmented habitat across the Chioglossa hybrid zone acts as a semi‐permeable barrier, with allele frequency shifts in two pockets positioned within the wider area that changed occupation by the two taxa. This scenario is similar in pattern and process to that of pockets of resistance (or ‘enclaves’) formed by one hybridizing taxon surrounding and replacing the other as a function of the environment. Such enclaves may be taken as an indication of a moving hybrid zone (Arntzen et al., 2021; López‐Delgado et al., 2021).

The Chioglossa case appears as a particularly intricate example of hybrid zone movement, illustrating the important evolutionary role of the environment/habitat in shaping present day patterns of introgression. Our findings highlight the potential confounding effects of (sub)species range border, hybrid zone dynamics and environmental factors in inferring the strength of reproductive isolation, while illustrating how two subspecies maintain their genetic cohesion in face of gene exchange, while undergoing distribution range shifts over extended periods of time. For a note on the taxonomic status of C. l. lusitanica and C. l. longipes, see Appendix S2.

ACKNOWLEDGEMENTS

We thank Miguel Vences for help with sampling and Arie van der Meyden for the transport of material. Tissue samples were collected under license 420/2020/CAPT from the Portuguese Institute for Nature Conservation and Forests (ICNF).

CONFLICT OF INTERESTS

The authors declare no competing interests.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/jeb.13982.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are openly available in Dryad at https://doi.org/10.5061/dryad.7h44j0zvr.

REFERENCES

Author notes