Distribution and evolution of the western European water frogs (genus Pelophylax) from Catalonia, northeastern Spain

Sampling

A total of 423 water frogs from 54 localities were collected between January to June 2018, at the peak of the breeding season of Pelophylax (Fig. 2 and Table S1). In addition, samples from the natural ranges of P. ridibundus kurtmuelleri (n = 1, Greece), P. lessonae (n = 1, France), and P. ridibundus (n = 1, obtained from an introduced population in Riyad, Saudi Arabia) were included. Sampling was designed to cover all main river basins within Catalonia, Spain, placing special focus in the closest basins to the French border, given the possible natural entrance of the klepton from southern France into the Iberian Peninsula. Populations connected through rivers were selected so that there was a minimum distance of 10 km between them. However, this minimum distance was not considered when localities were separated by a geographical barrier.

To ensure that the sampling effort was similar at all localities Iberian water frogs were captured over a 20-minute span. Individuals collected were identified to species level according to general morphological characteristics such as size, color, head shape, snout shape, extent of the interdigital membrane, size and shape of the metatarsal tubercle, and, especially, by the shape and distance between the proximal edges of the vomerine ridges (Crochet et al., 1995; Duguet & Melki, 2003; Pérez-Sorribes et al., 2018; Plötner, Becker & Plötner, 1994; Plötner, 2010). Tadpoles and doubtful specimens were classified as Pelophylax sp. All specimen identifications were conducted by the same people (CM-D and DF-G). For every collected specimen live stage (adult, juvenile, or tadpole) and sex (male, female, or undetermined; Table S1) were recorded, and a tissue sample was obtained by toe-clipping in case of adults or juvenile specimens, or by a five mm tail tip in case of tadpoles. Tissue samples were preserved in absolute ethanol and stored at −20 °C. To minimize blood loss, a tissue adhesive produced from enbucrilate (Histoacryl, B. Braun Surgical) was immediately applied to the wound. All collected individuals were released at the same collecting spot after manipulation. Manipulation of all amphibian specimens followed a strict biosecurity protocol (Protocol Sanitari d’Amfibis de la Generalitat de Catalunya) established to avoid the spread and infection of emergent diseases (Berger et al., 1998; Martel et al., 2013).

Collection and manipulation of the specimens included in this work was approved by the Servei de Fauna i Flora del Departament de Territori i Sostenibilitat de la Generalitat de Catalunya (permit numbers: SF/0040, SF/0041, SF/0042, SF/0043). A special permission was acquired for collecting inside the following protected natural areas: Parc Natural del Delta de l’Ebre; Parc Natural dels Ports; Parc Natural del Montgrí, les Illes Medes i el Baix Ter; Parc Natural dels Aiguamolls de l’Empordà; Parc Natural d’Interès Nacional de l’Albera; Parc Natural de la Zona Volcànica de la Garrotxa; Espai d’Interès Natural de l’Alta Garrotxa. Regarding the sampling, study design, and sample manipulation of wildlife specimens, this study adheres to the ARROW guidelines (Field et al., 2019).

Molecular species identification

Genomic DNA was extracted from the alcohol-preserved tissues following a standard high-salt protocol (Sambrook, Fritsch & Maniatis, 1989). A two-step molecular protocol was then performed to differentiate between pure P. ridibundus, pure P. perezi, and the klepton P. kl. grafi, as well as differentiate them from any other recognized extant species and kleptons of Western European water frogs.

