https://orcid.org/0000-0002-8497-8908
Abstract
The rich genetic and phenotypic diversity of species complexes is best recognized through formal taxonomic naming, but one must first assess the evolutionary history of phylogeographic lineages to identify and delimit candidate taxa. Using genomic markers, mitochondrial DNA barcoding, and morphometric analyses, we examined lineage diversity and distribution in the Iberian endemic frog Rana parvipalmata. We confirmed two deep phylogeographic lineages, one relatively homogenous genetically, found in Asturias and adjacent areas (T2), and one more fragmented and locally genetically impoverished, restricted to Galicia (T1). Analyses of their hybrid zone suggested a shallow transition characterized by far-ranging admixture, which was modelled by a wide geographic cline (~60 km for the genome average) and no obvious barrier loci (i.e. loci showing disproportionally restricted introgression). The relatively young T1 and T2 have thus remained reproductively compatible, which argues against their distinction as biological species, and we accordingly describe T2 as a new subspecies, Rana parvipalmata asturiensis ssp. nov. Remarkably, we highlight striking discordances between mitochondrial and nuclear distributions across their hybrid zones, as well as between their genetic and morphological differentiation. Our study illustrates how genomic-based phylogeographic frameworks can help decipher the high genetic and phenotypic variation of species complexes and substantiate the taxonomic assessment of candidate lineages.
INTRODUCTION
Delimiting and scientifically naming the diversity of species complexes, i.e. widely-distributed species composed of multiple phylogeographic lineages, is necessary to raise social awareness towards evolutionarily-unique populations and promote their legal protection (Freudenstein et al. 2017, Pollock et al. 2020, Liu et al. 2022). This is particularly relevant when threats differ between geographically-restricted lineages, but these lineages are all treated as a single widespread taxon in wildlife conservation assessments, hence appearing unthreatened at the global scale (Bland et al. 2015). To date however, delimiting geographic variation within species is complicated by methodological limits to detect lineages (Degnan and Rosenberg 2009), and controversies regarding the way to split these lineages into multiple taxa (Zachos 2015, 2016, de Queiroz 2020, Hillis 2020, 2022, Sukumaran et al. 2021, Dufresnes et al. 2023).
One of the most important issues in translating phylogeography into taxonomic classification is the over-reliance on molecular datasets built mostly from mitochondrial DNA (mtDNA) sequences (Rubinoff and Holland 2005, Chan et al. 2022, Dufresnes and Jablonski 2022). It has become clear that mtDNA may show discordant phylogeographic patterns compared to the rest of the genome (i.e. cyto-nuclear discordance), and may thus misinform on the evolutionary history of taxa and their geographic distributions (e.g. Cairns et al. 2021, Dufresnes et al. 2021a, Kimball et al. 2021). Moreover, forming a robust monophyletic group in any mtDNA or nuclear DNA tree does not always imply speciation in the biological sense: structured populations (e.g. diverging in allo-parapatry) are expected to form robust lineages at the genomic scale, even if they remained reproductively compatible (Sukumaran and Knowles 2017, Campillo et al. 2020, Dufresnes et al. 2020a). Describing any evolutionary lineage as a ‘species’ may contribute to the emerging problem of taxonomic inflation—poisoning biodiversity assessments and red lists (Zachos 2013). In a similar fashion, phenotypic differentiation, which is often used to complement phylogenetic evidence for new species (e.g. ‘integrative taxonomy’, Padial et al. 2010), must be regarded with caution in the case of phylogeographic lineages. Striking ecological and behavioural innovations may arise between closely-related populations due to processes (local adaptation, plasticity) other than evolutionary divergence and reproductive isolation (e.g. reproductive strategies, Dufresnes et al. 2022; mating call ‘accents’, Weaver et al. 2020), and reciprocally, the same phenotypes may remain conserved even across species boundaries (i.e. ‘cryptic species’, Bickford et al. 2007).
Genomic phylogeography (also called ‘next-generation phylogeography’) can help decipher and delimit the diversity of species complexes. Genomic methods like restriction-site associated DNA-sequencing (RAD-seq) or targeted sequence capture allow the genotyping of thousands of nuclear markers distributed all over the genome to identify candidate evolutionary lineages for taxonomic revision, and test whether these correspond to previously identified mitochondrial lineages and eco-morphotypes (Cruaud et al. 2014, Edwards et al. 2016). Providing a dedicated sampling of parapatric populations, such analyses can further inform on the specific vs. subspecific rank of candidate taxa through the study of hybrid zones (Dufresnes et al. 2020a, 2021b, Chambers et al. 2022). Narrow hybrid zones, where gene flow between the contacting lineages is weak and geographically restricted, may reflect partly unsuccessful hybridization and introgression due to genomic incompatibilities responsible for post-zygotic isolation, in which case these lineages match the definition of biological species (see e.g. Seehausen et al. 2008). Alternatively, the pervasive admixture typically observed across shallow hybrid zones implies that the lineages have remained genetically compatible and should rather be considered as subspecies (Dufresnes et al. 2021b). Genomic phylogeography can thus reliably evaluate both phylogenetic (genomic divergence) and biological (amount of contemporary gene flow) criteria to delimit species complexes in a way that will convince a broad community of taxonomists. Moreover, because hybrid zone width and genetic divergence covary (Singhal and Moritz 2013, Dufresnes et al. 2021b), one could eventually predict in a probabilistic framework (Kollár et al. 2022) whether lineages represent biological species based on their level of molecular differentiation (e.g. Dufresnes and Litvinchuk 2022). Enforcing the hybrid zone approach of species delimitation in additional pairs of candidate taxa is needed to fine-tune this relationship and make it usable across a wide array of organisms.
European common frogs (Rana cf. temporaria) feature tremendous phenotypic and genetic variability that has puzzled zoologists for centuries. Widespread from northern Spain to northwestern Asia, where they have colonized various temperate habitats from sea level to mountain tops, these cold-tolerant anurans exhibit strong variation in size, shape, and coloration, which has led to the descriptions of dozens of forms that were later synonymized or still occasionally considered as subspecies (Frost 2023). Phylogeographic analyses have identified no less than six Plio–Pleistocene mitochondrial lineages in Europe (labelled T1–T6 in Vences et al. 2013, 2017) yet with a general lack of correspondence to morphotypes (e.g. Veith et al. 2002), and some important discordances with nuclear markers (Veith et al. 2002, 2012, Vences et al. 2013, Dufresnes et al. 2020b).
The northwestern Iberian populations are particularly diverse. Four parapatric mitochondrial lineages segregate along the Spanish Atlantic coast, labelled T1 (Galicia), T2 (Asturias), T6 (Cantabria), and T4 (Basque Country, Pyrenees) in previous works (Vences et al. 2013, 2017, Dufresnes et al. 2020b). Using RAD-seq, Dufresnes et al. (2020b) showed that the sister lineages T1 and T2 were also significantly differentiated at the nuclear level, but that populations bearing the unrelated T6 and T4 matrilines were one and the same evolutionary lineage. Furthermore, the old nuclear divergence (c. 4 Mya) and the narrow hybrid zone (~25 km wide) between the T2 and T6 populations led to the recognition of their respective clades (T1/T2 vs. T4/T6) as different species. Galician and Asturian T1/T2 populations have subsequently been called Rana parvipalmataLópez-Seoane, 1885, and populations from the rest of Europe (T3–T6) remain provisionally assigned to Rana temporaria Linnaeus, 1758 (Speybroeck et al. 2020).
Despite the recent species distinction of R. parvipalmata, the diversity of Iberian common frogs may still remain under-appreciated. The deep phylogeographic structure recovered between the Galician (T1) and Asturian (T2) populations suggests that at least two candidate taxa might be hidden under this name. While their estimated divergence (Early Pleistocene, ~2.5 Mya) is on the subspecies rather than the species side of the anuran grey zone of speciation (Dufresnes et al. 2021b), a thorough taxonomic assessment requires locating their geographical transition zone and quantifying genetic introgression. In addition, these candidate taxa may not be so cryptic. The name ‘parvipalmata’ initially referred to a Galician morphotype characterized by reduced foot webbing (López-Seoane 1885) compared to other Iberian populations, together with call features and a smaller body size (Galán 1989, Vences 1992). Whether the phylogeographic structure associates to morphological differentiation thus becomes a relevant taxonomic and evolutionary question.
Finally, assessing the genetic diversity and structure throughout the species’ distribution range appears timely to appreciate conservation needs independently between phylogeographic lineages. In northwestern Spain, common frogs inhabit humid Atlantic forests and may be sensitive to global warming (Enriquez-Urzelai et al. 2019a, b, 2020). Yet, immediate threats could differ between the T1 Galician populations, which mostly occupy lowland humanized landscapes that become increasingly exploited for eucalypt plantations, known to be harmful to amphibians (Burraco et al. 2018, Iglesias-Carrasco et al. 2022), and T2, which spans a wider altitudinal gradient that could help to buffer the effect of climate change (Enriquez-Urzelai et al. 2020). Moreover, the southernmost edge populations sustain in isolated patches that could face a higher risk of extinction (Galán et al. 2010).
Following up on our previous work (Dufresnes et al. 2020b), we here report on a genomic (RAD-seq) phylogeography of R. parvipalmata, with a specific focus on the transition between the Galician T1 and Asturian T2 lineages. We were primarily interested in assessing their taxonomic rank (species vs. subspecies) based on the hybrid zone approach of species delimitation. Given their relatively recent divergence and the apparent lack of spatial barriers to dispersal, we predict that T1 and T2 admix over relatively large distances, in which case they may be distinguished as subspecies (Dufresnes et al. 2021b). Second, because discordances between nuclear and mitochondrial phylogeographic patterns have been reported in the complex, notably in Iberia (Veith et al. 2012, Dufresnes et al. 2020b), special attention was given to comparing the distributions inferred from RAD-seq and mtDNA loci. Third, we assessed whether the two R. parvipalmata lineages (T1 and T2) differ morphologically, also with respect to the genetically more distant (but geographically proximate) Cantabrian populations assigned to R. temporaria (featuring T6 mtDNA). We predict that T1 and T2 are more similar to each other than to T6 if the morphology follows the phylogeny and represents a solid cue for taxonomic distinction. Alternatively, only the Galician T1 may be differentiated if the remarkable ‘parvipalmata’ morphotype reflects local adaptation or phenotypic plasticity. Building on our results, we provide a taxonomic description to recognize the T2 populations of R. parvipalmata as a new taxon for the European herpetofauna.
MATERIAL AND METHODS
RAD-seq genotyping and analyses
Tissue samples were obtained from 67 R. parvipalmata specimens sampled in 35 localities (Supporting Information, Table S1). Samples included non-invasive buccal swabs for adult frogs (N = 11)—which provide adequate DNA for genomic analyses (Ambu and Dufresnes 2023)—or young larvae fixed in 96% ethanol (N = 56). DNA was extracted using a Qiagen DNeasy Blood & Tissue Kit and run on an 1.5% agarose gel to check for DNA yield and integrity.
A genomic library was prepared based on the double-digest RAD (ddRAD) protocol with SbfI and MseI detailed in Brelsford et al. (2016), modified as in Ambu et al. (2023). The modifications consisted of two replicate PCR amplifications of the ligated fragments (instead of four), and the size selection (400–500 bp) being carried on a PippinPrep (Sage Science) instrument (rather than on an agarose gel), keeping an overall design that is not expected to be sensitive to PCR duplicate errors (Rochette et al. 2023). The final library was purified and sequenced (paired-end) on an Illumina NextSeq 550 using a 2 × 150 bp High-Output Kit. Raw reads were demultiplexed with STACKS 2.59 (Catchen et al. 2013) and trimmed to 100 bp. We then used the denovo_map.pl pipeline of STACKS for constructing RAD loci, assembly, and cataloguing, with default ‐m, ‐n, and ‐M values, which provided a good balance between data quality and quantity (Paris et al. 2017). The STACKS catalogue contained 384 892 loci, with an average coverage of 15.8 reads (9.0–25.6, SD = 2.8). Loci were then filtered using the populations module for downstream analyses.
Patterns of phylogeographic structure were explored with two different datasets based on all the RAD loci present in all individuals. First, we exported a 144 916 bp alignment concatenating 812 RAD loci, which we analysed with a phylogenetic network in SplitsTree 4.18.3 (Huson and Bryant 2006), using the uncorrected P distances. Second, we exported a genotype matrix of 812 presumably unlinked SNPs by randomly keeping a single SNP (-wrs) from each RAD locus. We conducted a Principal Component Analysis (PCA) on individual allele frequency with the R packages ade4 and adegenet (Jombart 2008). Individual genotypes were assigned to genetic groups using the clustering algorithm of STRUCTURE 2.3 (Pritchard et al. 2000), testing five replicates each from K = 1–8 and using the ancestry model with correlated allele frequencies. Run statistics were visualized in STRUCTURE HARVESTER (Earl and vonHoldt 2012). We further estimated expected heterozygosity with the R package hierfstat (Goudet 2005).
To assess the structure of the T1/T2 hybrid zone, we used the R package hzar (Derryberry et al. 2014) to fit clines to the hybrid index of 30 populations located along a west-east geographic transect (locs 3–22, 24–32, 34). The STRUCTURE Q estimate obtained from the K = 2 runs, which correspond to the genome average population ancestry between T1 and T2, was considered as the hybrid index. Clines are defined by a centre (c) and a width (w), and may be flanked by introgression tails of various length (δ) and slope (τ). We tested cline models from two (no tails) to six (different tails on each side) parameters, fixing the minimum and maximum frequencies to 0 and 1, and selected the one with the lowest Akaike information criterion (AIC) score. Finally, to assess the heterogeneity of introgression across the genome, we fitted clines individually to the allele frequencies of lineage-diagnostic SNPs, i.e. those fixed between the edge populations 3–7 (pure T1) and 32–34 (pure T2).
Mitochondrial genotyping and analyses
The mitochondrial lineage appartenances of 168 samples from 42 localities were identified by amplifying and Sanger-sequencing ~550 bp of the cytochrome b gene (cyt-b) with the primers and protocol from Dufresnes et al. (2020b). Sequences were trimmed and aligned in SeaView 5 (Gouy et al. 2021) alongside previously published R. parvipalmata sequences taken from Galán et al. (2010), Vences et al. (2013), and Dufresnes et al. (2020b). The final alignment contained 400 sequences (501 bp) from 69 localities (Supporting Information, Table S2), and was analysed by a phylogenetic network with SplitsTree. For taxonomic diagnosis, we used MolD (Fedosov et al. 2022) as implemented in iTaxoTools (Vences et al. 2021) to determine diagnostic positions in the cyt-b alignment, using the full gene sequence of R. temporaria as a reference (GenBank accession NC_042226).
Morphometric analyses
Because anuran morphology is often sex-specific and males were more frequently caught in the field, we limited our morphometric comparisons to adult males. First, we compared the body size (snout-vent length, SVL) of 1947 individuals measured during the 1997–2021 period in 35 populations of R. parvipalmata (T1 and T2) and nearby R. temporaria from Cantabria (T6) (Supporting Information, Table S3). The effect of lineage was assessed by a mixed model using population as a random factor, as well as altitude as a covariate. The relationship between altitude and body size was further assessed by a linear regression on population-averaged log-transformed data. Pairwise differences between lineages were visualized by a neighbour-joining (NJ) tree of the distance matrix. For the purpose of diagnosis, we also tested differences between pairs of lineages with a Tukey Honestly Significant Difference (HSD) test.
Second, for a subset of 82 individuals representative of the three lineages (Supporting Information, Table S3), 17 morphological characters summarizing body shape were measured by only one of us (C. Dufresnes) with a digital caliper (precision to the nearest 0.1 mm), as follows: space between the dorsolateral ridges (DLT); head width (HW); eye diameter (ED); inter-eye distance (EE); eye-nostril distance (EN); inter-nostril distance (NN); nostril-snout distance (NS); tympanum diameter (TD); tympanum-eye distance (TE); femur length (FMR); tibia length (TIB); tarsus length (TAR); unwebbed portion of the first toe (WT1); unwebbed portion of the second toe (WT2); unwebbed portion of the third toe (WT3); height of the metatarsal tubercle (MTH); length of the metatarsal tubercle (MTL). A schematic drawing and additional details on these measurements are provided in Supporting Information, Fig. S1. To remove the effects of body size, we applied the allometric adjustment approach of Chan and Grismer (2021) implemented in the R package GroupStruct, which allows to account for structure between taxa and populations. Data from each lineage was independently scaled with the common within-group pooling method that considers each population separately (‘population1’ option), as recommended when body size differs between and within lineages (see Results). The dataset was then analysed by a PCA, and multivariate pairwise Euclidian distances between lineages were visualized by a neighbor joining (NJ) tree. To test for global differences between lineages, we performed two independent multivariate analyses of covariance (MANCOVA): one on the size-adjusted dataset, using population as a cofactor, and one on the raw measurements, using population and size as cofactors. We assessed statistical significance through a permutation approach by randomizing lineage and population assignments 10 000 times. Finally, to search for diagnosable characters between lineages, we computed simple body ratios (raw measurements divided by SVL) and tested for statistical differences by Tukey tests, adjusting significance thresholds for multiple testing. All analyses were performed in R 4.1.3 (R Core Team 2022).
RESULTS
Phylogeographic structure and diversity
Population genetics (812 SNPs) and phylogenetics (~145 kb) concurred that R. parvipalmata forms two divergent phylogeographic lineages, namely T1 in Galicia and T2 in Asturias (Fig. 1). They are discriminated by the first axis of the PCA (Fig. 1C), the phylogenetic network (Fig. 1D), and the STRUCTURE analysis (Fig. 2B). In the latter, K = 2 best summarizes the genetic variation—runs with K > 2 did not depict any additional spatially-explicit structure, with no improvement in the runs’ statistics (Supporting Information, Fig. S2).
Phylogeographic structure and diversity in R. parvipalmata based on RAD-seq analyses. The map (A) shows the location and numbering of sampled localities; the yellow shades outline the presumed distribution of the species. Barplots (B) and the PCA (C) display within-population expected heterozygosity He and between-sample genetic variation, both based on 812 unlinked SNPs. (D) The phylogenetic network shows genetic relationships among samples, based on ~155 kb sequences. Colour-coding follows a red-orange gradient that corresponds to the average genomic ancestries of samples/populations between the Galician R. p. parvipalmata T1 (red) and the Asturian R. p. asturiensis T2 (orange) lineages (Fig. 2). Triangles and squares highlight the isolated southernmost and westernmost samples/populations. The star shows the type locality and the holotype of R. p. asturiensis.
Characterization of the hybrid zone between the Galician R. p. parvipalmata T1 (red) and the Asturian R. p. asturiensis T2 (orange) lineages. The map (A) and the STRUCTURE barplot (B) show population and individual ancestries, based on 812 unlinked SNPs. The cline plot (C) models the geographic transition by fitting a sigmoid function to the genome-average ancestry data (STRUCTURE Q) and to 17 unlinked individual SNPs fixed between the lineages.
Genetic distances and diversity are relatively homogenous in the Asturian T2, but not in the Galician T1. For T1, southern (locs 1–2, triangles in Fig. 1) and western populations (locs 3–8, squares) show lower levels of expected heterozygosity compared to northern populations (Fig. 1B). They also feature some genetic isolation, as seen from the long branches on the phylogenetic network (Fig. 1D) and their contribution to the second PCA axis (Fig. 1C).
The analyses suggested extensive contemporary gene flow between the T1 and T2 lineages. Many geographically intermediate samples received intermediate scores on the first PCA axis (Fig. 1C) and branched in the middle of the phylogenetic network (Fig. 1D). STRUCTURE accordingly retrieved mixed ancestry in these samples (Fig. 2), especially between the 80 km that separate loc. 13 (90% T1) from loc. 28 (88% T2). Using STRUCTURE’s Q as a genome-average hybrid index, the transition was best modelled by a two-parameters geographic cline (Fig. 2C), with a width w of 61.4 km (95% confidence interval [CI]: 37.1–94.0) and a centre c of 202.5 km (95% CI: 190.4–214.5), located within the 14 km gap that separates loc. 19 and loc. 20. Among the 812 SNPs, only 17 were fixed between the edge populations of the two lineages and could be fitted locus-specific clines (Fig. 2C). Their widths varied from 42.3 to 264.3 km, with a median of 93.7 km; the centres varied from 149.8 to 224.1 km, with a median of 204.5 km (Supporting Information, Table S4).
The mitochondrial phylogeography partly contrasts with the nuclear genomic structure. The cyt-b network features three haplogroups (Fig. 3): T1a in Galicia, which corresponds to the Galician T1 nuclear cluster; T2 in Central and Eastern Asturias, which corresponds to the Asturian T2 nuclear cluster; and T1b, which is genetically closer to T1a, but segregates in Western Asturian populations that are either admixed or mostly assigned to the T2 cluster in the nuclear analysis (e.g. loc. 23). Accordingly, the geographic transition between the T1 and T2 nuclear clusters corresponds to the transition between the T1a and T1b mtDNA, while the transition between T1b and T2 mtDNA is located within nuclearly pure T2 populations. The MolD analysis flagged three cyt-b positions that always differentiate the T2 frogs (carrying either lineage T2 or T1b) from any other lineage of the R. cf. temporaria complex.
Mitochondrial phylogeography of R. parvipalmata. Dots on the map (A) show the distribution of the three main cyt-b haplogroups identified on the phylogenetic network (B). The respective ranges of the two nuclear lineages are illustrated by shaded backgrounds (T1: red; T2: orange; dashed line: putative centre of the hybrid zone). Note how the T1b matrilines (pink dots) are found across T2 populations from western Asturias (orange area).
Morphometric comparisons
Comparisons of body size (SVL) among 1947 adult males suggest that R. parvipalmata T2 frogs (68.4 ± 5.7 mm) were mildly smaller than R. temporaria T6 frogs (72.8 ± 4.6 mm), but clearly larger than the more closely related R. parvipalmata T1 frogs (59.4 ± 5.2 mm), as reflected by the NJ tree of pairwise size distances (Fig. 4). All three pairwise comparisons were significant in the Tukey test (P < 0.001). Lineage was a significant predictor of SVL in the mixed model (Χ2 = 10.1, P = 0.006, d.f. = 2), so was altitude (Χ2 = 3.9, P = 0.05, d.f. = 1), but not their interaction (Χ2 = 1.0, P = 0.60, d.f. = 2). Accordingly, SVL significantly increased with elevation (F = 6.1, R2 = 0.16; P = 0.02; Supporting Information, Fig. S3).
Comparison of body size in common frogs from northwestern Spain. (A) Violin plots showing variation in SVL among 1947 males of R. p. parvipalmata (T1, red), R. p. asturiensis (T2, orange), and R. temporaria from nearby Cantabria (T6, blue). (B) Neighbour-joining tree of pairwise size differences among these three lineages. Photo: T1 frog from La Garganta (credit: J. Ambu).
Based on 17 characters measured in 82 individuals, T1 and T2 R. parvipalmata frogs also differed in their general morphology, with non-overlapping scores on the first axis of the PCA (Fig. 5A). As for SVL, the T2 individuals were more similar to those sampled in Cantabria (R. temporaria, mtDNA T6) (Fig. 5B). Accordingly, the effect of lineage was statistically significant in both MANCOVAs, i.e. on adjusted measurements (F = 11.4, P <<< 0.0001, d.f. = 2), and on raw measurements with body size as a covariate (F = 106.2, P <<< 0.0001, d.f. = 2).
Comparison of body shape in common frogs from northwestern Spain. (A) PCA on 17 morphological characters (after allometric size correction) among 82 males of R. p. parvipalmata (T1, red circles), R. p. asturiensis (T2, orange circles), and R. temporaria from nearby Cantabria (T6, blue triangles); scatterplots show coordinates and circles show variable contributions on PC1–2 (top) and PC3–4 (bottom). (B) Neighbour-joining tree of multivariate pairwise differences among lineages.
Exploring diagnostic characters by comparing body ratios revealed significant differences (after Bonferroni corrections) between pairs of lineages for nine variables, but the data distributions always overlapped (Supporting Information, Fig. S4). The main differences pertained to the generally larger metatarsal tubercle of lineage T2 (variables MTH and MTL), as well as the lesser webbing (variables WT1, WT2, WT3), longer femur (variable FMR) and shorter tibia (variable TIB) of lineage T1 (Supporting Information, Fig. S4).
DISCUSSION
Previous investigations of Spanish common frogs had documented two distinct evolutionary lineages of R. parvipalmata, namely T1 in Galicia and T2 in Asturias, which were dated to the Early Pleistocene in multilocus nuclear (~2.5 Mya) and mitochondrial (~2 Mya) time trees (Vences et al. 2017, Dufresnes et al. 2020b). Combining genomic and morphometric analyses, we mapped the diversity and distribution of these lineages, analysed their hybrid zones, and examined whether they differed in morphology. As developed in the following, the case of R. parvipalmata provides insights into the delimitation of phylogeographic diversity in species complexes.
Our main result pertains to the hybrid zone between the T1 and T2 lineages, which we found to be rather smooth, featuring nuclear introgression over a hundred kilometres across the hilly landscapes of western Asturias (Fig. 2). The region seems to act as a meeting (and admixing) point for amphibian phylogeographic lineages associated with the humid Atlantic forests of northwestern Iberia, such as the salamander subspecies Salamandra salamandra gallaica and Salamandra salamandra bernardezi (Velo-Antón et al. 2021), the midwife toad subspecies Alytes obstetricans boscai and Alytes obstetricans pertinax (Dufresnes and Hernandez 2021, Ambu et al. 2023), and probably between evolutionary significant units of the palmate newt Lissotriton helveticus (Recuero and García-París 2011). The western edge of the Cantabrian mountains can thus be regarded as a local hotspot of intraspecific amphibian diversity.
Without spatial barriers to dispersal, the extent of admixture should reflect the degree of genetic compatibility between contacting lineages, and thus their level of intrinsic post-zygotic reproductive isolation (Barton and Gale 1993). To this end, sigmoid clines are useful to model the transition, and cline shape, summarized by the width parameter w, informs on how constrained genetic introgression is. In frog and toad hybrid zones, partly incompatible lineages feature a mosaic of barrier and neutral loci—loci that barely (w close to zero) and freely (larger w) introgress—for a genome average cline width that does not exceed 30 km (Dufresnes et al. 2021b). This is what we previously reported between R. parvipalmata (T2) and R. temporaria (T6) in Eastern Asturias, where the steep hybrid zone (w ~ 25 km) involved multiple genomic regions impermeable to gene flow (Dufresnes et al. 2021b). Providing some knowledge of the species’ dispersal rate σ, it is possible to compute intrinsic coefficient of selection against hybrids s* from w (Barton and Gale 1993), which Dufresnes et al. (2020b) estimated to s* ~ 0.07 for the T2/T6 interspecific hybrid zone. Using the same demographic parameters, here the much larger width for the T1/T2 genome average (w ~ 60 km) translates into s* ~ 0.01, and even the width of the least introgressed lineage-diagnostic locus remains above 40 km (Fig. 2C). The entire genome, at least the portions covered by our RAD-seq markers, thus seems to admix freely without obvious costs to hybrids, and pervasive introgression is consequently not constrained to the area of initial contact. There is thus little doubt that the T1 and T2 lineages can be considered the same biological species (= R. parvipalmata), and a subspecies rank would be most appropriate for their taxonomic distinction.
This conclusion verifies our initial hypothesis that the T1 and T2 lineages are conspecific based on their young evolutionary age. There is ample evidence that genomic divergence and reproductive isolation evolve in concert in anuran amphibians, as seen from the relationship between divergence time and hybrid zone steepness (Dufresnes et al. 2021b), or between genetic distance and offspring fitness in experimental crosses (Malone and Fontenot 2008). The trend is tangible in common frogs: the older R. parvipalmata/temporaria species pair (~4 Myr) admixes significantly less than the younger R. parvipalmata T1/T2 (~2.5 Myr). Predicting taxonomic ranks based on genetic divergence may thus come in handy when hybrid zones are not available, providing that inferences are calibrated from closely related lineages that share similar ecologies and dispersal capacities, and were assessed using the same molecular methods and analyses (Dufresnes et al. 2023).
Given that cyto-nuclear discordances are a common phenomenon in amphibians (Dufresnes and Jablonski 2022), we stress that genome-wide markers should become a gold standard for such analyses. In R. parvipalmata, most frogs sampled on the T2 side of the hybrid zone feature mtDNA haplotypes (T1b) related to those found in true T1 populations (T1a) (Fig. 3). Climatic conditions in Asturias were alternatively suitable and unsuitable for the species during glacial and interglacial stages (fig. 5 in Dufresnes et al. 2020b), thus the T1/T2 hybrid zone probably formed multiple times, offering demographic conditions that favour discordances. Here, the fact that T1a and T1b consist of distinctively diverged haplogroups pleads for two alternative scenarios. On the one hand, a past mitochondrial capture across Western Asturias during a previous episode of secondary contact would have led to the replacement of the original T2 mtDNA by the T1 ancestral mitogroup. On the other hand, the T1a/T1b divergence could reflect historical genetic structure within lineage T1 (whose range would have extended further east), and the T1b populations would have been assimilated by T2 through a westward movement of the hybrid zone. Hybrid zone movement is known to shift mitochondrial and nuclear transitions (Wielstra 2019), including the fixation of foreign haplotypes (here, T1b) at the front of the expanding wave (here T2) (Currat et al. 2008). In both cases, this situation illustrates how mtDNA can make a poor predictor of phylogeographic transitions (see e.g. Dufresnes et al. 2021a), and adds to the list of striking cyto-nuclear discordances found in Spanish common frogs (Dufresnes et al. 2020b). A genomic phylogeography extending to the entire range of the complex, notably to incorporate additional mtDNA lineages found in the Pyrenees (T3) and Eastern Europe (T5), may thus lead us to reconsider what we know of its evolution, diversity, and taxonomy.
Another striking discordance is the mismatch between genetic divergence and morphological differentiation among the northwestern Iberian lineages. Various evolutionary (local adaptation) or environmentally driven factors (phenotypic plasticity) could have contributed to conserving a similar morphotype across species boundaries in the distantly related T2 and T6 populations, whereas body shape and size diverged in the young T1 lineage (see also Galán 1989, Vences 1992). Among potential factors, elevation was a marginally significant predictor of SVL, and the smaller size of T1 frogs could be partly due to their occurrences at generally lower altitudes compared to the other Cantabrian lineages. Interestingly, the positive relationship retrieved between size and altitude is reminiscent of Bergmann’s rule (Supporting Information, Fig. S3), which remains a conundrum in amphibians due to their ectothermic physiology (Ashton 2002, Adams and Church 2008, Slavenko and Meiri 2015). Whatever the proximate causes, these morphological patterns have complicated the taxonomy of populations. The distinct T1 morphology corresponds to the Galician form described as R. p. parvipalmata (López-Seoane 1885), and the T2 populations were initially regrouped with R. t. temporaria as per their external resemblance, causing that latter taxon to be paraphyletic. Our study thus supports that phylogeographic lineages should be primarily delimited based on multilocus genetic (ideally genomic) evidence of divergence and reproductive isolation, which alone appears more reliable than integrating several lines of traditional evidence like genetic distances at mitochondrial barcoding genes, or singular phenotypic cues—unless the latter are explicitly shown to trigger reproductive isolation.
Weighting both phylogenetic (e.g. monophyly and divergence) and reproductive (e.g. hybridization and introgression) criteria to delimit species complexes is gaining momentum. In the herpetofauna of Europe, enforcing the hybrid zone approach of species delimitation (Dufresnes et al. 2020a, Chambers et al. 2022) has resulted in numerous splits and lumping in the latest species list (Speybroeck et al. 2020). Efforts need to be made, however, to acknowledge conspecific phylogeographic lineages beyond the technical concept of ‘evolutionary significant units’, which can only bring a limited spotlight for conservation and biodiversity stakeholders (Pollock et al. 2020). The concept of ‘phylogeographic subspecies’, i.e. monophyletic yet reproductively compatible evolutionary lineages that may (O’Brien and Mayr 1991) or may not differ phenotypically (Kindler and Fritz 2018, Dufresnes et al. 2023), appears timely to fill this role without contributing to species-level taxonomic inflation (see also Hawlitschek et al. 2012, de Queiroz 2020, Hillis 2020, 2022, Scherz et al. 2022). It is also highly relevant when different lineages of the same species are experiencing contrasting conservation situations, so each is given the appropriate attention rather than being treated as a single entity (e.g. Miralles et al. 2017). This might be the case in R. parvipalmata, where the structured western and southern T1 populations are genetically impoverished (Fig. 1), whereas the populations found along the Cantabrian mountains (T2) retained higher levels of neutral genetic diversity (see also Dufresnes et al. 2020b)—apart from a few high-elevation sites (e.g. Picos de Europa, Choda 2014). Management efforts should ideally be population-specific, especially when functional diversity fluctuates across the populations of the same taxon/lineage due to local adaptation, such as along altitudinal gradients (Cano et al. 2017).
In this respect, here we delimite the two lineages of R. parvipalmata as distinct subspecies that differ by evolutionary and phenotypic divergences: lineage T1, which corresponds to R. p. parvipalmata, and lineage T2, which we formally describe below as a new subspecies.
Description of Rana parvipalmata asturiensis ssp. nov.
(Fig. 6)
ZooBank registration:
urn:lsid:zoobank.org:act:89BF6DB4-36E4-4B3F-9735-D9D5D6F6D699.
Identity:
Both lineages T1 and T2 are currently assigned to the species Rana parvipalmata, a taxon originally described as Rana temporaria parvipalmataLópez-Seoane, 1885. The original description (López-Seoane 1885) reported this taxon from ‘the north-west of Spain, equally common from the level of the sea to an elevation of 1422 feet’, with the additional mention ‘Near Pontevedra, Tui and Ferrol, it occurs conjointly with Rana iberica’—all these explicitly mentioned localities are situated on the Galician coast. Although measurements of three specimens were given, no locality information for these was published and no type material is known to exist (Frost 2023). Salvador (1974) proposed to restrict the type locality to ‘La Coruña’ but did not designate a lectotype or a neotype. Given the localities and the diagnostic characters purported in the original description (parvipalmata means ‘poorly webbed’), it is unambiguous that the nomen parvipalmata is to be assigned to lineage T1. For lineage T2, no taxonomic name is available (see e.g. Frost 2023), and we therefore here describe it as Rana parvipalmata asturiensis ssp. nov.
Diagnosis:
Rana parvipalmata asturiensis is a member of the R. temporaria complex. It corresponds to phylogeographic lineage T2, distributed in Asturias and northern León (northwestern Spain), which has been recently attributed to the distinct species R. parvipalmata based on multilocus mitochondrial and nuclear DNA analyses (Dufresnes et al. 2020b). Its closest relative is Rana parvipalmata parvipalmata, corresponding to phylogeographic lineage T1 distributed in nearby Galicia, from which it differs by 1.2% of divergence at mitochondrial genes (~4.4 kb spanning six genes), including 4.7% at cyt-b and 0.9% at 16S, and by 0.07% at nuclear RAD-seq markers conserved across distant Rana species (all data from Vences et al. 2017 and Dufresnes et al. 2020b). In the western parts of its range, R. p. asturiensis features mtDNA haplotypes derived from R. p. parvipalmata (cyt-b mitogroup T1b), following hybridization. As identified by MolD, R. p. asturiensis can be robustly identified by a combination of diagnostic nucleotides: ‘T’ in the site 666, ‘C’ in the site 714, ‘G’ in the site 732 for specimens carrying T1b mtDNA: ‘T’ in the site 627, ‘A’ in the site 715, ‘T’ in the site 825 for specimens carrying T2 mtDNA (positions relative to the full cyt-b gene of Rana temporaria).
Externally, the new taxon generally resembles other common frogs of the R. temporaria complex. In western ranges, R. p. asturiensis can be confounded with R. p. parvipalmata, from which most individuals can be diagnosed by a combination of morphological characters (Fig. 5; Supporting Information, Fig. S4), notably body size (larger in R. p. asturiensis, up to 80 mm), feet webbing (more extensive in R. p. asturiensis), femur and tibia lengths (respectively shorter and longer in R. p. asturiensis), and the metatarsal tubercle (larger in R. p. asturiensis). These differences alone are not fully diagnostic, however, and field identification should take into account distribution, and ideally rely on genetic barcoding. The mating call and tadpole of R. p. asturiensis have not been studied. Vences (1992) noted that Galician frogs (R. p. parvipalmata) featured a lower number of pulses per call compared to populations from Central Europe (R. temporaria), but whether this pattern holds across the entire range of R. parvipalmata (i.e. including the new subspecies) is unknown. In eastern ranges, R. p. asturiensis can be confounded with the distinct species R. temporaria, from which it remains even more similar phenotypically. Finally, the new taxon shares its range (but rarely its habitats) with the Iberian frog Rana iberica but can be reliably distinguished based on a number of diagnostic criteria pertaining to size, shape, and coloration (Dufresnes 2019).
The taxonomic status of R. p. asturiensis is justified by the genuine evolutionary divergence and morphological differentiation from other R. parvipalmata populations, and the subspecies rank is justified by the recent Plio–Pleistocene origin, the weak sequence divergence at barcoding genes, and most importantly the extensive admixture with R. p. parvipalmata, suggestive of reproductive compatibility.
Holotype:
MNCN 50655, adult male found by J. Ambu and C. Dufresnes on 14 March 2021 on a forest track c. 500 m northeast of Alto de la Fumarea, Picu Fario area, Asturias, Spain (43.4306403°, -5.5710408°, 549 m a.s.l.), and deposited at the Museo Nacional de Ciencias Naturales (MNCN) in Madrid, Spain. The specimen is depicted both live and post-mortem in Fig. 6. Its identity was confirmed by mitochondrial (cyt-b) and nuclear (RAD-seq) analyses (Figs 1–3).
The holotype of Rana parvipalmata asturiensis (MNCN 50655), depicted live (top) and post-mortem before ethanol preservation (bottom). Credit: C. Dufresnes.
Description of the holotype:
Body size (SVL) of 64 mm for a wet weight of 24.1 g; head 16 mm wide: large protruding eyes (8 mm of diameter) spaced by 12 mm, with horizontal pupil; tympanum disk very distinctive and large (5 mm), located 3 mm behind the eyes; nostrils separated by 6 mm, roughly equidistant from the eye (5 mm) and the snout (6 mm); blunt snout with a rounded profile; dorsolateral lines distinguishable, spaced by 11 mm; vomerine teeth present, slightly pigmented, V-shaped and well separated; four fingers (Nos 1–4), including protruding thumbs (No. 1) with dark nuptial pads, and respective size (inner to outer) as follows: 3 > 1 > 4 > 2; slender hind legs almost as long as the body, with femur 26 mm and tibia 33 mm, reaching the nostrils when projected on the head (‘leg test’); tarsus 30 mm with five elongated toes (up to 31 mm) with relatively extensive webbing, of formula I1-1II1/2-1/2III2/3-2IV3-3V; oval-shaped metatarsal tubercles (3 mm long by 2 mm high). Smooth skin; colourful dorsum with diffused dark blotches on a light-brown background; hints of creamy yellow on the flanks, groin, and armpits; predominantly white/light grey ventrum with faded dark spots and fine reddish speckles; throat greyish with dark reflections; iris yellow; temporal markings not pronounced (mostly brown, getting darker near the armpits); upper lip with reddish speckles on a light background, whiter near the tympanum and armpit. All observations and measurements were made post-mortem before fixation in ethanol.
Paratype:
MNCN 50656, adult female found by A.G. Nicieza and A. Peláez on 12 March 2021 in a small pool near Lagunas de Valdecarrín, Puerto de Las Señales, León, Spain (43.07429°, -5.23967°), and deposited at the Museo Nacional de Ciencias Naturales (MNCN) in Madrid, Spain.
Etymology:
The new nomen, Rana parvipalmata asturiensis, makes reference to the region of Asturias where this taxon is mainly distributed. As vernacular names, we recommend using ‘Asturian common frog’ for R. p. asturiensis, and to restrict the name ‘Galician common frog’ to the subspecies R. p. parvipalmata. To avoid confusion, the species R. parvipalmata as a whole may be referred to as the ‘Western common frog’.
Diversity and distribution:
Rana parvipalmata asturiensis is distributed in the western part of the Cantabrian mountains and coastal lowlands that mainly encompass the province of Asturias. It also mildly penetrates into the neighbouring provinces of Castilla y León (south) and Cantabria (southeast). In the west, the subspecies forms an extensive hybrid swarm with R. p. parvipalmata across Western Asturias, from about La Espina to the Eo River (at the border with Galicia). In the east, R. p. asturiensis meets and locally admixes with R. temporaria near the border between Asturias and Cantabria in the north, and near the border between Cantabria and Castilla y León in the south (Montaña Palentina). Populations of R. p. asturiensis are genetically structured due to isolation by distance (Choda 2014, Dufresnes et al. 2020b), but they feature relatively homogenous patterns of genetic diversity at neutral markers (Choda 2014, Cano et al. 2017, Dufresnes et al. 2020b). However, lowland and highland populations differ in genes related to metabolic function and immune response (Cano et al. 2017).
Natural history:
Not specifically studied, but presumably similar to other lineages of the R. temporaria complex. In Asturias, common frogs are found across a variety of environments including wetlands, rivers, puddles and ditches in forest tracks, small retention dams, and any kind of low depth water bodies, such as mountain and cattle ponds. The breeding season extends from early autumn to early summer, starting progressively later with higher elevation (Álvarez et al. 2012, Gutiérrez-Pesquera et al. 2022). The fine-scale spatial distribution may differ between highland areas, where breeding nuclei are larger but more isolated, and the more human-impacted mid-elevation and coastal areas, where breeding occurs in small nuclei in small and often anthropogenic water bodies.
Conservation:
The Asturian common frog is abundant and widespread throughout its range, hence it is probably not threatened despite its narrow distribution (c. 10 000 km2). The higher risks may be related to climate change and timber exploitation. Analyses of thermal tolerance and mechanistic niche modelling suggest a high vulnerability to climate warming, and suitability area projections point to an upward retreat of R. p. asturiensis populations, and extinction in most of the southern distribution ranges (Enriquez-Urzelai et al. 2018, 2019a, b, 2020, Gutiérrez-Pesquera et al. 2022). Although forestry promotes suitable water bodies associated with forest tracks, eucalypt tree plantations may have a negative effect because eucalypt leachates are detrimental to the growth and development of larvae (Burraco et al. 2018, Iglesias-Carrasco et al. 2022), and ultimately to metamorph survival (D. Álvarez and A.G. Nicieza unpublished data). Differences in functional genetic diversity along the altitudinal gradient (Cano et al. 2017) suggest that conservation actions should be population specific.
CONCLUSION
By improving our ability to detect genetic diversity, structure, and introgression, population genomics are a powerful approach to define, map, and rank phylogeographic lineages in taxonomy, notably by inferring reproductive barriers vs. connections and to see whether these relate to documented phenotypic differences. We showed that the Iberian endemic frog Rana parvipalmata is composed of two genetically and phenotypically well-defined phylogeographic lineages that seemingly admix freely across a shallow hybrid zone, thus implying conspecificity. Because formal taxonomic naming, diagnosis, and description is vital to inventory and protect biodiversity, these two lineages warrant distinction as subspecies, so they can be treated independently by conservation stakeholders, especially given their contrasting conservation situations. Our study also emphasized that morphometric differences and genetic differences are often disassociated, and that mitochondrial DNA barcoding may be misleading in the mapping of lineage distributions.
[Version of Record, published on 2 December 2023; http://zoobank.org/urn:lsid:zoobank.org:pub:EE53C938-C726-40CB-B117-5B9939453DA7]
ACKNOWLEDGEMENTS
We thank Andrés Peláez Cueto for help during fieldwork, as well as Sven Gippner for curating the cyt-b data. Sampling was authorized for the following authorities: Principado de Asturias (2006/000223, 2008/000272, 2010/000371, 2016/001092, 2017/001208, 2017/019842, 2018/001076, 2018/007781, PA/356/2020/144, AUTO/2019/3373, DECO/2021/2321, AUTO/2021/2141), Parque Nacional Picos de Europa (CO/09/0032/2005, CO/09/0007/2006, CO/09/646/2006, CO/09/077/2009, CO/09/0571/2009, CO/09/041/2011, CO/09/121/2012, CO/09/0125/2013, CO/09/012/2014, CO/09/065/2015, CO/09/0316/2015, PNP-1096/17-SCN, CO/09/073/2018, PNP-471/2018-SCN), Junta de Castilla y León (EP/CYL/389/2007, EP/LE/428/2010, EP/P/428/2010, EP/CYL/31/2010, EP/P/426/2010, EP/CYL/292/2013, EP/CYL/625/2013, EP/CYL/725/2015, EP/CYL/112/2017, EP/CYL/514/2018), Gobierno de Cantabria (EST-275/2016-SEP, EST-81/2017-SEP, EST-75/2018-SEP, CAP-001/2021), Xunta de Galicia (39/2008, 61/2011, 560/2011, 24/2014).
FUNDING
This study was funded by the Taxon-Omics priority program (SPP1991) of the Deutsche Forschungsgemeinschaft (grant VE247/19-1 to M. Vences and C. Dufresnes), and projects from Ministerio de Ciencia, Innovación y Universidades, Gobierno de España (MICINN grants CGL2012-40246 and CGL2017-86924-P to A.G. Nicieza).
DATA AVAILABILITY
Raw demultiplexed RAD-sequencing reads were uploaded on NCBI SRA under BioProject PRJNA949685 (accessions SAMN37663490-SAMN37663556); cyt-b sequences were deposited on GenBank (accessions OR640725-OR640892); the datasets used in the analyses were archived on Zenodo (doi: 10.5281/zenodo.8388713).
