Genomic Evidence for Cryptic Speciation in Tree Frogs From the Apennine Peninsula, With Description of Hyla perrini sp. nov

Dufresnes, Christophe; Mazepa, Glib; Rodrigues, Nicolas; Brelsford, Alan; Litvinchuk, Spartak N.; Sermier, Roberto; Lavanchy, Guillaume; Betto-Colliard, Caroline; Blaser, Olivier; Borzée, Amaël; Cavoto, Elisa; Fabre, Guillaume; Ghali, Karim; Grossen, Christine; Horn, Agnes; Leuenberger, Julien; Phillips, Barret C.; Saunders, Paul A.; Savary, Romain; Maddalena, Tiziano; Stöck, Matthias; Dubey, Sylvain; Canestrelli, Daniele; Jeffries, Daniel L.

doi:10.3389/fevo.2018.00144

ORIGINAL RESEARCH article

Front. Ecol. Evol., 02 October 2018

Sec. Phylogenetics, Phylogenomics, and Systematics

Volume 6 - 2018 | https://doi.org/10.3389/fevo.2018.00144

Genomic Evidence for Cryptic Speciation in Tree Frogs From the Apennine Peninsula, With Description of Hyla perrini sp. nov

$\r\nChristophe Dufresnes,,*$ Christophe Dufresnes^1,2,3^*

Glib Mazepa^1,4

Nicolas Rodrigues¹

Alan Brelsford^1,5

Spartak N. Litvinchuk⁶

Roberto Sermier¹

Guillaume Lavanchy¹

Caroline Betto-Colliard¹

Olivier Blaser¹

Amaël Borzée⁷

Elisa Cavoto¹

Guillaume Fabre¹

Karim Ghali¹

Christine Grossen¹

Agnes Horn¹

Julien Leuenberger¹

Barret C. Phillips¹

Paul A. Saunders¹

Romain Savary¹

Tiziano Maddalena⁸

Matthias Stöck⁹

Sylvain Dubey^1,3

Daniele Canestrelli¹⁰

Daniel L. Jeffries¹

¹Department of Ecology and Evolution, University of Lausanne, Lausanne, Switzerland
²Department of Animal and Plant Sciences, University of Sheffield, Sheffield, United Kingdom
³Hintermann & Weber SA, Montreux, Switzerland
⁴Department of Ecology and Genetics, Evolutionary Biology, Uppsala, Sweden
⁵Department of Biology, University of California, Riverside, Riverside, CA, United States
⁶Institute of Cytology, Russian Academy of Sciences, Saint Petersburg, Russia
⁷Division of EcoScience and Department of Life Sciences, Ewha Womans University, Seoul, South Korea
⁸Maddalena & Associati sagl, Gordevio, Switzerland
⁹Leibniz-Institute of Freshwater Ecology and Inland Fisheries (IGB), Berlin, Germany
¹⁰Department of Ecological and Biological Sciences, Tuscia University, Viterbo, Italy

Despite increasing appreciation of the speciation continuum, delimiting and describing new species is a major yet necessary challenge of modern phylogeography to help optimize conservation efforts. In amphibians, the lack of phenotypic differences between closely-related taxa, their complex, sometimes unresolved phylogenetic relationships, and their potential to hybridize all act to blur taxonomic boundaries. Here we implement a multi-disciplinary approach to evaluate the nature of two deeply-diverged mitochondrial lineages previously documented in Italian tree frogs (Hyla intermedia s. l.), distributed north and south of the Northern Apennine Mountains. Based on evidence from mitochondrial phylogenetics, nuclear phylogenomics, hybrid zone population genomics, niche modeling analyses, and biometric assessments, we propose that these lineages be considered distinct, cryptic species. Both mitochondrial and nuclear data affirm that they belong to two monophyletic clades of Pliocene divergence (~3.5 My), only admixing over a relatively narrow contact zone restricted to the southeast of the Po Plain (50–100 km). These characteristics are comparable to similarly-studied parapatric amphibians bearing a specific status. Inferred from their current geographic distribution, the two Italian tree frogs feature distinct ecological niches (<15% of niche overlap), raising questions regarding potential adaptive components contributing to their incipient speciation. However, we found no diagnostic morphological and bioacoustic differences between them. This system illustrates the speciation continuum of Western-Palearctic tree frogs and identifies additional cryptic lineages of similar divergence to be treated as separate species (H. cf. meridionalis). We recommend combined approaches using genomic data as applied here for the future taxonomic assessment of cryptic diversity in alloparapatric radiations of terrestrial vertebrates, especially in controversial taxa. Finally, we formally described the northern Italian tree frogs as a new species, Hyla perrini sp. nov.

urn:lsid:zoobank.org:pub:C446334C-3D94-4820-A851-99BD95084F6D

Introduction

Three decades of molecular biogeography have shed light on unsuspected amounts of genetic diversity among terrestrial vertebrate species that were once considered monotypic. In the Western Palearctic, this diversity was shaped by the combined actions of geological and climatic events (e.g., Mountain chain uplifts, the Messinian Salinity crisis, Quaternary glaciations), which have promoted allopatric diversifications between and within Mediterranean regions (e.g., Iberian, Italian, and Balkan Peninsulas, Northwestern Africa, Anatolia), and isolated off-shore islands (e.g., Crete, Cyprus) (Schmitt, 2007; Dufresnes, 2019). Yet, the resulting lineages often remained remarkably cryptic, sharing similar ecologies and featuring no obvious phenotypic differentiation, be it morphology, coloration or behavior, despite millions of years of independent evolution (Bickford et al., 2007, e.g., Wielstra et al., 2013; Diaz-Rodriguez et al., 2017).

What taxonomic status do such cryptic lineages deserve? This question is a topical controversy among zoologists and one of the major challenges in modern phylogeography (Bickford et al., 2007). Mayr's universal “biological species concept,” referring to pre- or post-zygotic barriers to reproductive isolation (Mayr, 1942), does not necessarily involved measurable morphological changes. Fortunately, secondary contact zones between pairs of taxa represent natural laboratories where cryptic species are put to the test (Hewitt, 2011). Pending opportunity for dispersal (Waters et al., 2013), a hybrid zone mediated by selection against hybrids will be narrower than expected under dispersal alone, which can be measured through analyses of genetic introgression. In many instances however, and notably when species have non-overlapping distribution ranges and do not naturally meet, the amount of genetic divergence (e.g., divergence time) is rather used as a proxy for their level of reproductive isolation. However, this indirect approach has several shortcomings. First, nuclear markers may fail to accurately resolve phylogenetic relationships, because of incomplete lineage sorting (ILS) and/or past hybridization events (e.g., during former interglacial expansions), especially for recent divergences (e.g., Pabijan et al., 2012). Second, the time and mode along which genetic divergence accumulates may greatly differ between vertebrate radiations and should therefore be independently assessed for the groups of interest. For instance, parapatric taxa may stop hybridizing after few million years in some systems (e.g., <3 My in Bufotes, Colliard et al., 2010), but long after in others (e.g., >10 Mya in Salamandrina, Hauswaldt et al., 2011; Mattoccia et al., 2011). Third, molecular clocks are extremely sensible to the calibration used in phylogenies, biasing divergence times by millions of years between equally plausible evolutionary scenarios (e.g., Veith et al., 2016). A standardized approach accounting for these limitations is still missing to harmonize species limits in cryptic polytypic radiations.