The first step was based on the different cutting pattern of the PCR amplified region of the Recombination activating gene 1 (RAG1) with the restriction enzyme EcoO109I. The different number and disposition of cleavage sites between Pelophylax species and kleptons allow for the discrimination of Pelophylax into three groups: (1) One band in P. perezi, P. lessonae and P. bergeri; (2) Two bands in P. ridibundus and P. kurtmuelleri; (3) Three bands in P. kl. grafi, P. kl. esculentus, P. kl. hispanicus and PK (Fig. S1). The RAG1 nuclear region was amplified using the same primers as in Newman et al. (2012) MartFL1 (5′-AGC TGC AGY CAG TAC CAC AAA ATG-3′) and AMPR1rana (5′-AAT TCA GCT GCA TTT CCA ATG TC-3′) and PCR conditions were as follows: initial denaturation step at 95 °C for 5 min; 41 cycles of denaturation at 94 °C for 1 min, annealing at 59 °C for 1.5 min and elongation 72 °C for 1 min; final extension step at 72 °C for 10 min. Digestion of the amplified RAG1 bands was carried out in a final volume of 25 µL, containing 1x EcoO109I Buffer, 10 units of EcoO109I enzyme, eight µL of the amplified PCR product, and the remaining volume up to 25 µl of H2O. The digestion mix was incubated at 37° for 1.5 h and inactivated at 65° for 20 min. The final product was run on a 1.5% agarose gel for visualization and band identification.

Although, the first step already allowed for a reliable identification of all Pelophylax species present in Spain and large areas of Southern France (P. ridibundus, P. perezi, and P. kl. grafi), a second step was performed to ensure a final discrimination at the specific level for each specimen. This step relied on the different sizes of the PCR amplified fragment of the serum albumin intron-1 (SAI-1), which contains a non-LTR retrotransposon of the chicken repeat family (RanaCR1; Cuevas et al., 2022; Hauswaldt et al., 2012). While some studies (Cuevas et al., 2022; Dufresnes et al., 2024b) have discussed potential limitations of the SAI-1 marker, particularly for distinguishing P. lessonae from P. ridibundus, this should not affect our study, as the former is not expected to occur in the study area. The SAI-1 region was amplified using the same primers as in Hauswaldt et al. (2012) with the following PCR conditions: Initial denaturation step of 94 °C for 1.5 min; 35 cycles of denaturation at 94 °C for 30 s, annealing at 59 °C for 40 s and elongation at 72 °C for 100 s; final extension at 72 °C for 10 min. A total of three µL of PCR product was run on a 1.5% agarose gel. As a result of the small size difference between P. kurtmuelleri, P. ridibundus, and P. perezi, apart from the 100 bp ladder we also loaded as a size control a mixture of amplified SAI-1 of P. ridibundus, P. kurtmuelleri, and P. lessonae.

Mitochondrial DNA (mtDNA) amplification and sequencing

A region of the cytochrome b (cyt-b) mtDNA was amplified for a representative sample (194 specimens) of all the specimens identified as kleptons including at least one sample from each locality where they were identified, to investigate the origin of their maternal DNA and the levels of mtDNA genetic diversity (Table S1). The cyt-b region of 647 bp was amplified using primers CYTB_PeloF (5′-ATG ATG AAA CTT TGG CTC C-3′) and CYTB_PeloR (5′-TCT GGC TTA ATG TGG GGG-3′) with the following PCR conditions: Initial denaturation step of 94 °C for 5 min; 35 cycles of denaturation at 94 °C for 30 s, annealing step at 56 °C for 45 s and elongation at 72 °C for 80 s; final extension step at 72 °C for 5 min. Sequencing of the PCR products was carried out by Macrogen Europe.

Phylogenetic and network analyses

For all 194 sequences, chromatographs were checked manually, assembled, and edited using Geneious v.9.0.5 (Biomatters Ltd.). Sequences were translated into amino acids and no stop codons were observed. A haplotype network was inferred for the 194 cyt-b sequences obtained in the present work using the TCS statistical parsimony network approach (Clement, Posada & Crandall, 2000) implemented in the program PopART (Leigh & Bryant, 2015). A subset of 22 specimens was then selected representing each retrieved haplotype for phylogenetic analyses (see Results). The 22 representative sequences were aligned together with a set of sequences downloaded from GenBank, including most of the described species of Pelophylax and a representative of Lithobates catesbeianus that was used to root the phylogenetic tree (Table S2; sequences from Dubey & Dufresnes, 2017; Lymberakis et al., 2007). The assembled dataset was aligned using MAFFT v.6 (Katoh & Standley, 2013) with default parameters. Phylogenetic analyses were performed using maximum likelihood (ML) and Bayesian inference (BI) methods. A GTR+G substitution model was selected under the Akaike information criterion (Akaike, 1973) using jModelTest v.2.10 (Posada, 2008), and ML phylogenetic analyses were conducted with RAxML v.7.4.2 (Stamatakis, 2006) with 100 random starting trees and 1,000 bootstrap replicates. BI analyses were performed with MrBayes v.3.2.6 (Ronquist et al., 2012). Two simultaneous runs were performed with four chains per run for 2 × 10⁶ generations sampling every 200 generations. The standard deviation of the split frequencies between the two runs and the potential scale reduction factor (PSRF) diagnostic were examined and convergence and stationarity were ensured (ESS > 200) using Tracer v.1.6 (Rambaut & Drummond, 2013). We conservatively discarded the first 25% of trees as burn-in.

