Introgressive hybridization between the endangered native Bombina pachypus and the introduced B. variegata in a protected area in central Italy

. Amphibians are experiencing an ascertained global decline, which causes include the introduction of alien species and the (anthropogenic) hybridization between native and exotic taxa. Detecting introductions and assessing their impact on populations of native species is crucial for amphibian conservation. We used mitochondrial and nuclear markers to reveal introgressive hybridization between the native Bombina pachypus and the exotic B. variegata (probably introduced from Albania) in a population from a protected area of central Italy. Almost all genotyped individuals were genetically admixed, showing a larger proportion of the allochthonous genome. To our knowledge, this study provides the ﬁrst evidence of successful hybridization between the two species (we found both putative F1 and backcrosses), hence representing a new threat to the conservation of the endangered, Italian-endemic B. pachypus .

The introduction of alien species is a major concern for biodiversity conservation. Negative consequences of such phenomenon include the extensive hybridization between exotic taxa and their native relatives, which may result in a serious threat to the persistence of autochthonous taxa (Rhymer and Simberloff, 1996;Allendorf et al., 2001;Todesco et al., 2016). Both ecological and genetic factors may promote the spread of exotic genes, ultimately yielding, in some cases, to the rapid replacement of the native genome (Fitzpatrick et al., 2010;Todesco et al., 2016). Alternatively, gene-pool contamination through introgression may cause deleterious effects, among them, an overall reduced fitness of natural populations (Rhymer and Simberloff, 1996). Thus, identifying early-stage introductions and hybridization events between exotic and native species is essential in a conservation perspective (Frankham, Briscoe and Ballou, 2010), especially when endangered species are involved.
Alien species also play an important role in the worldwide amphibian decline. Both intentional and accidental introductions/translocations of amphibians are frequent (Smith and Sutherland, 2014), often as a consequence of human activities (e.g. pet trade, farming, human consumption). However, the impact of hybridization between autochthonous and allochthonous taxa on native populations may differ case by case, with local or regional effects observed (Beebee and Griffiths, 2005). In Europe several human-mediated introductions have been reported, sometimes involving the hybridization between native and exotic amphibians (e.g. Dufresnes et al., 2015Dufresnes et al., , 2016Meilink et al., 2015;Palomar, Vörös and Bosch, 2017).
Bombina pachypus (Bonaparte, 1838) is an endangered anuran, which distribution spans the Italian peninsula, south the Po river valley. Populations of B. pachypus are experiencing recent and severe declines probably caused by multiple factors, among which the spread of chytridiomycosis (Andreone et al., 2009;Canestrelli, Zampiglia and Nascetti, 2013). However, hybridization with exotic relatives has not caused concern to date. Indeed, the European congeneric species -i.e. B. variegata and B. bombina, which naturally hybridize in their contact zone Barton, 1986, 1991) -are geographically isolated from B. pachypus and introductions of related species within B. pachypus range have never been reported.
Here we carried out a genetic assessment of a Bombina population in a protected area in central Italy where we found toads showing intermediate morphological traits between the native B. pachypus and the exotic B. variegata -the two species are morphologically similar, although differences in the ventral colour pattern have been described (Vaccaneo, 1931). Therefore, we aimed to assess the taxonomic status of the population and to investigate the genetic background of suspected hybrids. We used mitochondrial (cytochrome b) and nuclear (recombination-activating gene 1, microsatellites) markers to achieve our purposes.
During the 2014 breeding season, 21 Bombina samples (buccal swab and tail tips from 17 toads and four tadpoles, respectively) were collected from the supposed hybrid population in the Regional Park of Lucretili Mountains (municipal district of Licenza), hereinafter referred to as LIC ( fig. 1, supplementary information S1). Genomic DNA was extracted using the GenElute™ Mammalian Genomic DNA Miniprep Kit (Sigma-Aldrich, St. Louis, MO, USA), following the mammalian tissue protocol.
We used fragments of the mitochondrial cytochrome b (cytb; 1096 bp) and the nuclear recombination-activating gene 1 (RAG-1; 1061 bp) to define the genetic background of LIC individuals and trace their origin. Indeed, the geographic distribution of haplotypes/alleles at those markers is available for European Bombina species (data from Hofman et al., 2007;Fijarczyk et al., 2011). Additionally, two pure B. pachypus from a locality six km away from LIC (i.e. pop 21; supplementary information S1) and 25 pure Bombina variegata (i.e. 17 individuals from 11 localities in northern Italy and eight individuals from a site in Albania; supplementary information S1) were genotyped for both markers and included in the reference dataset of sequences. Amplification, sequencing and genotyping procedures are described in supplementary information S2.
LIC individuals and 25 B. variegata samples were genotyped for nine microsatellites loci -F2, B14 and F22 (Stuckas and Tiedemann, 2006); 8A, 9H, 12F, 1A, 5F and 10F (Hauswaldt, Schroder and Tiedemann, 2007 (Rousset, 2008). The admixture in LIC individualsnamely the genome proportion belonging to B. pachypus or B. variegata -was inferred in STRUCTURE v2.3.3 (Pritchard, Stephens and Donnelly, 2000). We run the USEPOPINFO model (analysis settings: K = 2; allele frequencies correlated and updated based on reference samples; 1 000 000 of Monte Carlo replicates after a burn-in period of 200,000) which uses pre-defined reference groups (i.e. genotypes of putative parental species) to assist admixture estimation for individuals of unknown origin (i.e. 21 LIC individuals). Our two reference groups, i.e. 25 B. variegata and 149 B. pachypus, were unambiguously validated by a preliminary STRUCTURE analysis (supplementary information S3). Additionally, the Snapclust function in the adegenet Rpackage (Jombart, 2008) was used to calculate membership probabilities of LIC individuals to each of five classes (pure parentals, F1-hybrid and first-generation backcrosses) according to a combined approach of geometric/fast likelihood optimization (Beugin et al., 2018). An individual was univocally assigned to a class if the most likely class inferred by Snapclust was consistent with its cytb/RAG profile. Otherwise, it was considered as unclassified. The analysis involved microsatellite genotypes of overall 195 Bombina: pure B. pachypus and B. variegata were given as two baseline groups to improve the assignment accuracy of LIC individuals. Finally, we assess and visualize genetic relationships among all samples performing a microsatellite-based principal component analysis (PCA) in adegenet.
We found two cytb haplotypes (A3, BW21) and three RAG-1 alleles (R3, R6, R14) in LIC (supplementary information S1). All of which were detected in previous studies and effectively distinguished between B. variegata and B. pachypus. Specifically, BW21, R3 and R6 were found among B. variegata populations from the Balkan peninsula, co-occurring in Albania only (i.e. PUK; fig. 1). By contrast, A3 and R14 were widespread across B. pachypus range ( fig.  1), although the latter was shared with two B. variegata populations from Bulgaria. However, R14 occurred even in localities close to LIC (i.e. in pop 21 and 135; fig. 1b), hence conservatively, it was considered a local B. pachypus allele. Fourteen LIC individuals showed sequences from both species ( fig. 2a, supplementary information S1), indicating their hybrid nature. The remaining seven individuals exhibited pure B. variegata mitochondrial-nuclear genotypes.
All microsatellite loci were polymorphic in LIC (2 to 4 alleles per locus). We found only a weak tendency of heterozygosity excess in seven out of nine loci, but none departed significantly from the Hardy-Weinberg equilibrium after applying the Bonferroni correction for multiple testing (P > 0.05). The microsatellite-based STRUCTURE analysis revealed that all LIC individuals were genetically admixed ( fig. 2a, supplementary information  S4), with a prevalent proportion of the B. variegata genome (range = 34.8%-80.3%; mean ± SD = 60.2 ± 11.8%). Applying the crossvalidation, Snapclust allowed to confidently assign 18 out of 21 LIC individuals to a putative hybrid/parental class: a single pure B. variegata, three F1 hybrids and 14 first-generation backcrosses (three F1 × B. pachypus, 11 F1 × B. variegata) ( fig. 2b; supplementary information S4). Finally, B. pachypus, B. variegata and LIC samples formed three distinct groups in PCA, with first and second PCA-axis distinguishing variation among and within species, respectively ( fig. 2c). Most of LIC individuals showed affinity to B. variegata samples, consistently with STRUCTURE and Snapclust results, particularly those from Albania. Molecular data, according to morphological observations, unequivocally revealed that almost all individuals from LIC are hybrids of various classes between local B. pachypus and exotic B. variegata. Given the known distributions of allochthonous sequences found in LIC (BW21 for cytb; R3 and R6 for RAG-1), we hypothesized the Balkans, specifically Albania, as the most plausible source region of introduced B. variegata. However, the circumstance of such a human-mediated introduction remains unknown. Anthropogenic introductions of amphibians and their hybridization with native taxa have been reported frequently in Europe, even within protected areas. For instance, Palomar and collaborators (2017) demonstrated the introduced origin of an Ichthyosaura alpestris population in the Guadarrama National Park (Spain). Similarly, evidence for introgression in Hyla arborea and Triturus cristatus by human-introduced Hyla intermedia and Triturus carnifex, respectively, was found in Switzerland (Dufresnes et al., 2015(Dufresnes et al., , 2016 and Netherlands (Meilink et al., 2015).
Our results imply the incomplete reproductive isolation between B. pachypus and B. variegata and suggest fertility of their hybrids (both first-generation backcrosses were inferred). Also, the occurrence of maternallyinherited mitochondrial haplotypes of both species in hybrids implies bidirectional hybridization. This is not particularly surprising since the two taxa are closely related: although they were previously considered a single species, B. pachypus is now recognized as a genetically isolated taxon, which evolved in allopatry for 2-3 million years ago (Pabijan et al., 2013). Besides, even the more distantly related B. variegata and B. bombina (estimated divergence = 6-9 million years ago; Pabijan et al., 2013) extensively hybridize in the wild (Szymura and Barton, 1986Barton, , 1991. However, to our knowledge, successful hybridization between B. pachypus and B. variegata has not yet been reported. Oddly, we did not detect any pure native B. pachypus in LIC according to mitochondrial/nuclear sequence patterns, while the Snapclust analysis inferred a single B. variegata. At least two hybrid generations (F1 and firstbackcrosses) co-occurred in LIC, and most of the hybrids showed a prevalent component of B. variegata genome ( fig. 2). A possible explanation of this pattern is that introduced B. variegata individuals may have outnumbered the native ones, therefore B. pachypus genes could have been diluted through successive generations of backcrossing -note that we were not able to confidently explore multi-generation backcrosses in Snapclust because of the limited number of loci. Besides, we cannot rule out that introduced individuals were only hybrids of various generations. Alternatively, B. variegata and/or hybrids between the two species may perform better compared to the local B. pachypus, hence favouring the introgression of exotic genes. For example, a reduced resistance against the emerging chytridiomycosis has been hypothesized for B. pachypus compared to B. variegata based on the different conformation of major histocompatibility complex-encoded proteins (Talarico et al., 2019). On the other hand, Szymura and Barton (1986) demonstrated lower viability of B. bombina × B. variegata embryos as opposed to pure embryos. Although fascinating, the hypothesis of selection favouring hybrids and/or B. variegata remains a conjecture without experimental support. Moreover, selection efficacy is expected to be reduced in small populations, such as LIC, and even alleles with different coefficients of selection may behave as neutral. However, the abovementioned explanations are not mutually exclusive.
Hybridization likely occurred a number of years before sampling (>7, according to earliest morphology-based observations; Tiberi, pers. comm.), as suggested by the occurrence of multiple hybrid generations, the almost-lack of pure individuals, and the non-significant excess of heterozygosity in microsatellite loci. The latter could be an artefact due to the small sample size, anyway a population is predicted to return to Hardy-Weinberg equilibrium after a single generation of random mating (in nonoverlapping generations), or a relatively small number of generations (in overlapping generations).
The LIC population is a new case of hybridization in the wild within Bombina, raising interesting questions about the competitive ability of hybrids/invaders as compared to that of native relatives. Nevertheless, in a conservation perspective, this hybrid population poses a threat to the genetic integrity of native populations because of its anthropogenic nature and the massive introgression from B. variegata, a process known as hybrid swarm (Allendorf et al., 2001;Frankham, Briscoe and Ballou, 2010; but see Pabijan et al., 2020 for opportunities related to hybridization for conservation purposes). Conservatively, our results stress the need for a conservation strategy aimed at preventing dispersal of hybrids and restoring the local gene pool. Specifically, we suggest the removal of hybrids, possibly followed by translocations of pure toads, tadpoles or clutches from suitable, closely-related populations of local B. pachypus (e.g. pop 21). Contextually, we strongly recommend genetic monitoring of native populations in the surrounding areas, in order to verify their genetic integrity. Genetic monitoring of wild populations is crucial for the conservation of endangered species, even to detect unpredictable cases of hybridization with cryptic exotic species, as this study highlighted. Ultimately, our findings contribute to understanding genetic consequences of anthropogenic introductions on populations of autochthonous amphibians.