This issue is particularly relevant for Western Palearctic amphibians (Shaffer et al., 2015), which form so called “species complexes” that diversified throughout the late Miocene (11.0-5.3 Mya), Pliocene (5.3-2.6 Mya) and early Pleistocene (2.6-1.8 Mya), giving birth to evolutionary lineages of different ages, resulting in unresolved taxonomic situations. In recent years, genetic studies have proposed that some of these taxa be considered new, distinct species (e.g., Hyla felixarabica, Gvoždík et al., 2010; Pelodytes hespericus and P. atlanticus, Diaz-Rodriguez et al., 2017; Bufo spinosus, Arntzen et al., 2013; Lissotriton graecus, L. schmidtleri, L. lantzi and L. kosswigi, Pabijan et al., 2017; Triturus karelinii, T. anatolicus and T. ivanbureschi, Wielstra et al., 2013) or, on the contrary, disproved their potential specific status (e.g., Bufo (Bufotes) variabilis, Dufresnes et al., 2018). Many lineages have yet to receive a similar assessment, upon which their taxonomic fate will depend.

With at least eight taxa distributed around the Mediterranean Basin, the Western Palearctic radiation of hylid tree frogs (Hyla arborea species complex) is an ideal system for the study of cryptic speciation (Stöck et al., 2008b, 2012). Most of these lineages are remarkably similar in appearance, mating calls and ecology, despite their Mio-Pliocene origins and signs of reproductive isolation as revealed by narrow transitions at their parapatric margins (see Discussion). Cryptic taxa were further suspected within the north-African populations of H. meridionalis (Recuero et al., 2007; Stöck et al., 2008b, 2012), as well as across the Apennine Peninsula in the Italian endemic H. intermedia (sensu lato, Canestrelli et al., 2007a,b), pending in-depth characterization.

Here we focus on the latter, H. intermedia s. l., which ranges from Sicily to the southern slopes of the Alps in northern Italy and the Swiss canton of Ticino. This taxon features a typical phylogeographic pattern for Italian amphibians, with distinct mitochondrial clades distributed north and south of the Northern Apennine Mountains (clades N and C/S in Canestrelli et al., 2007b), tentatively split as H. intermedia sensu stricto and Hyla “new taxon 2” by Stöck et al. (2008b, 2012). Both mitotypes are found in sympatry along the Apennine foothills from central Italy (Emilia-Romagna), over a transition that appears no wider than tens of kilometers (Canestrelli et al., 2007b). Their nuclear identity, however, has remained cumbersome. Population genetics using allozyme markers accordingly suggested latitudinal population differentiation, yet with massive allele sharing, indicative of ancestral polymorphism or widespread admixture (Canestrelli et al., 2007a). Similarly, subsequent phylogenies based on slowly-evolving nuclear sequences showed little support for monophyletic clades (Stöck et al., 2008b, 2012), supposedly because of past or present gene flow (Gvoždík et al., 2015) or purifying selection on gene markers (e.g., rag-1, Stöck et al., 2008b). Therefore, their evolutionary and taxonomic situation is still to be determined.

In the present study, we provide several lines of evidence that Italian tree frogs belong to two separate species, and we formally describe the northern one as Hyla perrini sp. nov. (see specific section in the Discussion). We combined population genomics and nuclear phylogenomic inferences to show that both species are monophyletic, deeply-diverged, and only admix across a well-defined hybrid zone restricted to their post-glacial contact in Emilia-Romagna, north-central Italy. We further implement analyses of ecological niche modeling, morphometrics, and bioacoustics to explore phenotypic variation between these two cryptic species.

Methods

DNA Collection

A total of 171 tissue samples of tree frogs were obtained from non-invasive buccal swabs (adults) or tail tips (tadpoles). Some were already available from our previous genetic studies (Canestrelli et al., 2007a,b; Stöck et al., 2008b, 2011, 2012; Dufresnes et al., 2011, 2013, 2015b,c, 2016b,c) while others were collected during fieldwork campaigns in Northern Italy and Southern Switzerland in the springs of 2013 and 2018. DNA was extracted using the Biosprint robotic workstation (Qiagen) or the DNeasy Blood and Tissue Kit (Qiagen).

RAD-Sequencing Genomic Libraries

Two separate genomic libraries were prepared using the double digest RAD (ddRAD) protocol detailed in Brelsford et al. (2015). Briefly, genomic DNA (5–15 ng/μL, 30–90 ng in total) was digested using Sbf1 and Mse1 restriction enzymes. Adapters containing PCR primer sites and individual barcodes (4–8 nucleotides) were then ligated to the Sbf1 end of the fragment. Barcoded fragments were purified and amplified using PCR with Illumina primers. The amplified libraries were then size selected to enrich for fragments between 400 and 500 bp using gel electrophoresis (2.5% agarose gel run for 2 h 30 at 80 V) and finally, the libraries were purified once more and sequenced (single-end only) on an Illumina Hi-Seq 2500.