Field sampling and preliminary morphological identification

Our sampling effectively covered the entire known distribution range of the genus Pelophylax in Catalonia (Fig. 2), including all major river basins. The final dataset comprised 423 water frogs spanning through 54 localities (Table 1). From all 423 Pelophylax, 125 were males, 65 females, 164 juveniles, and 69 tadpoles (Table 2). According to morphological identification, 192 specimens were classified as P. perezi, 79 as P. kl. grafi, and 152 could not be identified (mainly juveniles and tadpoles; Table 2 and S1). Genetic analyses revealed frequent misidentification of P. kl. grafi based on morphology in both sexes (28 adult P. kl. grafi males and 23 females were morphologically assigned to P. perezi). Similarly, 42 juveniles of P. kl. grafi were misidentified as P. perezi. Conversely, morphological assignment of P. kl. grafi was more reliable, with only four P. perezi adults (three males, one female) misidentified as P. kl. grafi (Table 2).

Molecular identification and geographic distribution

The two-step molecular identification method allowed for the unambiguous identification of all 423 sampled water frogs, including all tadpoles (Table 2 and Fig. 3, Figs. S2, and S3). Interestingly, despite the high number of specimens analyzed from a wide range of localities across Catalonia, only two species of water frogs were found, P. perezi and P. kl. grafi. As shown in Table 1, 23 out of the 54 populations sampled presented both P. perezi and P. kl. grafi in sympatry, with 17 populations including only representatives of P. perezi and 14 populations only representatives of P. kl. grafi. While both species are widely distributed across Catalonia, a relatively large region between Barcelona and the Ebro Delta seems devoid of P. kl. grafi (i.e., populations P4, P12, P14-P17, P21-P23 add up to 53 P. perezi and no P. kl. grafi; Fig. 3). Localities with exclusively P. kl. grafi tended to have small sample sizes (1–5 specimens), with only three localities (P25, P39, and P46) having larger sample sizes of 8, 16, and 16 specimens, respectively. No clear geographic pattern emerges in the distribution of the 14 populations containing only P. kl. grafi (Fig. 3), and the presence of P. perezi at low frequencies cannot be ruled out due to the limited sample sizes.

Mitochondrial DNA (mtDNA) analyses

A total of 194 specimens of P. kl. grafi were selected for molecular analyses, including at least one representative from each one of the 37 localities with P. kl. grafi specimens representing 78.5% of the total number of P. kl. grafi specimens identified in this work (Table 1). The selected specimens were sequenced for the cyt-b mitochondrial gene in order to know: (1) The mitochondrial composition of the kleptons; (2) the level of mtDNA genetic variability, and (3) the geographic distribution of the different haplotypes.

The results of the phylogenetic analyses on the cyt-b indicate that all 194 P. kl. grafi analyzed had P. perezi mtDNA (Fig. 4). A total of 22 unique haplotypes (H) were identified, with H1, H16, and H22 being the most frequent, represented by 91, 18, and 44 specimens, respectively. As shown in Table 1 and Fig. 5, the most common haplotype (H1) is distributed throughout the entire study area, except for the northwesternmost part, where only H22 was found. Aside from the extreme northwest, H22 is also widely distributed around Barcelona stretching to the north but seems absent from the southern half of Catalonia (15 specimens analyzed from three different populations, P18-P19, P43). The third most common haplotype was restricted to two geographically close populations (P6, P53) from the Montnegre i el Corredor Natural Park. While some populations with a relatively high number of sequenced specimens presented a single haplotype (P39, P11, P24, P46), other populations presented higher haplotype diversity with the highest being P36 (six haplotypes), P31 (five haplotypes), P5 (four haplotypes), and P32-P33 and P37 (three haplotypes). Interestingly, and in contrast to the situation in northwestern Catalonia, all populations with three or more haplotypes were restricted to northeastern Catalonia (Fig. 5).

