Molecular methods

Hula painted frogs were collected in the Hula Valley during daytime by probing the banks of ponds and sifting through dead vegetation. Captured animals were measured, sampled for DNA and released. Dead individuals were preserved in ethanol for further analyses (Tel Aviv University collection; Am2572-4). DNA was sampled from toe clips of six L. nigriventer, and tissue samples from a single representative of each of the following taxa: D. galganoi, D. jeanneae (here treated as species but often seen as subspecies of D. galganoi14), D. montalentii, D. sardus, D. scovazzi, A. cisternasii, A. dickhilleni, A. maurus, A. muletensis, B. bombina, B. maxima and B. orientalis. DNA was extracted using proteinase K followed by standard phenol/chloroform protocol. We used published primers for three mitochondrial genes (mitochondrial DNA (mtDNA)): 12S rRNA24,25, 16S rRNA24 and cytochrome-b (ref. 26), and three nuclear genes (nuclear DNA (nDNA)): histone H3 (ref. 27), rhodopsin exon 1 (ref. 27) and recombination activating protein 128 (details in Supplementary Table 4). PCR products were sequenced with an ABI 3100 automatic sequencer (Applied Biosystems). Sequence data for all other species were taken from published sources29. Overall, we compiled a total of 2,503 bp for 17 species (Fig. 2), and 1140, bp of mtDNA for 21 species (Supplementary Fig. S1, Supplementary Note 4).