Library 1 was designed to reconstruct the nuclear phylogeny of this radiation. It included 51 samples representative of all Western-Palearctic taxa (plus Hyla (Dryophytes) japonica as an outgroup; Table S1B). Particular focus was given to H. perrini sp. nov. (n = 12, including the holotype and two paratypes) and H. intermedia s. s. (n = 7); 3–6 individuals for every other species (Table S1B). Raw data from one sequencing lane was processed as follows. Raw, single end Illumina sequences were quality checked using FastQC v0.10.1 (Andrews, 2010) and the individuals of the library were demultiplexed by individual sample barcodes with the process_radtags.pl module of STACKS v 1.42 (Catchen et al., 2013). The resulting individual fastq files were mapped on the draft genome assembly of H. arborea from Brelsford et al. (2015) with bowtie2 v 2.3.1 (Langmead and Salzberg, 2012; using –local –highly-sensitive specifications. Resulting sam files were subsequently filtered for mapping quality >20 with samtools v 0.1.19 by Li (2011) and the bam files were used as input for ref_map.pl module of STACKS with default parameters. The STACKS catalog obtained consisted of 181,040 loci merged across all individuals (mean of 67,000 loci/individual with average coverage of 31 reads/locus), which was then whitelisted down to 621 loci (biallelic per individual) present across at least 50 individuals. The populations module of STACKS was then run with –phylip_var_all to obtain individual-based alignment (using a population map comprising all 51 individuals as separate “populations”) of 42,834 bp. This alignment was then translated into nexus format with Geneious v8.1.9 (https://www.geneious.com/). Finally, in order to compute genome-wide diversity and divergence between H. perrini sp. nov. and H. intermedia s. s., we outputted a second alignment containing all tags sequenced in both species (2.9 Mb from 29,294 loci).

Library 2 was designed for population genomic analyses of the hybrid zone between the two Italian lineages and included 120 samples from 12 populations sampled across the Apennine Peninsula (Table S1C). Raw, single end illumina sequences from two lanes were again quality checked using FastQC v0.10.1 and demultiplexed using the process_radtags module of STACKS. RAD loci were constructed and SNPs were called using STACKS. We first ran initial tests of parameter values for the core STACKS parameters [Ustacks: -M (2-4), -m (3-6) and cstacks: -n (2-4)] to test for any major effects on the number of loci, average coverage, and number of heterozygous genotype calls (a proxy for genotype call accuracy). The default values (-M 2, -m 3 and -n 2) provided a good balance between data quality and quantity according to the criteria set out in Paris et al. (2017). Modules ustacks, cstacks, sstacks, and populations were then run separately using these parameter values. The STACKS catalog contained 670,396 loci; these were filtered using the populations module and loci were retained if they had a minor allele frequency filter (–min_maf) of >0.05, a maximum observed heterozygosity filter (–max_obs_het) < 0.75 (to accounts for cases of “over merging” of paralogous loci), if they were present in at least 80% of individuals from all 12 populations and had a coverage of at least 6 reads per locus. The resulting dataset contained 2,431 SNPs from 1,639 RAD tags, with an average coverage of 30.45 reads (SD = 11.77).

Mitochondrial Data

For phylogenetic purposes, we harvested published sequences for the mitochondrial genes cytochrome b (cyt-b, 840 bp) and cytochrome c oxidase subunit I (COI, 766 bp) in 33 individuals representative of the Western Palearctic diversity of tree frogs (2–6 individuals per taxa), as well as 3 outgroups (H. japonica, H. cinerea, H. squirella). Details on their origin can be found in Table S1. Importantly, we only selected individuals with sequences available for both genes that could thus be concatenated for analysis (total alignment length = 1,606 bp).

Furthermore, distribution data of the two Italian mitochondrial lineages (H. perrini sp. nov. and H. intermedia s. s.) were gathered from previous studies (Canestrelli et al., 2007a,b; Verardi et al., 2009; Stöck et al., 2012; Dufresnes et al., 2015c), totaling 45 populations.

Phylogenetic Analyses

For each mitochondrial (1,606 bp from two partitioned genes) and genomic (42,834 bp from 621 RAD tags concatenated in a supermatrix) alignments separately, we independently performed time-calibrated Bayesian phylogenetic reconstructions in BEAST 2.4.8 (module StarBEAST; Bouckaert et al., 2014), with a relaxed lognormal molecular clock and a birth-death tree prior. Site substitution models were GTR+G (nuclear), TN93+G+I (COI) and HKY+G (cyt-b), selected with bModelTest (Bouckaert and Drummond, 2017). We used three calibration points with normally distributed priors. The first two were inferred from multiple fossil calibrations, independently applied by Smith et al. (2005) and Li et al. (2015): the split between the H. japonica group and Western Palearctic Hyla (36 Mya ± 6 My), and the time since the most recent common ancestor of Western Palearctic tree frogs (20 Mya ± 3 My). In order to stabilize the clock, we used a third, younger calibration: the Messinian divergence of the Tyrrhenian endemic H. sarda (set to 5.33 Mya ± 2 My), which is the most plausible origin for this insular species, based on phylogeographic (Stöck et al., 2008b, 2012; Bisconti et al., 2011) and fossil evidence (Kotsakis, 1980). The chains were run for 50 million (nuclear) and 250 million iterations (mtDNA), sampling every 5,000 (nuclear) and 25,000 iterations (mtDNA), and checked for stationarity and effective sample sizes (ESS) of parameters >200 with Tracer 1.5. The first 10% of trees were discarded as burnin, and the remaining 9,000 were analyzed by TreeAnnotator 2.4.8, and visualized by DensiTree 2.2.6 (Bouckaert and Heled, 2014).

Furthermore, we explored the relationships between Western-Palearctic hylids with a Principal Component Analysis (PCA) on the nuclear genomic dataset (haplotype VCF file analyzed in adegenet R package), discarding the early-diverged H. japonica and H. meridionalis. Finally, we estimated nucleotide diversity (Π) for Hyla perrini sp. nov. and H. intermedia, as well as the percentage of fixed differences between them.

Population Structure Analyses

In order to characterize the nuclear population structure between Hyla perrini sp. nov. and H. intermedia s. s. and within the contact zone, we performed STRUCTURE (Pritchard et al., 2000) analyses using the SNP dataset from library 2. STRUCTURE was run from a k of 1–12, with an MCMC chain length of 10,000,000 (burn-in 1,000,000) and three iterations were performed for each value of k. STRUCTURE outputs were fed to STRUCTURE HARVESTER (Earl and vonHoldt, 2012) to identify the most appropriate value of k using the Evanno method (Evanno et al., 2005). We further computed pairwise genetic distance (F_st) between pure populations (loc. 1 and 2 for H. perrini sp. nov.; loc. 12 for H. intermedia s. s.), as well as their respective diversity (observed and expected heterozygosity H_o and H_e) with VCFtools (http://vcftools.sourceforge.net/).

Additionally, we analyzed the modality of the Hyla perrini sp. nov./H. intermedia s. s. transition trough cline analyses, which are powerful tools to understand the dynamic of hybrid zones in the geographic space, through the estimation of cline width (w) and center (c). To this purpose, we fitted sigmoid clines to mitochondrial frequency data (n = 10 populations, totaling 70 individuals) and nuclear hybrid index (STRUCTURE admixture coefficient; n = 12 populations, totaling 120 individuals) along a transect running from Southern Switzerland (pure H. perrini sp. nov.) to Abruzzo, Italy (pure H. intermedia s. s.), using the R package hzar (Derryberry et al., 2014). We performed model selection between two-parameter clines (w and c) and models with exponential tails (up to eight parameters) and kept the model with the best AIC score.

Niche Modeling Analyses

To examine potential ecological differences between both Hyla species, as well as to obtain an approximate age for the contact zone, we modeled their present and past (Last Glacial Maximum) distributions with MaxEnt (v. 3.3.3k; Phillips et al., 2006). We combined 84 localities of H. intermedia and 154 localities of Hyla perrini sp. nov. (Figure S1), comprising our own and previously published records, as well as localities from the Global Biodiversity Information Facility database information (http://www.gbif.org/). Duplicated localities were removed by ENMTools v. 1.3 (Warren et al., 2010).

Nineteen bioclimatic layers were extracted from the WorldClim 1.4 database (http://www.worldclim.org). An additional nine uncorrelated variables were included: Altitude (WorldClim 1.4), Global Aridity and Potential Evapo-Transpiration (aridity index; http://www.cgiar-csi.org/data/global-aridity-and-pet-database), GlobCover 2009 (Global Land Cover Map; due.esrin.esa.int/globcover/), spatial heterogeneity of global habitat (EarthEnv, http://www.earthenv.org/texture.html), percent tree cover (Github, https://github.com/globalmaps/gm_ve_v1), as well as four topographic layers (aspect, exposition, slope, and terrain roughness index) calculated in QGIS v. 2.18.1 (http://www.qgis.org/).

To eliminate predictor collinearity prior to generating the models, we calculated Pearsons's correlation coefficients for all pairs of bioclimatic variables using the ENMTools. We excluded the variable of a correlated pair with |r| > 0.8 that we considered to be the less biologically important of the pair, based on known preferences of Hyla species. The resulting dataset hence contained 10 uncorrelated bioclimatic variables, plus the nine additional variables mentionned above. We used a jackknife analysis for estimating their relative contributions to the MaxEnt model (Table S2) and applied a mask that extends from 35° to 48° N and 5° to 20° E.

These 19 variables were implemented in the final models for both species (30 arc-seconds resolution). We used 70% of the occurrence localities as training data, and the remaining 30% were kept for testing the resulting models. Model performance was measured using the Area Under the Curve (AUC) derived from the Receiver Operating Characteristic (ROC) plots. The plots represent the model's ability to discriminate species locations from pseudo-absences by plotting sensitivity against 1 - specificity. AUC values range from 0.5 to 1.0, with 0.5 indicating no greater fit than expected by chance and 1.0 indicating a perfect model fit. Models with test AUC values above 0.75 are considered useful and above 0.90 very good (Swets, 1988; Elith, 2000). To properly parameterize the model, we evaluated the performance of various combinations of ten regularization multipliers (from 0.5 to 5, in increments of 0.5). The best-fit models were parameterized with a regularization multiplier of 0.5 for H. intermedia and 1.5 for H. perrini sp. nov. We used the default auto feature setting within the MaxEnt software to generate the model. The niche overlap between species was then estimated using Schoener's D distance in ENMTools (Schoener, 1968), which is the best suited to compute niche overlaps from potential distributions derived from species distribution models (Rödder and Engler, 2011).

The ten selected uncorrelated bioclimatic layers (2.5 min resolution, extracted from WorldClim 1.4) were further used to predict the Last Glacial Maximum (LGM) distributions. Two general atmospheric circulation models were used to generate LGM climate scenarios: the Community Climate System Model (CCSM; http://www.cesm.ucar.edu/models/ccsm4.0/ccsm/) and the Model for Interdisciplinary Research on Climate (MIROC; Watanabe et al., 2011). Models were generated by MaxEnt using the above described settings.

Biometric Analyses

Morphometrics

We tested whether the newly described taxon H. perrini sp. nov. exhibits significant differences in morphology and breeding calls compared to its sister species H. intermedia s. s. and its other parapatric relative H. arborea. In spring 2018, male tree frogs from two Hyla perrini sp. nov. (n = 24, including the type locality), four H. intermedia s. s. (n = 15) and two H. arborea populations (n = 22), were captured, photographed and measured for nine variables commonly used in amphibian taxonomy: snout-vent length (SVL), hind leg length (HL), tarsus + tibia length (TTL), head width (HW), eye diameter (ED), tympanum diameter (TD), inter-eye distance (EE), inter-nostril distance (NN) and eye-nostril distance (EN). We used a 1 mm accurate ruler with a blocking end for the first three variables, and a digital caliper (0.1 mm accuracy) for the remaining six variables. For the sake of consistency, all individuals were measured by the same person (CD), except for six of the H. intermedia s. s. individuals (DC). Measures were corrected by the relative size of individuals (i.e., divided by SVL) and compared between the three species by means of Principal Component Analyses (PCA) and Multivariate Analysis of Variance (MANOVA), which allowed for the inclusion of between-population intraspecific variation as a cofactor. Variables showing significant difference between species were further compared by non-parametric Tukey tests.

Bioacoustics

Mating calls were recorded directly at the breeding sites mentionned above using a portable digital recorder (EVISTR L57, setup to 48,000 Hz). A total of 11 different males of Hyla perrini sp. nov., 10 of H. intermedia and 10 of H. arborea could be recorded. We pre-edited the calling sequences of each individual with the software Audacity (http://audacity.fr/), and used the R package Seewave (Sueur et al., 2008) for downstream processing and data analyses.

The calls of most European hylids, including Italian taxa, consist of series of short monotonous notes, each made of partially-fused pulses (Rosso et al., 2006). For each male, we visualized and analyzed a set of ten consecutive notes in the middle of a call series, using the built-in functions of Seewave (spectro, oscillo, a.timer, and fpeaks). We averaged note duration (ND), duration of pause between two consecutive notes (PD), pulse rate (PUL) and measured the frequency peaks of the two main energy bursts (fundamental frequency F_f at ~1 kHz and dominant frequency F_d, at ~2 kHz). The actual band-width of the call is wider but cannot be reliably analyzed because of background noise and/or because the signal was too weak in some recordings. Temperature, known to affect call properties (Köhler et al., 2017), was recorded. Other important abiotic variables such as time of the day (early evening: 21−23 h), moment in the breeding season (main peak of activity: mid-April 2018 and early May 2018 for south-alpine and north-alpine populations, respectively), and motivation context (chorus of >15 males) were all similar and should have limited effect on our comparisons. In the same ways as for morphometric analyses, we explored and tested species differences with PCA and MANOVA, considering populations and temperature as co-variables for the latter.

Collection of Type Specimens

Two males and two females were collected in April and May 2018 from the type locality of H. perrini sp. nov, in the Swiss canton of Ticino (see details in Discussion). A third male, found dead but in a good state of preservation, was also included in the type series. All five were measured and DNA was collected with buccal swabs; the holotype and two paratypes were included in the phylogenomic dataset (Table S1B). Specimens were subsequently anesthetized and euthanized by overdose of MS-222 0.015% (Sandoz) buffered in 0.03% sodium bicarbonate, and inventoried in the herpetological collection of the Cantonal Zoological Museum of Lausanne (MZL42340-42344). The type series was later examined post-mortem, including photography and 14 standard morphometric measurements with a digital caliper (0.1 mm accuracy), made by a single one of us (GM).

Data Accessibility

The raw sequence reads from both RAD libraries are available on the NCBI SRA archive under BioProject PRJNA485248. Phylogenomic, population genomic, ecological niche modeling and biometric data files used in the analyses are provided in Data Sheet S1.

Results

Phylogeny of Western-Palearctic Tree Frogs

Time-calibrated phylogenies of mitochondrial (1,606 bp from two genes) and nuclear sequences (42,834 bp from 621 RAD tags) congruently recovered all known Hyla taxa from the Western Palearctic (Figures 1, 2). The nuclear tree is fully resolved (Figure 2), suggesting no signs of past admixture or incomplete lineage sorting, while the mtDNA tree differs slightly at a few unresolved nodes (notably H. sarda) and larger uncertainty in the molecular dating (Figure 1). Both trees fully support the split of H. intermedia s. s. and Hyla perrini sp. nov. as two monophyletic clades, with divergence pointing to the upper Pliocene (nuDNA: 3.5 Mya, 95% HPD: 5.0-2.2; mtDNA: 2.5, 95% HPD: 4.2-1.0). The first two axes of the PCA (based on individual genotypes) recovered most known species, including H. intermedia and H. perrini sp. nov. which form adjacent yet distinct clusters (Figure S2). These two Italian taxa featured 9.8% of fixed difference at COI and 8.7% at cyt-b (9.2% for the concatenated alignment). At the nuclear level, they differ by 0.3%, based on 2.9 Mb. Nucleotide diversity (Π) was estimated at 0.0061 for H. intermedia s. s. and 0.0069 for H. perrini sp. nov.

FIGURE 1

Figure 1. Time-calibrated phylogeny of mitochondrial sequences from cyt-b (840 bp) and COI (766 bp). Clock rates were nearly identical between markers (2.5e⁻² substitutions/site/My). Diffuse lines show branches of each of the sampled trees, summarized by the root canal (blue framework). For major nodes, the Bayesian posterior probabilities and confidence intervals of divergence times are indicated. Species distributions use the color codes shown on the tree. The map was built in QGIS 2.18 using map layers from ESRI and the IUCN red list (adapted by Dufresnes, 2019). Hyla japonica, H. cinerea and H. squirella were used as outgroups (Table S1A).

FIGURE 2

Figure 2. Time-calibrated phylogeny of nuclear loci (42,834 bp from 621 RAD tag sequences). Clock rate averaged 3.8e⁻⁴ substitutions/site/My. Diffuse lines show branches of the sampled trees, summarized by the root canal (blue framework). For major nodes, the Bayesian posterior probabilities and confidence intervals of divergence times are indicated. Species distributions use the color codes shown on the tree. Species distributions use the color codes shown on the tree. The map was built in QGIS 2.18 using map layers from ESRI and the IUCN red list (adapted by Dufresnes, 2019). The polytypic H. japonica was used as outgroup (Table S1B).

Hybrid Zone Between Apennine Tree Frogs

MtDNA data in across the distribution clearly delineates the ranges of Hyla perrini sp. nov. (Northern Italy and Southern Switzerland) and H. intermedia s. s. (Central-Southern Italy and Sicily) (Figure 3), which are mostly kept apart by the Apennines, but locally meet in southeastern Emilia-Romagna. Population genomics (RAD data) yielded a similar nuclear signal (based on 2,431 SNPs from 1,639 RAD tags, Figure 4): STRUCTURE analyses with k = 2 (the most likely number of k) perfectly recover the two gene pools, with a spectrum of intermediate coefficients for frogs sampled at the transition, indicative of genetic admixture (loc. 4–11 in Figure 4).

FIGURE 3

Figure 3. Distribution of the mitochondrial lineages of H. perrini sp. nov. (green) and H. intermedia s. s. (blue) across their natural ranges, based on records from previous studies. Dashed lines delineate range margins. The map was built in QGIS 2.18 using map layers from ESRI.

FIGURE 4

Figure 4. Genetic structure across the H. perrini sp. nov. / H. intermedia s. s. transect, inferred from Bayesian clustering of individuals (n = 120 from 12 populations) in k = 2 groups (STRUCTURE), based on 2,431 loci. Horizontal barplots show individuals assignments to the H. perrini sp. nov. (green) or H. intermedia s. s. (blue) gene pools; pie charts feature the population averages, with size proportional to sample size. The map was built in QGIS 2.18 using map layers from ESRI.

Accordingly, steep clines of less than 100 km (Figure 5) were inferred for the nuclear hybrid index (w = 96 km; 95% CI: 64–153) and the frequencies of mitotypes (w = 41 km; 95% CI 20–92). Both clines remarkably coincided, with centers at the northeastern foothills of the Apennines (between population 7 and 8, Figure 5), shifted by only 2 km (for nuclear: c = 378 km along north-south transect; 95% CI: 365–394; for mtDNA: c = 380 km: 95% CI: 365–392). No trace of genetic introgression was found apart from the contact zone.

FIGURE 5

Figure 5. Geographic clines fitted to nuclear hybrid index data (STRUCTURE admixture coefficient, left panel) and mitochondrial frequency data (right panel) along the H. perrini sp. nov./H. intermedia s. s. transect. The best selected models featured only two parameters (width and center), while parental frequencies could be set to 0 and 1. Nuclear cline: width w = 96 km and center c = 378 km; mitochondrial cline: w = 40 km and c = 380 km.

Pure populations of H. perrini sp. nov. featured nearly identical observed (H_o) and expected (H_e) heterozygosity (loc. 1, the type locality: H_o and H_e = 0.18; loc. 2: H_o and H_e = 0.17). Similar values were obtained for H. intermedia s. s. (H_o = 0.16, H_e = 0.15). Pairwise genetic distance between the two species was high (loc. 1–12, F_st = 0.49; loc. 2–12, F_st = 0.46), compared to between the two H. perrini sp. nov. populations (loc. 1-2, F_st = 0.09).

Ecological Niche Modeling

The present time MaxEnt models for H. intermedia and Hyla perrini sp. nov. received robust evaluation metrics (AUCtest = 0.899 ± 0.023 and 0.939 ± 0.017 respectively) and showed significance for the binomial omission test, indicating a good performance of the model. The predicted niches under the current climate conditions are shown in Figure 6, and the contribution of variables shown in Table S2. These remarkably differ between the two species (D = 0.092), and only overlap on the northern slopes of the Apennine Mountains.

FIGURE 6

Figure 6. Projection of the ecological niches for H. perrini sp. nov. (Left) and H. intermedia s. s. (Right) under current (top) and Last Glacial Maximum (LGM; 21'000 years ago) environmental conditions from the Community Climate System Model (CCSM, middle) and the Model for Interdisciplinary Research on Climate (MIROC, bottom). Parameter contributions are available in Table S2. The maps were built in QGIS 2.18 using map layers from ESRI.

The CCSM and MIROC models both suggested that the LGM environmental conditions were most suitable along the Tyrrhenian and southern Adriatic coasts for H. intermedia s. s., while the inlands seem unfit (Figure 6). In contrast, within current range limits, the predicted distribution of glacial Hyla perrini sp. nov. populations is essentially restricted along the Alpine edge of Po Plain (Figure 6). Additional potential areas highlighted by the models (especially the CCSM), like the Tyrrhenian islands and the Balkan coast, are irrelevant given that these are inhabited by other Hyla species (H. sarda and H. arborea). The current area of contact between the two species does not seem to have hosted suitable LGM conditions.

Biometrics

Morphometric Variation

Analyses of nine variables (corrected by SVL) showed little morphological differentiation between Hyla perrini sp. nov., H. intermedia s. s., and H. arborea (Figure 7A, Figure S3). Multivariate analyses of variance (MANOVA), accounting for intra-population variation, nevertheless highlighted a significant species effect for five variables after Bonferroni correction: HW, TD, EN, NN, and HL (Figure S3). Among them, Hyla perrini sp. nov. differs from H. intermedia s. s. by a smaller head (HW) and shorter hind legs (HL) (Tukey tests). Both Italian taxa differ from H. arborea by a smaller tympanum (TD) and a shorter eye-nostril distance (EN) (Figure S3) (Tukey tests).

FIGURE 7

Figure 7. Principal Component Analyses (PCA) on biometric variables of H. perrini sp. nov. (green circles), H. intermedia s. s. (blue triangles) and their close relative H. arborea (black crosses). (A) PCA on morphometric data. The photos illustrate the nine measurements included (see section Methods). (B) PCA on bioacoustic data. The lower graphs illustrate three notes of a call from H. perrini sp. nov., showing the spectral of temporal features included. Note that the temperature was lower (12°C) for the six H. intermedia s. s. with a longer pause duration and lower pulse rate (negative values on axis 2), while all other individuals were recorded at >15°C.

Bioacoustic Analyses

Call structure and parameters were globally similar between species (Figure 7B). Among the five bioacoustic variables measured, two showed a significant species effect (F_f and PD), while accounting for populations and temperature as co-variables (which were both significant). Species-by-species comparisons failed to find a significant pair for F_f; pause duration (PD) was longer in H. intermedia s. s. than Hyla perrini sp. nov. (Figure S4). This is illustrated by the second axis of the PCA (Figure 7B), although the six males driving this axis were all recorded on a colder night (11°C) compared to the rest of the samples (>15°C).

Discussion

Cryptic Speciation in Apennine Tree Frogs

Two lines of evidence indicate that the Italian tree frog Hyla intermedia s. l. consists of distinct evolutionary entities, distributed across the southern and the northern parts of the Northern Apennine Peninsula, hereby named Hyla intermedia s. s. and Hyla perrini sp. nov., respectively (see its description below). First, the mitochondrial clades previously identified (Canestrelli et al., 2007a,b; Stöck et al., 2008b, 2012) are fully supported by our nuclear loci (Figure 2). Both lineages are monophyletic and diverged at Pliocene times (Figures 1, 2). Second, population genomic analyses unambiguously distinguished the two gene pools across the Italian Peninsula, and only detected admixture over less than 100 km along a well-defined corridor between the northeastern Apennines and the Adriatic coast (Figures 3–5). These results imply that the weak nuclear identity reported by previous studies (Canestrelli et al., 2007b; Stöck et al., 2008b, 2012; Gvoždík et al., 2015) stemmed from low marker resolution and persistent ancestral polymorphism, rather than frequent gene flow and widespread introgression.

Italian Hyla most likely differentiated due to the continued uplift of the Apennines since the lower Miocene, some 20 Mya (Picotti and Pazzaglia, 2008). Moreover, what is now the Po Plain was flooded as part of the Apennine foredeep after the Messinian (Pellen et al., 2017), which may have contributed to the initial stages of divergence. Thermophilic species like tree frogs survived the ensuing Pleistocene glaciations in isolated refugia. Phylogeographic analyses identified a single one for Hyla perrini sp. nov. (northern Italy) but at least two for H. intermedia s. s. (Sicily/Calabria and Central Italy; Canestrelli et al., 2007b). According to our ecological niche modeling, these were probably nested in subalpine valleys (Hyla perrini sp. nov.) and along the paleo-coasts (H. intermedia s. s.; Figure 5). Their secondary contact in north-central Italy is thus mostly likely post-glacial, following demographic expansions after the LGM (Canestrelli et al., 2007b). Tree frogs are a text-book example of the phylogeography of peninsular Italy, where the Apennines triggered various north/south divergences in many widespread amphibians, eventually resulting in different species/subspecies (e.g., Pelophylax lessonae / bergeri, Canestrelli and Nascetti, 2008; Bombina variegata / pachypus; Canestrelli et al., 2006b; Salamandrina perspicillata/terdigitata, Canestrelli et al., 2006a; Triturus carnifex, Canestrelli D Salvi et al., 2012; Lissotriton vulgaris meridionalis, Maura et al., 2014; Rana dalmatina, Canestrelli et al., 2013).

Our complementary genetic approaches suggest that these two taxa are worthy of a specific status, in light of the current taxonomic divisions of other European amphibians. First, their deep, congruent mtDNA and nuclear divergence (averages: 2.5 and 3.5 Mya) is comparable to or older than that of several already established species, such as the Iberian and Oriental tree frogs, H. molleri and H. orientalis (~3 Mya; Stöck et al., 2012, this study); the Parsley frogs Pelodytes hespericus, P. atlanticus and P. punctatus (2–3 Mya; Diaz-Rodriguez et al., 2015, 2017); the green toads Bufo (Bufotes) viridis, B. siculus and B. balearicus (2–3 Mya; Stöck et al., 2008a; Dufresnes et al., 2014); the pool frogs Pelophylax lessonae and P. bergeri (1–2 Mya, Canestrelli and Nascetti, 2008), the Cyprus water frog P. cypriensis (2–3 My, Lymberakis et al., 2007, but 5.3 Mya in Plötner et al., 2012), the Caucasian toad Bufo verrucossisimus (1–2 Mya; Recuero et al., 2012); it is barely younger than many Messinian taxa, e.g., H. arborea (5–6 Mya, Stöck et al., 2012, this study), Middle Eastern Triturus newts (Wielstra and Arntzen, 2011). Note that molecular dating is inherently sensitive to the marker(s) and calibration(s) used (e.g., Veith et al., 2016) and should obviously be interpreted cautiously. Here the mitochondrial phylogeny is not fully resolved, particularly for H. sarda (one of our calibration points); hence the nuclear estimate is probably more trustworthy than the mitochondrial one, which also stands out from their respective confidence intervals.

Second, the restricted introgression along a well-defined contact zone and clear-cut clines argue for limited gene flow between the two gene pools. While the Apennine Mountains kept the two lineages apart over most of the range, the area of contact features no obvious barriers to dispersal. Under a model of neutral diffusion (Barton and Gale, 1993; see Dufresnes et al., 2015a), and assuming a dispersal distance σ of 1.5 km per year (Vos et al., 2000), the cline width w measured from our nuclear (w = 96 km) and mitochondrial data (w = 41 km) would be exceeded in about 650 and 120 generations since first contact, respectively. Accounting for a generation time of 1–2 years, these timing estimates are shorter than the presumed post-glacial contact, suggesting partial reproductive barriers, and that the hybrid zone has reached equilibrium. Selection coefficients s for an averaged nuclear gene, estimated from w and σ (w ∞ σ/√s; Barton and Gale, 1993) were relatively weak, i.e., s = 0.002 under intrinsic selection (heterozygote disadvantage) and s = 0.001 under extrinsic selection (ecological disadvantage). These modest estimates reflect the fact that most genes probably face little or no selection against hybrids, and so that reproductive isolation may involve only a few genes (that have lost linkage with neutral ones) since the post-glacial contact. Male-biased dispersal, expected in polygynous, lek-breeding animals such as tree frogs (Lampert et al., 2003; and documented in European Hyla, Vos et al., 2000), could account for the wider nuclear transition compared to the mitochondrial one. Because both clines coincide, asymmetric fitness between crosses appears unlikely. Analyses of sex-linked vs. autosomal markers may shed some light on the nature of potential hybrid incompatibilities with respect to sex-biased dispersal (Sciuchetti et al., 2018). In parallel, estimation of selection at the locus level should inform on the genomic architecture of speciation and identify the genes/genomic regions involved.

This neat transition falls within the range of those reported between other interspecific lineages that still hybridize at their parapatric margins, i.e., Bufo (Bufotes) viridis and B. balearicus (50–100 km, Dufresnes et al., 2014), Bufo bufo and B. spinosus (50–100 km, Arntzen et al., 2018), Triturus ivanbureschi and T. anatolicus (<50 km, but with hundreds of kilometers of introgression due to hybrid zone movement, Wielstra et al., 2017). In European hylids, it is similar to that of the Messinian diverged taxa, H. arborea and H. orientalis (40 km in the Balkans but >100 km in Poland, Dufresnes et al., 2015a, 2016a), H. arborea and H. molleri (>100 km, Stöck et al., 2012, CD and MS unpublished data), suggesting that the speciation continuum in tree frogs spans from the Pliocene to the upper Miocene (Dufresnes et al., 2015a). Given such a timeframe, the Algerian/Tunisian lineage H. cf. meridionalis (“new taxon 1” in Stöck et al., 2008b, 2012) also deserves taxonomic attention, according to our phylogenomic analysis (Figure 2).

The clear-cut Hyla perrini sp. nov./intermedia s. s. transition might be bounded by differential ecological adaptation. The ecological niches of the two taxa are highly disruptive, with projected distributions overlapping only at range margins (Figure 6). A similar pattern was evidenced for other parapatric amphibians meeting in Central Italy (e.g., Lissotriton vulgaris and L. italicus, Iannella et al., 2017). Yet, these differences might also stem from their narrow realized niches, bounded to different geography, i.e., H. intermedia s. s. inhabiting the southern, dry Mediterranean Italy while Hyla perrini sp. nov. is restricted to the more temperate Po Plain. European Hyla species can be quite generalist, e.g., the genetically homogeneous Hyla arborea ranges from southern Greece to Northern Germany (Stöck et al., 2012; Dufresnes et al., 2013). Populations from Hyla perrini sp. nov. also persisted for decades north of the Alps following illegal introductions (Dufresnes et al., 2015c, as H. intermedia).

Are these two species really cryptic? Unsurprisingly, we quantified little morphometric and bioacoustic differences (Figure 7). Amphibians are known for remarkable resemblances between sibling taxa, hence the need for molecular tools to identify them (Dubey et al., 2018). Morphological variation is rather associated to ecological conditions such as climate and geography (as in H. intermedia s. l., Rosso et al., 2004, and other Hyla, Gvoždík et al., 2008), stemming from plasticity or convergent adaptation (e.g., the “spinosus” morphotype found in southeastern European Bufo bufo, Garcia-Porta et al., 2012; Recuero et al., 2012), and therefore not reflecting genetic divergence and taxonomic levels. Surprisingly however, we found the tympanum of Italian taxa to be smaller than that of H. arborea, while the inverse was initially reported (Nascetti et al., 1995; based on Balkan populations of H. arborea). Given the overlapping morphological variation between these species, we discourage the use of any external criteria for identification purposes. Further investigations should be conducted at the contact zone, where character displacement, especially in breeding calls, could have evolved following reinforcement. More generally, our study should motivate future assessment of the phenotypic variation (notably bioacoustics) between Hyla perrini sp. nov. and H. intermedia s. s. species across their respective ranges, to test whether their genetic differentiation is also adaptive.

We thereby conclude that Italian tree frogs consist of two cryptic species, and describe the northern one, Hyla perrini sp. nov., in the next section. We recommend such multidisciplinary approaches for future taxonomic evaluation of cryptic lineages, as also conducted elsewhere (e.g., Diaz-Rodriguez et al., 2017). Particularly, we stress the need for genomic resources to resolve nuclear divergence and assess patterns of introgression, especially in young sibling species where traditional genetic markers (e.g., allozyme, microsatellite, or few nuclear sequences) might lead to erroneous conclusions.

Description of Hyla perrini sp. nov.

Identity and Diagnosis

Populations of this species have previously been considered as H. intermedia, due to phenotypic similarity and phylogenetic relationships. Based on mtDNA variation, H. perrini sp. nov. was previously called “clade N” by Canestrelli et al. (2007a,b) and “new taxon 2” by Stöck et al. (2008b, 2012). Hyla perrini is a medium-sized tree frog, morphologically similar to the parapatric H. intermedia s. s. (its sister species) and H. arborea. According to the populations analyzed here, it subtly differs from H. intermedia s. s. by a narrower head and shorter hind legs, yet without diagnostic differences. Like most European hylids, H. perrini displays a green or brownish back coloration and an immaculate belly, separated by a lower black and upper white lateral stripe ending near the groin in an inguinal loop of varying shapes. Other typical characteristics include a horizontal pupil, adhesive disks on fingers and toes, and a yellowish vocal sac on the male throat. The advertisement calls have a similar pattern and structure as H. intermedia s. s. and H. arborea: singular notes made of 7–9, partially fused pulses, totaling ~60–80 ms, interrupted by pauses of 80–120 ms, emitted over the 0.5–5 kHz range, with main bursts of energy at 0.95–1.05 kHz (fundamental frequency) and 2.0–2.4 kHz (dominant frequency); a sonogram and a spectrogram of H. perrini from the type locality is provided in Figure 7B, as well as an audio file (Audio S1). To our knowledge, the tadpoles have not been systematically compared but are supposedly similar between these species (Nöllert and Nöllert, 2003): spiracle sinistral; cloaca dextral; pointy tail tip and high upper fin reaching the eyes, with a single line of dark patches on the tail muscle; Labial Tooth Row Formula (LTRF) 2/3. Hyla perrini has evolved since the upper Pliocene, 2.5–3.5 Mya, and can be distinguished from its relatives, notably H. intermedia s. s., by concordant differences at genetic markers, i.e., ~9% of mitochondrial divergence (COI, cyt-b), and ~0.3% of nuclear divergence (RAD tags).

Distribution

Hyla perrini is endemic to the Po Plain and adjacent valleys in Northern Italy and the Swiss canton of Ticino (Figure 3). Its range borders H. arborea in northeastern Italy/Slovenia, where the two species meet in a narrow contact zone along the Isonzo River (Verardi et al., 2009); it borders and hybridizes with H. intermedia s. s. on the northeastern edge of the Apennines in Emilia Romagna. A now-extinct population was introduced in the 1950s from Ticino to the eastern shores of Lake Geneva (Dufresnes et al., 2015c; Dubey et al., 2018). The species is found from sea level up to about 1,500 m in Italy and 1,100 m in Switzerland.

Holotype

MZL42340, adult male collected by Christophe Dufresnes on April 7th 2018 at Piazzogna, canton of Ticino, Switzerland (46.1361°N, 8.8206°E), subsequently deposited at the Cantonal Museum of Zoology of Lausanne, Switzerland (Figure 8). Full measurements are available in Table S3. SVL of 43.3 mm, head narrower than body, with a rounded snout. Tympanum conspicuous, circular, with diameter 55% of the eye diameter. Nostrils slightly protuberant, located at the same distance to the snout than to the eye. Vomerine teeth present, at the same level than choanae. Long legs, 1.4 × the size of the body. Five long, half-webbed, toes, bearing oval discs; relative lengths from inner to outer toes: 4 > 3 > 5 > 2 > 1; protruding metatarsal tubercle, rounded; subarticular tubercles distinct. Arms slender, with four unwebbed fingers bearing oval discs; subarticular tubercles distinct; thenar tubercle present, elongated. Ventral skin granular; vocal sac folds on the throat; dorsal skin smooth. Coloration in life: dorsum olivaceous green, ventrum white, gular folds of the vocal sac yellowish; prominent black and white lateral stripe running from the nostril to the groin, forming a well-marked inguinal loop. Changes of coloration in alcohol: dorsum bluish gray, vocal sac whitish, black lateral stripe diffused.

FIGURE 8

Figure 8. The holotype of Hyla perrini sp. nov. in life and deposited at the Cantonal Museum of Zoology of Lausanne, under voucher MZL42340.

Paratypes

The type series includes four paratypes, two males (MZL42341 and MZL42342) and two females (MZL42343 and MZL42344), collected on April 7th 2018 and May 8th 2018, respectively, by the same collector at the same locality. Their measurements and photos can be found in Table S3 and Figure S5.

Variation

Males body size (SVL) averaged 40.7 mm (35–47 mm) from the two populations of H. perrini surveyed here (Swiss Canton of Ticino and Lombardy, Italy, Table S1D). Rosso et al. (2004) reported mean SVL of 30.7–38.4 mm across five populations of this species. Females are generally larger in tree frogs; few females measured during our 2018 fieldwork were as large as 55 mm SVL. As in related species, back coloration was found to vary from light green to brownish gray, with the tint physiologically changing between environments (CD pers. obs.). The lateral stripe greatly differs between individuals and could be used for individual recognition in mark-capture-recapture studies (e.g., Pellet et al., 2007, on H. arborea). The inguinal loop, especially, can take various shapes, from weak and interrupted dots to heavily marked and elongated commas.

Karyotype

Like other European hylids, Hyla perrini has two sets of 12 chromosomes (2n = 24). A karyotype and idiogram, most likely from this species (specimen K51), is shown by Anderson (1986; p. 61, as “H. arborea”; cf. Anderson, 1991), who analyzed eight metaphases of voucher HGD74640 (Florida State Museum), listed as “? Po River” (p. 347) yet elsewhere detailed “unknown locality in Italy”. Importantly, Anderson (1986) found a secondary constriction in the long arm of chromosome 10 (Anderson, 1986: Figure 2; see Anderson, 1991), interpreted as the position of the nucleolous organizing region (NOR), which is similar to that in the closely related H. arborea (Schmid, 1978). In H. perrini, sex determination is under genetic control by chromosome 1 in an XY manner (LG1, Stöck et al., 2011; Dufresnes et al., 2015b; as H. intermedia), most likely at the dmrt1 gene (Brelsford et al., 2016).

Nomenclatural History

The nomenclatural history of the new species is tightly intertwined with that of H. intermedia Boulenger (1882, p. 381) first described Hyla arborea var. intermedia based on two syntypes: BMNH 52.12.11.55 from Bononia (i.e., Bologna, at the south-eastern edge of the range of the new species), and BMNH 82.7.13.6 from Palermo, Sicily. However, by lectotype designation of BMNH 82.7.13.6, Dubois (1995) made this latter specimen the only name bearer of the epitheton intermedia, and also restricted its type locality to “Palermo,” making intermedia inapplicable to North-Italian tree frogs. In the same work, Dubois (1995) synonymized H. intermedia (i.e., the South-Italian taxon) with H. arborea (Linnaeus, 1758). Using allozymes, Nascetti et al. (1995) found tree frogs from the entire Italian Peninsula to be distinct from H. arborea, and described them as Hyla italica. However, they conducted two additional nomenclatural acts. First, they stated for another old name coined for Italian tree frogs, Hyla variegata (Rafinesque, 1814; type specimens not designated or known to exist; cf. Frost, 2018), an homonym of Rana variegata (Bonnaterre, 1789), hence making it unavailable. Second, they designated a neotype (MZUF 20365) with its type locality near Palermo (Sicily; ca. 37° 00′N, 14° 20′E) and called it H. italica as a replacement name (nomen novum) for H. variegata. By this action and type locality, H. italica became a junior subjective synonym of H. intermedia (as stated by Dubois, 1995), and thus inapplicable to our new species.

Therefore, to the best of our knowledge, there has been no available scientific name for North Italian tree frogs (see also Frost, 2018). In this region, several other authors (cf. Bonato et al., 2007) have previously used either H. arborea or, like De Betta (1857, p. 265), Hyla viridis (= Rana arborea Linnaeus, 1758; cf. Stejneger, 1907), for the herpetofauna of Veneto and Tirol Provinces. Phylogenetic studies did not apply any name other than H. intermedia (Canestrelli et al., 2007a,b). Stöck et al. (2008b, 2012) considered the northern taxon (preliminary called “Hyla new taxon 2”; “n. t. 2”), as a subspecies of H. intermedia (Stöck et al., 2008b, 2012).

To conclude, the name Hyla intermedia applies for the southern lineage (H. intermedia s. s. in this paper), and we have hereby proposed a new name for the tree frogs from Northern Italy.

Etymology

We name this new species Hyla perrini after Nicolas Perrin, professor of ecology, who recently retired from the University of Lausanne (Switzerland). He made countless contributions to the field of evolutionary biology, especially on Palearctic tree frogs, and led numerous studies on phylogeography, population genetics and sex chromosome evolution in this species group, several of them involving this taxon (as H. intermedia or H. n. t. 2). Beyond his scientific excellence, we acknowledge his exceptional mentoring abilities and dedicated support throughout many years of great collaborations.

As a vernacular name in english, we propose H. perrini to be called the “Po's tree frog,” which reflects its geographic distribution and translates as “Raganella padana” (Italian), “Rainette du Po” (French), “Po-Laubfrosch” (German), “Рахкавка з По/Rakhkavka z Po” (Ukrainian), “Квакшаиз По/Kvaksha iz Po” (Russian), “Rela do Pó” (Portuguese), “Trädgroda från Po” (Swedish), “Rainette du Pôô” (Vaudois), and “ yes /Po cheong-keguri” (Korean).

Natural History and Threats

Hyla perrini is a thermophilic amphibian breeding in April-May in shallow standing waters such as small temporary ponds or flooded meadows. Its lifestyle is globally similar to related European Hyla for which detailed information are available in naturalist publications (e.g., Nöllert and Nöllert, 2003; Dufresnes, 2019, where it is treated). In the latest assessments, Italian tree frogs (H. intermedia s. l.) are currently listed as “Least Concern” across the whole range (Andreone et al., 2009) as well as across Italy (Rondinini et al., 2013). Separate re-evaluations for each taxon will be needed, especially for H. perrini which might suffer from agriculture pollution and habitat loss across the heavily-impacted Po Plain, its main range. This new species appears on the Swiss red list of amphibians as “Endangered” (Schmidt and Zumbach, 2005), with relatively stable populations but a restricted distribution in the country (TM pers. comm.). For the time being, we recommend H. perrini to receive the IUCN red list status “Data Deficient”, pending further specific assessment.

Ethics Statement

This study was carried our in accordance of the recommendations of the Ufficio della natura e del paesaggio (Ticino, Switzerland) and the Direzione generale per la protezione della natura e del mare (Italy). The protocol was approved by these committees.

Author Contributions

CD, NR, and DJ designed the study. CD, GM, NR, AlB, DC, and DJ conducted fieldwork. RSe conducted labwork. CD, GM, AlB, SL, SD, and DJ analyzed the data. CB-C, OB, AmB, EC, GF, KG, CG, AH, GL, JL, TM, BP, PS, RSa, and MS contributed samples and ideas. CD drafted the main part of the manuscript, supplemented by SL, GM, MS, and DJ, and subsequently improved by all co-authors.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The first and main person to thank is Nicolas Perrin, for his scientific support and positive influence, which guided our research for years, on Hyla and other systems. We wish him successful fieldtrips to observe the new species. We are also grateful to E. Geffen, R. Rouag, and N. Bruyndonckx for unpublished samples included in the phylogeny; M. Poujai for supplying bioacoustic equipment; Z. Daeppen for help on the field. Olivier Glaizot and André Keiser for their help at the Cantonal Museum of Zoology of Lausanne. Frogs were captured under collecting permits from the Ufficio della natura e del paesaggio to CD (Ticino, Switzerland) and the Direzione generale per la protezione della natura e del mare to DC (Italy). CD was funded by a fellowship from the Swiss National Science Foundation (fellowship P2LAP3_171818).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2018.00144/full#supplementary-material

Figure S1. Localities used in the niche modeling analyses.

Figure S2. PCA on the phylogenomic dataset.

Figure S3. Distribution of morphometric data for H. perrini, H. intermedia s. s. and H. arborea.

Figure S4. Distribution of bioacoustic variables for H. perrini, H. intermedia s. s. and H. arborea.

Figure S5. Ventral and dorsal view of the paratypes of H. perrini.

Table S1. Samples included in this study.

Table S1A. Sequences included in the mitochondrial phylogeny.

Table S1B. Samples included in the nuclear phylogeny.

Table S1C. Localities included in the hybrid zone analyses.

Table S1D. Populations sampled for the biometric analyses.

Table S2. Variables used in the ecological niche modeling and their contributions.

Table S3. Measurements of the type series.

Data Sheet S1. Data files analyzed in this study.

Audio S1. Chorus of H. perrini from the type locality.

References

Anderson, K. (1986). A Cytotaxonomic Analysis of the Holarctic Tree Frogs of the Genus Hyla. Ph.D. dissertation, New York University, New York, NY.

Anderson, K. (1991). “Chromosome evolution in Holarctic Hyla tree frogs,” in Amphibian Cytogenetics and Evolution, eds D. M. Green and S. K. Sessions (New York, NY: Academic Press). 299–322. doi: 10.1016/B978-0-12-297880-7.50017-2

CrossRef Full Text

Andreone, F., Schmidt, B., Vogrin, M., Corti, C., Sindaco, R., and Romano, A. (2009). Hyla Intermedia. The IUCN Red List of Threatened Species 2009:e.T55517A11323168.

Andrews, S. (2010). FastQC: a Quality Control Tool for High Throughput Sequence Data. Babraham Bioinformatics. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc

Arntzen, J. W., McAtear, J., Butôt, R., and Martinez-Solano, I. (2018). A common toad hybrid zone that runs from the Atlantic to the Mediterranean. Amphib. Reptil. 39, 41–50. doi: 10.1163/15685381-00003145