Identification and distribution of Pelophylax in Catalonia

In the present study, we have generated a comprehensive and exhaustive sampling of Pelophylax water frogs across Catalonia. This region serves as a natural connection between the Iberian Peninsula and the rest of Europe due to the lower mountain ranges found in the easternmost side of Pyrenees. Not surprisingly, this area has been a key corridor for the dispersal of several amphibian genera into Europe, such as Salamandra (Gippner et al., 2024), Alytes (Ambu et al., 2023; Dufresnes & Hernandez, 2021; Dufresnes & Martínez-Solano, 2020), or Pelodytes (Dufresnes et al., 2020), but it may have also facilitated the entry of other groups, including Triturus marmoratus (Kazilas et al., 2024) or Discoglossus pictus (Dufresnes et al., 2020), into Iberia. Given its potential role as a natural dispersal route for Pelophylax species and kleptons from France into the Iberian Peninsula, this region provides an ideal setting to investigate the dynamics of water frogs in Iberia (Arano et al., 1995; Arano & Llorente, 1995). Here, we have provided reliable identification for 423 specimens through the implementation of a two-step, cost-effective molecular protocol and showed that, despite previous reports on the presence of P. idibundus and the hybrid P. kl. esculentus in the Iberian Peninsula (Arano et al., 1995; Holsbeek & Jooris, 2010), the only species found in this work were the autochthonous species P. perezi and the hybrid P. kl. grafi of the P-G hybridogenetic system (Fig. 1A).

The usage of a molecular protocol has been vital for the proper identification of P. perezi and P. kl. grafi, as the morphological identification of both species has proven unreliable in adults (incorrect identification was as high as 26% in males and 38% in females), and nearly impossible in juveniles (26% of misidentifications and 46% of undetermined specimens; Table 2). Most identification errors stem from P. kl. grafi being incorrectly identified as P. perezi, making the primary diagnostic character used to differentiate between P. perezi and P. kl. grafi (the distance between the vomerine ridges) useful only in clear cases—either when the vomerine ridges are distinctly separated, as in P. perezi, or when they are clearly in contact, as in P. kl. grafi. This situation contrasts with the L-E system of P. kl. esculentus in which, morphological identification of the klepton is more reliable. For instance, Plötner (2010) was able to assign in the field 90% of adult individuals correctly to either parental (P. lessonae or P. ridibundus) and the hybrid (P. kl. esculentus) based on morphological characters, particularly the size and shape of the metatarsal tubercle. Additionally, the usage of this molecular protocol has allowed for the identification of tadpoles (n = 69), for which there are no morphological diagnostic characters available.

Arano & Llorente (1995) reported the presence of P. kl. grafi in six localities widely distributed across Catalonia. In the present study, we observed an increase in the abundance of P. kl. grafi in those previously reported localities (Table 1) and expanded its known distribution to nearly the entire study area, with occurrences in 37 localities (Table 1). The only exception was the southeastern Catalonian coast, where no P. kl. grafi were recorded despite intensive sampling efforts (Fig. 3 and Table 1). The absence of kleptons in these populations was previously proposed to be explained by a relatively recent entrance of the complex from Southern France into the Iberian Peninsula (Arano & Llorente, 1995). However, the relatively wide distribution of P. kl. grafi already reported 23 years ago (Arano & Llorente, 1995), the high dispersal capability of these frogs (as evidenced by the widely distributed mitochondrial haplotypes; Fig. 5), and the high abundance of P. kl. grafi in population P19 further south (accounting for 90% of all sampled specimens; Fig. 3 and Table 1) suggest that other factors may play a role in the absence of P. kl. grafi from this large area of Catalonia. One such factor could be the lack of waterborne corridors (such as rivers or streams) reaching that region, which may impose a geographic barrier impeding the dispersal of the klepton.