We used maximum likelihood and Bayesian frameworks to reconstruct the phylogenetic position of L. nigriventer after determining from 28 possible substitution models the one best fitting our data (FindModel; http://darwin.uvigo.es/software/modeltest.html) according to the Akaike Information Criterion (AIC). For the nDNA and mtDNA data sets, the best substitution model was General Time Reversible (GTR) plus Gamma (ΔAIC=25.8 between best and second best models). Maximum likelihood trees were inferred in MEGA 5.0530, and 1,000 bootstrap replicates were calculated to assess the support of nodes. For the Bayesian phylogeny, we used MrBayes 3.1 (ref. 31). For all Bayesian analyses, we ran 106 Markov chain Monte Carlo (MCMC) generations and sampled at a frequency of 100 (10,000 samples and 25% burn-in). We ensured that the s.d. of split frequencies was <0.01 and that potential scale reduction factor was about 1.

Alternative analyses partitioning the combined molecular data set by either gene or codon resulted in topologies identical to the unpartitioned analysis and high support for all nodes. Specifically, partitioned Bayesian analyses were performed under two different partition schemes: (1) each gene as a separate partition except for 12S and 16S, which were defined as a single partition and (2) a set of seven partitions defining first, second and third codon positions of the three nuclear genes as three partitions (each merged for the three genes), the three codon positions of the mitochondrial cytochrome-b gene as one partition each, and the merged 12S and 16S sequences as one partition. Substitution models were separately determined for each partition. MrBayes analyses under these settings resulted in topologies identical to those of the unpartitioned analysis and had high support values for all nodes.

Divergence time estimates

Reconstruction of molecular time-trees was based on two different relaxed-clock approaches, either of which allows substitution rates to vary independently over clades (ICR; BEAST, with 50 (ref. 6) generations and a conservative burn-in of 106 generations after examination of convergence of likelihood values and mixing of chains32), or in an autocorrelated fashion (ACR: Multidivtime33). All analyses were carried out separately for the concatenated mtDNA data set and the concatenated nDNA data set, and composite 95% credibility intervals were compiled from the two nDNA analyses34. We used substitution models as in the phylogenetic analyses, except for the nDNA-based BEAST analysis. In this latter analysis, we deviated from the complex GTR+I+G model suggested by AIC. Instead, we applied a simple Jukes–Cantor model to obtain a more realistic time estimate for some of the shallow splits, where only very few nDNA substitutions separated the terminals, and to achieve a better convergence of the posterior values. However, additional runs with GTR+I+G led to similar age estimates for the focal node, that is, the Latonia–Discoglossus split. Preferred analyses ran with all time constraints activated (as reported in results above), but additional runs were carried out to cross-validate these.

The time constraints used for estimating divergence times of Discoglossus frogs were, on one hand, secondary constraints from a previous study35 (with upper and lower bounds from previous 95% credibility intervals), and biogeographic constraints from sister group taxa occurring on opposite sites of the strait of Gibraltar. Time constraints (calibrations) in the following are named C0–C4:

C0: We set the root prior in our analysis (split of Ascaphus from other frogs) at 240 mya35.

C1: We used the constraint of 152 mya (129–179 mya)35 for the split between Bombinatoridae and Alytidae in the form of an upper and lower bound in ACR (129 and 179 mya), and as a prior of 152 mya with normal distribution and a s.d. of 25 mya in ICR.

C2 and C3: Because our combined analysis placed in both the genus Alytes and the genus Discoglossus those species occurring on the opposing sides of the Strait of Gibraltar into monophyletic groups (D. scovazzi versus D. galganoi/jeanneae, and A. maurus versus A. dickhillenii), we concluded that these two species pairs diverged by vicariance after the last direct land connection between Iberia and North Africa had severed, that is, after the Messinian Salinity Crisis at about 5 mya. This is in full agreement with a previous time-tree analysis of Alytes, which inferred a diversification of the A. dickhillenii/maurus/muletensis clade between 2.7 and 5.6 mya36. We therefore constrained the divergence between the two trans-Mediterranean pairs of sister species, each with lower bounds of 5 mya and higher bounds of 7 mya (ACR) and with a normal prior at 5 mya with 2 mya s.d. (ICR).

C4: Because our first exploratory analyses reconstructed unrealistically old ages for the genus Bombina compared with previous studies16, we constrained the split between B. orientalis and Bombina variegata at 5–17 mya in ACR, and with a normal-distributed prior at 10 mya with a large s.d. of 10 mya in ICR35.

Morphological methods

To study Hula painted frog osteology, we used X-ray microtomography and three-dimensional (3D) visualization. Four L. nigriventer specimens, preserved in ethanol, were scanned by a microtomograph (Micro-XCT 400; Xradia, California, USA) at the Weizmann Institute of Science, Israel. One of the specimens was collected in 1955 (HUJ-R-544; National Natural History Collections at the Hebrew University of Jerusalem), and others collected during 2011 (Am2572-4; National Collections of Natural History at Tel Aviv University). Two specimens were fully scanned (HUJ-R-544 and Am2572), and head scans were performed for the other two (Am2573 and Am2574). The samples were sealed in a square poly (methyl methacrylate) box, with 2 mm width windows. At the bottom of the box, a cloth soaked in ethanol ensured a saturated atmosphere around the sample that prevented it from drying out during the scan.

The microtomograph includes a sealed micro-focused source, a detector mounted on a revolving head allowing the choice of linear magnifications and a CCD (charge-coupled device) camera with 2,048 × 2,048 pixels. The micro-computed tomography scans were conducted at 40 kV, 200 μA (8 W of power) and a spatial resolution of 35 μm. The data set for an acquisition consists of nine tomograms centred on different areas of the sample, each of them including 500 projections taken over 180°, with an exposure time of 5 s per projection. Under this setup, at a linear magnification of × 0.5, scanning each individual took about 30 h. The volumes of the nine regions were reconstructed separately with specialized volume reconstruction software (Xradia, California, USA). The 3D images were generated with 16-bit precision and subsequently converted into 8 bit for visualization. Three-dimensional processing and rendering were obtained after semi-automatic segmentation of the skeleton using surface and volume rendering in AVIZO 7.01 (VSG, SAS, Merignac, France; http://www.vsg3d.com).

We studied 22 osteological characters, identifiable on non-articulated specimens of extinct and extant Alytidae taxa (Supplementary Table S2). Character choice was based on Clarke’s19 comparative morphology of extant Discoglossus, with some adjustments made based on personal observations and previous studies18,20,37,38,39.

Character states for the osteological traits were determined by us for extant and fossil L. nigriventer, and based on literature records and museum specimens for the other fossil and recent Alytidae (Eodiscoglossus, Callobatrachus, Wealdenbatrachus, Kizylkuma, Paradiscoglossus, Paralatonia, Latoglossus, Latonia, Discoglossus, Alytes, Bombina and Barbourula; Supplementary Note 1). Because D. montalentii had some distinct characters from the other Discoglossus species, we defined it as separate taxon for phylogenetic analysis, whereas other Discoglossus species were merged because of their osteological similarity. Phylogenetic reconstruction on osteological characters used Bayesian inference31 (MrBayes, 506 generations with 106 burn-ins) and maximum parsimony40 (PAUP 4.0, all characters unordered, treating multiple states in one or several taxa as polymorphism and calculating 1,000 bootstrap replicates).

We also used logistic regression to calculate the probability of presence of an additional set of nine osteological characters common in extant L. nigriventer, extant Discoglossus species, fossil L. nigriventer, Latonia gigantea and Latonia ragei (Supplementary Table S3). The matrix of presence/absence of morphological characters, composed from skeletons of recent species and non-articulated bones of fossil species, was fitted to a fully factorial logistic model with species and character as the main effects.