A signature of dynamic biogeography: enclaves indicate past species replacement

Understanding how species have replaced each other in the past is important to predicting future species turnover. While past species replacement is difficult to detect after the fact, the process may be inferred from present-day distribution patterns. Species with abutting ranges sometimes show a characteristic distribution pattern, where a section of one species range is enveloped by that of the other. Such an enclave could indicate past species replacement: when a species is partly supplanted by a competitor, but a population endures locally while the invading species moves around and past it, an enclave forms. If the two species hybridize and backcross, the receding species is predicted to leave genetic traces within the expanding one under a scenario of species replacement. By screening dozens of genes in hybridizing crested newts, we uncover genetic remnants of the ancestral species, now inhabiting an enclave, in the range of the surrounding invading species. This independent genetic evidence supports the past distribution dynamics we predicted from the enclave. We suggest that enclaves provide a valuable tool in understanding historical species replacement, which is important because a major conservation concern arising from anthropogenic climate change is increased species replacement in the future.


Background
Species that meet in nature but geographically exclude one another have parapatric distributions [1]. While the position of parapatric range boundaries may be dynamic, their actual movement is a protracted process and has necessarily been recorded over shallow time frames only [1,2]. We here argue that the current spatial arrangement of parapatric species can be informative about distributional shifts in the more distant past. In particular, we highlight the insight provided by enclaves. An enclave is formed when part of the distribution of one member of a pair of parapatric species is isolated inside the range of the other ( figure 1). An enclave could originate via colonization-analogous to a species establishing a peripatric population on an oceanic island [3,4]. Alternatively, the enclave might be a relict distribution patch of a previously broader distribution, left after incomplete species replacement-akin to rising sea levels disconnecting a continental island population from the mainland [5,6]. The likelihood of either scenario depends on the interplay between the mean and shape of a species's dispersal kernel [7] and the enclave's distance from the main distribution range [8]. We assert that enclaves of species with low vagility, relative to the distance between the enclave and the main range, provide strong a priori support for species replacement. We provide proof of concept for an enclave system involving newts.
We take advantage of the fact that contact zones for parapatric species tend to correspond to hybrid zones [1]. When, upon secondary contact, a hybrid zone is first established, and one of the species possesses a competitive edge over the other one, this would result in a non-equilibrium phase of species replacement, with the hybrid zone moving as a consequence [1,4,5]. Furthermore, for a hybrid zone initially stabilized at an environmental gradient, this equilibrium phase could be disrupted if climate change shifts the balance in favour of one of the species [9]. If the species involved show introgressive hybridization, a moving hybrid zone is predicted to leave a trail of unlinked, selectively neutral alleles in its wake, derived from the displaced species and inside the territory newly claimed by the expanding one (initially based on theoretical principles [10], later supported by simulation [11] and recently backed up by empirical findings [12]). Under this rationale, the hypothesis of enclave formation can be tested independently in hybridizing species, based on the geography of genome-wide interspecific gene flow: an enclave formed by incomplete species replacement is expected to be accompanied by a genomic footprint of hybrid zone movement, while an enclave formed by colonisation is not consistent with this pattern (figure 1).
We examine an enclave observed in crested newts (Amphibia: Triturus), located in central Serbia, on the Balkan Peninsula [13]. Part of the range of T. ivanbureschi (blue in figure 2a) is detached from the main distribution, because the range of another species, T. macedonicus (green), intervenes. The ranges of two additional species, T. cristatus (red) and T. dobrogicus (orange), border that of T. ivanbureschi in the north. Both T. ivanbureschi and T. macedonicus colonized the central Balkan Peninsula postglacially, as they expanded from discrete glacial refugia, distant from the enclave [15]. Triturus cristatus was already in place as it had a glacial refugium in the southern protrusion of its current range, where it borders T. ivanbureschi today [16]. Although T. dobrogicus expanded its range postglacially [17], it is, uniquely among crested newts, a specialized lowland species [18], and its distribution in the south is bounded by an elevational ecotone (figure 2b). Therefore, while T. cristatus or T. dobrogicus represent a biological barrier north of the enclave, it is unlikely that either species invaded the range of T. ivanbureschi over a substantial area.
As the distance between the enclave and the main range of T. ivanbureschi (approx. 83 km) is over 20 times larger than lifetime dispersal distance (approx. 3.7 km [12]), incomplete species replacement (rather than long-distance colonization) is the likely explanation for the enclave's origin. Because crested newt species ranges meet at hybrid zones that facilitate introgression [19], we are in a position to test for a genomic footprint of hybrid zone movement around the enclave. We predict that such a genomic footprint was left by T. ivanbureschi, where we assume it was superseded by T. macedonicus. We do not predict a pronounced asymmetry in introgression between T. ivanbureschi and either T. cristatus or T. dobrogicus, as we expect T. ivanbureschi's hybrid zones with these species to have been relatively stable.