Origin and maintenance of P. kl. grafi

The origin and age of arrival of Pelophylax kl. grafi into the Iberian Peninsula is still unknown. Arano & Llorente (1995) considered two entry paths for the P-G system in the Iberian Peninsula, one across the central Pyrenees following the Segre-Ebro River system and another one across the eastern Pyrenees, following the Llobregat River. In the present study, we have also found populations with P. kl. grafi (P46, P48-P50; Fig. 3) in northwestern Catalonia, suggesting that colonization through the river basins directly connecting with France through the Val d’Aran (Aran Valley) could have also been possible. Further studies evaluating the nuclear P. ridibundus genome in P. kl. grafi are needed to better understand how and when this genome entered the Iberian Peninsula.

The lack of P. ridibundus in Catalonia may give some clues on how many times this particular hybridogenetic system has appeared. In the L-E system, mattings between two P. kl. esculentus almost invariably leads to P. ridibundus tadpoles with developmental difficulties, resulting in death before reaching sexual maturity (Berger, 1967; Graf & Polls-Pelaz, 1989; Vorburger, 2001a; Vorburger, 2001b; Vorburger, 2001c) (Fig. 1B). This is the consequence of the accumulation of recessive deleterious mutations due to the Muller’s ratchet (Muller, 1964) in the hemiclonaly inherited P. ridibundus genome (Holsbeek & Jooris, 2010; Lehtonen et al., 2013). However, as demonstrated by Vorburger (2001b), experimental crossings between P. kl. esculentus pairs with different hemiclonal P. ridibundus genomes produce viable P. ridibundus offspring with no signs of reduced fitness. Moreover, it has also been reported rare cases (3%) of all-female offspring viability between crossings of P. kl. esculentus bearing the same clonal genome (Vorburger, 2001b). In such cases, recombination of the clonal genomes they inherited from the two parental hybrids seems to occur and, after mating with syntopic P. lessonae, produce new hemiclones with novel combinations of alleles, thus providing an escape from the evolutionary dead-end they were trapped in. The lack of P. ridibundus in the present study suggests that hybrids in the P-G system are not able to produce viable offspring, pointing towards a single origin for the hemiclonal P. ridibundus genome transmitted by P. kl. grafi. Multiple origins of the P-G hybridogenetic system would involve more diversity and a higher probability for successful offspring resulting from a crossing between two hybrids.

In populations where hybrids and P. perezi coexist, P. kl. grafi tends to be more abundant than P. perezi (Fig. 3). This result aligns with the reproductive strategy of hybridogenetic complexes, where offspring between hybrids and pure specimens will always result in more hybrids. Additionally, hybrid females are larger than P. perezi females, providing a further reproductive advantage over P. perezi. However, this increased reproductive success of P. kl. grafi poses an evolutionary conundrum that should, theoretically, lead to the collapse of water frogs of the P-G system (Vorburger & Reyer, 2003). If the reproductive advantage of the klepton over the parental species were to drive P. perezi to extinction, it would also result in the demise of the klepton due to the non-viability of hybrid-hybrid crossings. Nevertheless, the persistence of hybrid populations over at least 30 years (Arano et al., 1995) suggests that some form of equilibrium might exist within the system. This ecological stability could be achieved through sexual preferences in mating, similar to those observed in the L-E system (Bove, Milazzo & Barbuti, 2014; Skierska et al., 2023). Empirical and theoretical studies have shown that the stability of L-E populations can be maintained by specific female mate choice, with P. lessonae females preferentially mating with P. lessonae males over P. kl. esculentus (Bove, Milazzo & Barbuti, 2014; Skierska et al., 2023). Moreover, female mate choice can facilitate the diffusion of deleterious mutations within the system (Bove, Milazzo & Barbuti, 2014). However, further experimental studies and long-term monitoring of water frog populations are necessary to determine whether this ecological stability exists or if P. perezi is gradually declining. Future monitoring should especially focus on localities where only P. kl. grafi has been reported to assess the viability of these populations and to better understand how full P. kl. grafi populations can persist.

Since the origin of the klepton remains unknown, it is not possible to consider the P. kl. grafi as either a native taxon or the result of anthropogenic intervention. However, the stability of the P-G system could be significantly threatened by the recent introduction of other Pelophylax species to the system. It has been shown that the introduction of pure P. ridibundus specimens can break the stability of the L-E system leading to the collapse of the entire population (Bove, Milazzo & Barbuti, 2014). Several stable populations of introduced P. ridibundus, P. kurtmuelleri, P. lessonae, and P. kl. esculentus as well as the recently identified klepton PK (i.e., P. perezi x P. kurtmuelleri) have already been reported in southern France, very close to the Spanish border (Fig. 1; Dufresnes et al., 2017; Dufresnes et al., 2024a; Demay et al., 2023). While the Pyrenees seem to be acting as a geographic barrier to the dispersal of these species into the Iberian Peninsula (i.e., no evidence of any of those species has been reported in this study), it is only a matter of time until some of these species reach the Iberian Peninsula, probably leading to the destabilization of the P-G system.

Origin and genetic variability of the mtDNA of P. kl. grafi in Catalonia

According to the reconstructed phylogenies (Fig. 4), all 194 sequenced P. kl. grafi specimens present P. perezi mtDNA. This could indicate that the first hemiclone occurred between a female P. perezi and a male P. ridibundus and therefore the P. kl. grafi never had mtDNA of P. ridibundus or P. lessonae (i.e., a possibility if the origin was from the P. kl. esculentus, which has been shown to have mainly P. lessonae mtDNA; Plötner et al., 2008). However, studies conducted on the L-E system of P. kl. esculentus suggest that, due to size-related behavioral factors, natural primary hybridizations between P. ridibundus and P. lessonae are more likely to occur between large P. ridibundus females and small P. lessonae males, rather than the reciprocal pairing (Berger, 1970; Borkin & Tikhenki, 1979). Another possibility is that mattings occur between P. perezi females and P. kl. grafi males. Such matings produce P. kl. grafi lineages carrying P. perezi mtDNA. In the absence of the parental species P. ridibundus or P. kl. esculentus (which may have P. ridibundus or P. lessonae mtDNA), the transfer of P. perezi mtDNA into the hybridogenetic lineage becomes irreversible. As previously suggested in the L-E system (Plötner et al., 2008), the mtDNA of the local P. perezi populations may offer a selective advantage over P. ridibundus mtDNA in withstanding the often hypoxic and polluted environments where P. perezi and P. kl. grafi populations are commonly found. Similar to P. ridibundus tadpoles introgressed with P. lessonae mtDNA, it is plausible that P. kl. grafi individuals carrying P. perezi mtDNA are better adapted to survive under hypoxic conditions.

Judging by the high levels of mtDNA genetic variability found in P. kl. grafi (22 different haplotypes; Fig. 5), crossings between male P. kl. grafi and female P. perezi must be a frequent phenomenon, facilitating the continuous incorporation of new P. perezi mitochondrial genomes into the hybrids. The large distribution range of certain haplotypes, such as H1 and H22, suggests a high dispersal capability, with gene flow being possible even between very distant localities (Fig. 5). This could be favored by the fact that Pelophylax is a genus of highly aquatic frogs capable of establishing populations in a variety of habitats and tolerating a wide range of environmental conditions, including salinity, contamination, human activities, and eutrophication. The only significant limiting factor for their presence is the permanent absence of water (Arano & Llorente, 1995; Barbadillo et al., 1999). It is noteworthy that some populations have up to six different haplotypes, most likely acquired through in situ crossings within sympatric P. perezi. However, since we did not sequence any P. perezi mtDNA we cannot ensure that the geographic distribution of the genetic variability found in P. kl. grafi (Fig. 5) mirrors the phylogeography of P. perezi in Catalonia.