Material and methods (a) Sampling
We included 1 -3 individuals (mean 2.645) from 251 Triturus localities; 664 individuals in total (figure 2a; electronic supplementary material, table S1). For part of these individuals the full set of sequence data (n ¼ 308) or (mt) DNAs (n ¼ 283) were available from previous studies [12,17,[20][21][22][23][24], while additional individuals (n ¼ 73) from Serbia and Bulgaria are studied here for the first time. We selected 15 reference individuals per species-three individuals from five localities positioned away from contact zones and presumed to be unaffected by interspecific gene flow-to determine whether individual markers are diagnostic (electronic supplementary material, table S1). We created a Voronoi diagram in ARCGIS 10 (www.esri.com), where each locality is represented by a Thiessen polygon that contains the area that is closer to that particular locality than to any other one (figure 2c).

(b) mtDNA sequencing and analysis
We Sanger sequenced a 658 bp mtDNA fragment for 625 individuals (following [17]) and were able to sequence a 110 bp internal fragment for 23 individuals following [12]. The clear-cut geographical distribution of mtDNA allowed us to infer mtDNA type for the 16 remaining individuals (electronic supplementary material, table S1). The 658 bp mtDNA fragments were added to the mtDNA haplotype database of [17] and collapsed into haplotypes with MACCLADE 4.08 [25]. To assign new haplotypes to species, we constructed a neighbour-joining phylogeny with 1000 bootstrap replicates in MEGA 5 [26], with the Pyrenean newt Calotriton asper and the marbled newt T. marmoratus (from [27]) as outgroups (electronic supplementary material, figure S1; table S2). Internal 110 bp mtDNA fragments were aligned with the set of 658 bp haplotypes and this set was then trimmed accordingly. By removing redundancy in MACCLADE, the internal fragments could be allocated unambiguously to species-diagnostic mtDNA type.

(c) Sequencing nuclear DNA
For all 664 individuals, we sequenced 52 nuclear markers following [22]. Details are in electronic supplementary material, text S1.

(d) Bayesian clustering analysis
We used STRUCTURE 2.3.3 [28] to estimate the fraction of ancestry for each individual derived from the four parental species based on nuclear DNA data. We used the admixture model in combination with the correlated allele frequency model with 1 000 000 iterations, after 250 000 iterations of burn-in, and ran 10 replicates. An initial STRUCTURE analysis on the set of reference individuals, in which we allowed the number of genepools k to vary from 1-20 (with the upper limit defined by the total number of localities included), confirmed k ¼ 4 as the most likely number of genepools under  Figure 1. Two hypotheses on enclave formation in a pair of parapatric species. In (a) a blue species founds a population within the range of a green one via long-distance colonization, by leapfrogging across a stable contact zone. In (b) a relict population of a receding blue species persists locally within the range of a superseding green species, behind a moving contact zone, which was initially positioned at the grey dotted line. If the two species hybridize, a genomic footprint of hybrid zone movement would be expected under the moving contact zone scenario (in the hatched area in (b)), but not under the stable contact zone scenario. (Online version in colour.) rspb.royalsocietypublishing.org Proc. R. Soc. B 284: 20172014 Evanno's Dk criterion [29], as implemented in CLUMPAK [30]. Next we conducted a STRUCTURE run on the full set of individuals, in which we fixed k to 4. The lowest Q-value with which a reference individual was allocated to its respective species was Q ¼ 0.9794. Localities allocated to two (or in some cases three) species with Q . 0.0206 (i.e. 1-0.9794) were considered genetically admixed. STRUCTURE scores (with replicates summarized in CLUMPAK) are presented in electronic supplementary material, table S1 and locality averages are plotted in figure 2b.

(e) Hybrid index, heterozygosity and ancestry, and test for asymmetric introgression
We determined the evolutionary origin of alleles based on our 15 reference individuals per species (see section on sampling). We found 23 nuclear DNA markers to be diagnostic, i.e. exhibiting fixed allelic differences between T. ivanbureschi (the species with the enclave) and the other three species. We determined the proportion of diagnostic T. ivanbureschi alleles present at each locality, i.e. the mean hybrid index (electronic supplementary material, table S1), and plotted this on a map (figure 2c). For pairwise species comparisons of T. ivanbureschi versus the other three Triturus species, we determined individual heterozygosity (the fraction of markers heterozygous for alleles from each parental species) and ancestry (the fraction of alleles derived from each parental species) using the R [31] package 'HIest' [32]. Locality averages (electronic supplementary material,  determined the 0.5 contour for the hybrid index using the 'akima' package [33] in R. We calculated the minimum straight-line distance of each locality from that contour in ARCGIS, and gave distances for T. ivanbureschi a positive sign and distances for other species a negative sign (electronic supplementary material, table S3). For the resulting one-dimensional transect [34] we then fitted a geographical cline to the hybrid index in the R package 'HZAR' [35]. Frequency was fixed at the ends of the cline to 1 (pure T. ivanbureschi) and 0 ( pure for the other species). The hybrid index is the proportion of diagnostic alleles, calculated from 23 diploid markers (minus a few missing genotypes) in one to three individuals. Therefore, it was fitted with binomial error distribution and sample size equal to the total number of genotypes (loci Â individuals) per locality. Although markers are unlikely to be closely linked, they may not behave independently. Using the number of genotypes rather than the number of alleles makes our test conservative, despite this constraint. We tested five models, differing in complexity: no tail, left tail only, right tail only, mirror tails or both tails estimated separately. We used the lowest Akaike information criterion score corrected for small sample size to select the best one; in each case the score of the chosen model was at least two likelihood points lower than that of the next best fitting model ( figure 4; electronic supplementary material, table S6).   The asymmetry in nuclear DNA introgression from T. ivanbureschi into T. macedonicus is statistically significant, based on either the presence (Fisher's exact test, p ¼ 0.000) or the intensity (one-tailed Mann-Whitney U test, U ¼ 1423.5; Z ¼ 23.11; p ¼ 0.000) of introgression. There is no such asymmetric nuclear DNA introgression into T. cristatus ( p ¼ 0.554; To remove the influence of a potential bias in sampling, we collapse the two-dimensional sampling into a one-dimensional transect, based on the shortest distance of localities to a 0.5 hybrid index contour. Onedimensional geographical clines, fitted to a hybrid index derived from the 23 diagnostic nuclear DNA markers for T. ivanbureschi versus each of the other three species, differ in the models that best fit the data. Only for T. macedonicus was an asymmetric cline model preferred, with a substantial introgression of T. ivanbureschi alleles into T. macedonicus, but not in the opposite direction (figure 4).

Discussion
Guided by the enclave observed in our crested newt system, we hypothesized a past parapatric range shift, in which T. macedonicus intersected the range of T. ivanbureschi and excised the enclave. As the two species must have hybridized in the process, we could test this scenario by exploiting the pattern of interspecific gene flow. In line with predictions for moving hybrid zones based on theory [10] and data simulation [11], we exposed a genomic footprint of hybrid zone movement, left by T. ivanbureschi, where it has been replaced by T. macedonicus. While cytonuclear incompatibilities or dispersal effects might cause asymmetry in mtDNA introgression and selection might cause asymmetry for any given locus (while not necessarily in the same direction), a genome-wide shared asymmetry in introgression over an extensive area strongly favours a scenario of species displacement with hybridization [12]. Hence, the crested newt case illustrates the predictive power of enclaves for inferring past species replacement.
Enclaves are not often reported, but, considering that three are known in Triturus newts alone (see [36,37] for two other species pairs), this does not necessarily mean that they are rare. Enclaves may well be transient, with species replacement not prevented but merely postponed locally. Indeed, all individuals from the T. ivanbureschi enclave show some degree of genetic admixture with other crested newt species. The ephemeral nature of enclaves is illustrated by a previously documented enclave (established from museum specimens) having been taken over by an invader over a period of roughly 50 years [38]. Furthermore, 'islands of alleles' derived from a displaced species that was genetically swamped out by a competitor, could be regarded as genetic fossils of a past enclave [39]. In general, being small and isolated from the main distribution makes enclave populations more susceptible to the extinction risks associated with habitat fragmentation [40].
Once an enclave is identified, it can provide crucial insights into historical biogeography. While hybrid zones have been tracked at the scale of decades [4], they are thought to stabilize quickly at environmental density troughs [10,41]. Enclaves, however, indicate hybrid zone movement over extended periods of time-an opportunity that has as yet received little-to-no attention in the hybrid zone literature [4,11]. Understanding past shifts in the mutual range boundaries of parapatric species is particularly relevant in light of the expected future increase in such shifts in response to manmade habitat alteration and climate change [9,42]. Data accessibility. The data underlying this paper are archived in the Dryad Digital Repository (https://doi.org/10.5061/dryad.9k1h0) [43].