Understanding the causes underlying heterogeneity of molecular evolutionary rates among lineages is a long-standing and central question in evolutionary biology. Although several earlier studies showed that modern frogs (Neobatrachia) experienced an acceleration of mitochondrial gene substitution rates compared to non-neobatrachian relatives, no further characterization of this phenomenon was attempted. To gain new insights on this topic, we sequenced the complete mitochondrial genomes and nine nuclear loci of one pelobatoid (Pelodytes punctatus) and five neobatrachians, Heleophryne regis (Heleophrynidae), Lechriodus melanopyga (Limnodynastidae), Calyptocephalella gayi (Calyptocephalellidae), Telmatobius bolivianus (Ceratophryidae), and Sooglossus thomasseti (Sooglossidae). These represent major clades not included in previous mitogenomic analyses, and most of them are remarkably species-poor compared to other neobatrachians.
We reconstructed a fully resolved and robust phylogeny of extant frogs based on the new mitochondrial and nuclear sequence data, and dated major cladogenetic events. The reconstructed tree recovered Heleophryne as sister group to all other neobatrachians, the Australasian Lechriodus and the South American Calyptocephalella formed a clade that was the sister group to Nobleobatrachia, and the Seychellois Sooglossus was recovered as the sister group of Ranoides. We used relative-rate tests and direct comparison of branch lengths from mitochondrial and nuclear-based trees to demonstrate that both mitochondrial and nuclear evolutionary rates are significantly higher in all neobatrachians compared to their non-neobatrachian relatives, and that such rate acceleration started at the origin of Neobatrachia.
Through the analysis of the selection coefficient (ω) in different branches of the tree, we found compelling evidence of relaxation of purifying selection in neobatrachians, which could (at least in part) explain the observed higher mitochondrial and nuclear substitution rates in this clade. Our analyses allowed us to discard that changes in substitution rates could be correlated with increased mitochondrial genome rearrangement or diversification rates observed in different lineages of neobatrachians.
It has been long acknowledged that character change in evolution occurs at different rates, which can vary extremely between different lineages . Although initially described for morphological characters, among-lineage rate heterogeneity also occurs at the molecular level, and evolutionary biologists have been long interested in quantifying molecular evolutionary rates as well as determining which are the underlying mechanisms that trigger their acceleration or slowdown in different lineages . However, uncovering the causes of lineage-specific rate variation has proven to be challenging, and previous studies reached different conclusions, which attempted to explain rate heterogeneity through correlation with species body size, generation time, population dynamics, metabolic rates, or habits (e.g., [2–4]). Furthermore, molecular evolutionary rates have also been correlated with diversification [5–7], but given the multiple factors that shape diversification patterns, the generalization of this correlation is elusive, and therefore, the cause-effect between rates of genome evolution and cladogenesis remain largely unknown .
At the molecular level, the fixation of mutations in an evolutionary lineage (i.e., substitution events), is a complex dynamic process determined by the interaction between evolutionary (selection) and demographic (drift) forces . Comparative studies of relative substitution rates have been particularly useful in providing insights into particular constraints of specific genetic systems, such as e.g. mitochondrial (mt) genomes . Mitochondrial DNA has been widely used as a marker in molecular systematics during past decades [10, 11]. As data accumulated, it has become apparent that animal mt DNA evolves at a rate 5 to 10 times faster than single-copy protein coding nuclear genes, although this varies extremely across genes and taxa [10, 12, 13]. Mitochondrial DNA suffers from high mutational pressure  likely due to the inaccuracy of its DNA repair system , the absence of histone-like proteins , the particular replication model with single-strand intermediates , and the presence of reactive oxidative compounds produced in the mitochondria . This high mutational pressure, together with a reduced population size  and the absence of substantial recombination  (but see ), leads to an increase of substitution rate in mt DNA. The comparison of evolutionary rates among lineages permits the identification of events of acceleration and slowdown of rates, which can be further studied to uncover the underlying process (or processes) that produced the observed patterns .
Furthermore, phylogeneticists are particularly interested in understanding evolutionary rate variation because the unequal substitution rates among lineages are a well-known source of phylogenetic artifacts [21, 22]. Rapidly evolving lineages may appear closely related (and often placed close to outgroups) regardless of their true evolutionary relationships (long-branch attraction; ), whereas short branches may also attract each other because of the “leftover” similarity of symplesiomorphic states that “eroded” away in rapid-evolving lineages .
Previous studies have shown that rates of molecular evolution, both for mt DNA and some nuclear genes are unequally distributed among lineages of frogs . Neobatrachian frogs exhibit higher mt substitution rates compared to their non-neobatrachian relatives [25–29]. Yet, it is neither clear when the shifts in substitution rates occurred during the evolutionary history of modern frogs nor whether rate changes are exclusive to the mt genome. Moreover, the heterogeneous distribution of mt substitution rates, together with the high genetic divergence between frogs and their closest living sister taxa (i.e., salamanders; ) are the source of several phylogenetic artifacts in previous studies, such as the monophyly of non-neobatrachian frogs (“Archaeobatrachia”: [29, 31, 32]) or the incorrect phylogenetic placement of Neobatrachia due to long-branch attraction effects . The unequal distribution of mt substitution rates across the anuran tree has also been suggested to yield considerably older time estimates for divergences among neobatrachians .
Neobatrachia (modern frogs) is acknowledged as the most derived lineage of frogs [34, 35] and the sister group of Pelobatoidea [36, 37]. Modern frogs constitute an evolutionarily highly successful clade that contains over 96% of the overall species diversity of extant amphibians [38–40]. Most of this diversity is concentrated in two major clades: Ranoides (= Ranoidea), which comprises three well-supported monophyletic groups (Afrobatrachia, Microhyloidea, Natatanura), and Nobleobatrachia [39, 41]. Neobatrachia` also includes the following species-poor families: Calyptocephalellidae, Heleophrynidae, Limnodynastidae, Myobatrachidae, Nasikabatrachidae, and Sooglossidae, whose relative position is still a contentious issue in anuran phylogeny [39, 41]. Most of the diversity of Ranoides and Nobleobatrachia is located in the Old World and the Neotropics, respectively , whereas the abovementioned species-poor neobatrachian families show a relict distribution . It was suggested that the shift in mt substitution rates in Neobatrachia could be related with the higher diversification rates observed in Ranoides and Nobleobatrachia, provided that further data and analyses could possibly assign short branches to species-poor lineages . To validate such possibility, it is necessary to precisely delimit the node in the anuran phylogeny at which the shift in evolutionary rates took place.
In this study, we newly determined the mt genomes and partial sequences of nine nuclear genes of several key representatives of species-poor neobatrachian families outside Ranoides and Nobleobatrachia. These new sequence data together with previously available orthologous sequence data from other anurans were used to infer a robust phylogeny and a timetree of major lineages of frogs. The new sequence data and the phylogeny were used to (i) confirm that mt substitution rates are (statistically) significantly higher in the different neobatrachian lineages; (ii) localize the rate shift in the phylogeny; (iii) explore whether substitution rates of nuclear genes are also accelerated in modern frogs; and (iv) determine whether changes in substitution rates could be correlated with life-history traits, an increase in mt genome rearrangement, or an increase in diversification rates.
Taxon sampling and DNA sequencing
Taxon sampling in this study was designed to represent all major groups of frogs with particular emphasis on neobatrachian lineages outside Ranoides and Nobleobatrachia. For phylogeny reconstruction, we used Leiopelma and Ascaphus as outgroup taxa. For timetree inference, the outgroup included three salamanders, three caecilians, a lizard, a bird, and a mammal. These extra outgroup taxa were included to provide additional calibration points, and thus to infer more accurate timetree estimates. In-group frog species were largely chosen based on available complete mt genomes, and to allow direct comparison and combination of mt and nuclear data. Available frog mt genomes were expanded with newly determined complete sequences for one pelobatoid (Pelodytes sp., from the southwest of the Iberian Peninsula, probably representing an undescribed species currently under study; in this paper, we name it P. punctatus following current taxonomy) and the following neobatrachians: Heleophryne regis (Heleophrynidae), Lechriodus melanopyga (Limnodynastidae), Calyptocephalella gayi (Calyptocephalellidae), Telmatobius bolivianus (Ceratophryidae), Sooglossus thomasseti (Sooglossidae). The nearly complete sequence of another Sooglossidae, Sooglossus sechellensis, was also determined. A nuclear DNA data set of partial sequences of nine protein-coding genes was compiled, including the recombination-activating genes 1 and 2 (rag1 and rag2), brain-derived neutrotrophic factor (bdnf), proopiomelanocortin (pomc), chemokine receptor type 4 (cxcr4, exon 2), members 1 and 3 of the solute carrier family 8 (slc8a1, exon 2, and slc8a3, exon2), rhodopsin (rho exon 1), and histone 3 (h3a). The nuclear matrix was generated expanding a recent dataset  with new data for the aforementioned species. Since our phylogenetic analyses were focused mainly on the family level or above, and in order to maximize the completeness of the nuclear data set, sequences from congeneric anuran species (for which there is strong evidence for the monophyly of the genus) were merged in few cases see (e.g.,[42–44]). Similarly, chimerical sequences were also used to represent some non-anuran major evolutionary lineages in the outgroup of the timetree analysis. Detailed information on the studied species and the corresponding GenBank accession numbers can be found in Additional file 1.
Total DNA was prepared from ethanol-preserved muscle tissue by proteinase k digestion, phenol-chloroform extraction, and ethanol precipitation . The complete mt genome of Pelodytes was amplified in several overlapping fragments by PCR using the primers and conditions reported in . The remaining mt genomes, corresponding to neobatrachians, were partially amplified using the same set of primers (from 5’-trnF to 3’-cox3) (abbreviations of mt genes follow ). Due to the gene rearrangements found in neobatrachians (see Results and discussion), the remaining halves of the mt genomes (from 5’-cox3 to 3’-trnF) were amplified using the primers and conditions reported in . Partial sequences of nuclear genes were amplified using the primers and conditions reported in the literature: rag1, rag2[25, 49], slc8a1, bdnf and pomc, rho, h3a. In all cases, PCR cycling conditions were experimentally adjusted from those reported in the original publications. Specific primers were designed when general primers did not work (mainly for control regions), and to sequence long-PCR products by primer walking (available from authors upon request). PCR reactions of fragments up to 1500 bp were carried out with 5PRIME Taq DNA polymerase (5PRIME GmbH, Hamburg, Germany), and longer fragments were amplified using LA Taq polymerase (TaKaRa Bio Inc., Otsu, Shiga, Japan), following manufacturer’s instructions. PCR amplicons were purified by ethanol precipitation  or directly from electrophoresis gels using the Speedtools PCR clean-up kit (Biotools B&M Labs. S.A., Madrid, Spain). The long-PCR products containing the control region of L. melanopyga and T. bolivianus were digested with PstI at 37°C for four hours, obtaining two fragments from each of the original amplicons. These four fragments as well as all other PCR products containing the control regions of the remaining species were cloned into pGEM-T vectors (Promega, Madison, WI, USA). PCR fragments and positive recombinant clones were cycle-sequenced with the ABI Prism BigDye Terminator v3.1 cycle sequencing ready reaction kit (Applied Biosystems, Foster City, CA, USA) using PCR and M13 universal primers, and following manufacturer’s instructions. Cycle sequencing products were run on ABI Prism 3700 and 3130xl DNA Analyzers (Applied Biosystems, Foster City, CA, USA).
The new mt sequences were annotated by comparison with other reported frog mt genomes using DOGMA . In this web-based tool, genes are identified by BLAST searches, open reading frames of protein-coding genes are translated using the appropriate genetic code (vertebrate mt code), and transfer RNA (tRNA) genes are further identified based on their putative cloverleaf secondary structure. The gene arrangements of the new mt genomes were compared against the Mitozoa database release 7.1 .
Mitochondrial and nuclear protein-coding genes were analyzed at the nucleotide and amino acid levels. Alignments were generated using TranslatorX . First, amino acids of deduced proteins were aligned using MAFFT L-INS-i  and default settings. Then, ambiguously aligned positions were removed using Gblocks, v.0.19b  and the following settings: minimum number of sequences for a conserved position 29, minimum number of sequences for a flanking position 29, maximum number of contiguous non-conserved positions 8, minimum length of an initial block 5, minimum length of a block 5, allowed gap positions with half. Finally, trimmed protein alignments were used as guide for a codon-based alignment of nucleotide sequences. Mitochondrial tRNA gene nucleotide sequences were aligned manually based on their putative secondary structure, whereas mt ribosomal RNA (rRNA) gene nucleotide sequences were aligned with MAFFT L-INS-i  and corrected by eye for any obvious misalignment. Ambiguously aligned positions in both mt tRNA and rRNA gene alignments were excluded with Gblocks v.0.19b  using the following settings: minimum number of sequences for a conserved position 31, minimum number of sequences for a flanking position 36, maximum number of contiguous non-conserved positions 5, minimum length of a block 10, allowed gap positions with half.
Previous studies (e.g.,) showed long-branch attraction artifacts in phylogenetic reconstruction of the anuran tree due to high rates of evolution of mt DNA in neobatrachians. In order to avoid such phylogenetic inference artifacts, especially those caused by possible saturation in the fast-evolving lineages, first and third codon positions of mt protein-coding genes, and third codon positions of nuclear protein-coding genes were excluded from the corresponding nucleotide alignments. For phylogenetic and dating analyses, protein-coding gene alignments at the nucleotide or amino acid level were combined together with tRNA and rRNA genes into two matrices, hereafter the combined nucleotide and amino acid data sets, respectively.
The combined nucleotide and amino acid data sets were analyzed by maximum likelihood (ML; ) using RAxML v.7.0.4 , and by Bayesian inference (BI; ) using MrBayes v.3.1.2 [60, 61]. ML searches used the rapid hill-climbing algorithm  starting from 100 randomized maximum-parsimony trees. RAxML optimized parameters of the GTR+I+Γ model in all partitions independently. For BI, two independent runs were performed, each with 4 simultaneous Markov chains for 20 million generations, sampling every 1000 generations. Convergence was checked a posteriori using Tracer v.1.5  and the online tool AWTY . The first 1 million generations were discarded as burn-in to prevent sampling before the Markov chains reached stationarity. Support for internal branches was evaluated performing 1000 replicates of non-parametric bootstrapping  (ML) and by posterior probabilities (BI). The genera Leiopelma and Ascaphus were used as outgroups, as they are confidently identified as the sister taxa of all other extant frogs [26, 39, 41, 66, 67].
The Akaike information criterion (AIC; ) was used to select the best partition schemes as well as best-fit nucleotide and amino acid models for each partition. The best partitioning scheme for the combined nucleotide data set was determined with PartitionFinder  and included five partitions: (i) second codon positions of all mt protein-coding genes, (ii) mt rRNA genes, (iii) mt tRNA genes, (iv) first codon positions of all nuclear genes, and (v) second codon positions of all nuclear genes. For the combined amino acid data set, partitions had to be determined separately for amino acid (PartitionFinderProtein; ) and nucleotide (PartitionFinder) data. The best partition scheme favored individual protein-coding gene partitions, a single mt rRNA partition and a single mt tRNA gene partition.
Estimation of divergence times
We used BEAST v.1.6.1  to estimate divergence times among major frog lineages based on molecular data. This program implements a Bayesian dating method, and assumes a relaxed uncorrelated clock in which the rate for each branch is drawn independently from an underlying lognormal distribution . We used the combined nucleotide data set, and constrained the tree topology to the best ML tree (Figure 1) by removing the operators that act on tree topology from the .xml file. The Yule process  was used to describe cladogenesis, and independent GTR+I+Γ models were applied for each of the five data partitions. The final Markov chain was run twice for 100 million generations, sampling every 10,000 generations and the first 1 million was discarded as part of the burn-in process, according to the convergence of chains checked with Tracer v.1.5. . The effective sample size of all the parameters was above 200 .
Seven calibration points were used as priors for divergence times of certain splits, using a lognormal distribution of prior probability. Calibration points were chosen based on previous literature and the online resource Lisanfos KMS v.1.2  that compiles data on amphibian fossils. Fossils provided hard minimum bounds (offset) and mean and standard deviations (SD) were chosen so that the 95% probability limit corresponds to a soft maximum bound. Details on fossil dates and prior distribution parameters for each calibration point are provided in Additional file 2.
Evolutionary rate heterogeneity analyses
To examine among-lineage rate heterogeneity patterns, we reconstructed phylogenetic trees based on mt and nuclear data, separately, and compared their topological congruence and branch length patterns. In these comparative analyses, we were not interested in recovering the species phylogeny (as above) but in evidencing the empirical differences in evolutionary rates between neobatrachian and non-neobatrachian lineages as inferred based on the two types of molecular markers. Therefore, all codon positions of protein-coding genes were included in the analyzed nucleotide data matrices.
Phylogenetic analyses were based on the concatenation of the nucleotide sequences of all mt single-gene alignments (hereafter the mt nucleotide data set) or all nuclear single-gene alignments (hereafter the nuclear nucleotide data set), as well as on the same alignments but with open reading frames translated into proteins (hereafter the mt amino acid and the nuclear amino acid data sets). Trees were reconstructed under ML as explained in the phylogenetic reconstruction section above. Protein-coding genes were partitioned by codon position in the mt and nuclear nucleotide data sets, or by gene in the mt and nuclear amino acid data sets. The mt nucleotide and amino acid data sets had two additional partitions corresponding to combined rRNA and tRNA genes, respectively.
In order to compare substitution rates of mt and nuclear genes, relative-rate tests (RRTs; ) were performed with the program RRTree . This program extends the method of Li and Bousquet  and compares mean rates between lineages relative to the outgroup, taking phylogenetic relationships into account by topological weighting . RRTs were performed for each single-gene alignment (note that all tRNA genes, due to their short length, were concatenated, and here considered like a single-gene alignment for analytic purposes), as well as for the mt and nuclear data sets at both the nucleotide and amino acid levels. In the case of protein-coding genes at the nucleotide level, tests were carried out (i) taking into account all codon positions and (ii) without first and third-codon positions of mt genes and third positions of nuclear genes. Genetic distances were estimated with K2P  and Poisson  models at the nucleotide and amino acid levels, respectively. We defined the three salamander species used in the timetree analysis as outgroup, and divided frogs into three assemblages: species-rich clades within Neobatrachia (Ranoides and Nobleobatrachia), species-poor neobatrachian lineages (Heleophryne, Calyptocephalella, Lechriodus and Sooglossus), and non-neobatrachian relatives (Amphicoela, Discoglossoidea, Pipoidea, and Pelobatoidea). RRTs were performed among these three assemblages, as well as between all neobatrachians versus non-neobatrachians.
Branch length comparisons
We optimized model parameters and branch lengths separately for the mt and nuclear data sets both at the nucleotide and amino acid levels using RA × ML v.7.0.4 , and constraining the topology to the preferred ML tree as recovered based on the combined nucleotide and amino acid data sets. In order to compare the branch-specific bias of mt versus nuclear branch lengths, the ratios between branch lengths (mt / nuclear) were calculated for each individual internal and terminal branch in the nucleotide- and amino acid-based trees, separately. The estimated ratios were subjected to one-way ANOVAs, after being log-transformed to meet the assumptions of normality and homogeneity of variance. The first ANOVAs compared neobatrachians against non-neobatrachians. For the second ANOVAs, neobatrachians were further divided into species-rich and species-poor lineages (following the scheme of RRT analyses; see above), and orthogonal contrasts were used to compare the three groups. All statistical analyses were performed with IBM SPSS Statistics, release 22.214.171.124.
Detecting changes in selective pressure
Simulation studies have shown that analyses of selection coefficients are rather robust to sequence divergence  (as is the case in the present study), having been successfully used in various studies with highly divergent species (e.g.,). In order to understand whether acceleration of evolutionary rates in neobatrachians is due to changes in selective pressure, we tested alternative models with different assumptions about ratios of non-synonymous/ synonymous substitution rates (ω). The software PAML v.3.15  was used to estimate the likelihood and the ω values of different models derived from the preferred topology (Figure 1) and sequence information from single-gene alignments with all codon positions, as well as the mt and nuclear nucleotide data sets. Branch lengths were first optimized for each data set assuming a single ω for the whole tree, and they were fixed when all other parameters were estimated under alternative models. The null model had a single ω value for all branches, and it was compared against four alternatives, which allowed a second ω value on (i) the stem branch of Neobatrachia, (ii) all neobatrachian branches, (iii) Ranoides, or (iv) Nobleobatrachia (including their stem branch). Given that the alternative hypotheses nest the null model, a likelihood ratio test (LRT) was used to determine their significance . To gain insight into the obtained results, we additionally allowed ω to vary on main non-neobatrachian lineages (including their stem branch): (v) Amphicoela, (vi) Discoglossoidea, (vii) Pipoidea, (viii) Pelobatoidea, or (ix) the stem branch of Pelobatoidea. These additional five models were compared against the null model by LRT, and, in addition, all 10 (non nested) models were compared simultaneously using the AIC .
Functional analysis of neobatrachian amino acid synapomorphies
Ancestral character states were reconstructed using MrBayes v3.1.2  based on single-gene protein alignments, for (i) Neobatrachia, (ii) its closer sister group (Pelobatoidea), and (iii) their last common ancestor. The three hypothetical ancestral sequences were compared against each other to identify synapomorphic amino acid changes in Pelobatoidea and Neobatrachia, taking only into account the sites with reliably reconstructed states (we empirically found that p > 0.75 offered a good balance between the number of predictions and their corresponding reliability). A two-sided binomial test was used to assess the asymmetrical distribution of molecular synapomorphies between the two clades. In addition, to further understand if molecular synapomorphies of neobatrachians (or pelobatoideans) were associated to particular regions of the proteins, we predicted the accessibility to solvent and the occurrence at the different trans-membrane regions for each of the identified synapomorphic sites. Solvent accessibility was calculated through BLAST searches against the PDBFINDER2 database , and trans-membrane helices of genes were predicted with TMHMM v.2.0 .
The newly determined mt (JF703228-34) and nuclear (JF703235-51) sequences are available at NCBI (http://www.ncbi.nlm.nih.gov/genbank/). The data sets used in this study (combined nucleotide, combined amino acid, mt nucleotide, mt amino acid, nuclear nucleotide, and nuclear amino acid data sets) as well as the .xml file used for the BEAST analyses can be accessed in the Dryad Digital Repository (doi:10.5061/dryad.3qd54).
Results and discussion
New mitochondrial genomes and nuclear data
We newly determined the complete nucleotide sequence of the light strand of the mt genome of one pelobatoid (Pelodytes punctatus; JF703231) and five neobatrachian species (Helephryne regis, JF703229; Lechriodus melanopyga, JF703230; Calyptocephalella gayi, JF703228; Telmatobius bolivianus, JF703234; Sooglossus thomasseti, JF703233), as well as the nearly complete mt genome of Sooglossus sechellensis (JF703232). We also determined partial sequences of several nuclear genes (see Additional file 1 for details) in order to construct a nuclear data matrix that complemented the mt sequence data set. Main structural features of the newly sequenced mt genomes are highlighted below, and other features can be found in Additional file 2.
The gene order of the mt genome in Pelodytes punctatus follows the consensus of vertebrates and other reported pelobatoids [27, 47] (Figure 1). Calyptocephalella gayi, Telmatobius bolivianus, Sooglossus thomasseti and S. sechellensis conform to the consensus mt gene order of neobatrachians (Figure 1). The neobatrachian and vertebrate consensus mt gene orders differ in the relative position of the trnL(CUN), trnT and trnP genes, which in neobatrachian mt genomes are found next to the trnF gene (forming the LTPF tRNA cluster) downstream the control region  (Figure 1). The mt genome of Lechriodus melanopyga follows the consensus neobatrachian gene order with the only exception of the location of the putative origin of replication of the light strand in a 218 bp-long intergenic spacer between the genes trnY and cox1 (Figure 1). Heleophryne regis departs from the consensus order of neobatrachians in two regions. The trnM gene is pseudogenized (Figure 1; ΨM) in its ancestral location (IQM tRNA gene cluster) because the anticodon has a deletion. The functional trnM appears within the WANCY tRNA gene cluster, which is rearranged as ANCMWY, without changes in the coding strands (Figure 1). The putative origin of the replication of the light strand (Figure 1; OL) is located in a 165 bp-long intergenic spacer between trnW and trnY genes (Figure 1). Interestingly, the two described new mt gene rearrangements are associated with origins of replication, which are considered hot spots for gene order change in the vertebrate mt genome [26, 86, 87].
These two newly reported gene rearrangements are consistent with the tandem duplication-random loss model [88, 89], which is considered the main mechanism of gene order change in vertebrate mt genomes [26, 87]. The tandem duplication-random loss model is reinforced by the presence of the trnM pseudogene, which remains in the ancestral location of trnM in H. regis (Figure 1). Other trnM pseudogenes have been reported in mt genomes of several members of the frog family Mantellidae [90, 91], and several fishes including parrot fishes of the family Scaridae , Carapus bermudensis (Ophidiiformes; ), and Diaphus splendidus (Myctophiformes; ). It has been suggested that these pseudogenes are maintained because they are needed for the posttranscriptional processing of the nad2 gene [91, 92].
The mt gene order of Heleophryne regis is key for understanding the evolution of mt genome rearrangements in neobatrachians given that this species is basal in this clade (see reconstructed phylogeny in Figure 1). The mt gene order of H. regis has the characteristic LTPF tRNA gene cluster, which can, thus, be confidently considered a molecular synapomorphy of all neobatrachians. We speculate that the long intergenic spacer between trnP and trnF genes found in H. regis could be a remnant of the ancestral tandem duplication and random loss event by which the trnL(CUN), trnT and trnP genes moved to the LTPF tRNA cluster in the origin of Neobatrachia .
Due to reported long-branch attraction effects in reconstructed anuran phylogenies based on mt sequence data , the concatenated data set of all mt and all nuclear single-gene alignments was trimmed to only retain most conserved positions (combined nucleotide data set: 11,136 positions). Based on the combined nucleotide data set with five partitions, both ML (-lnL=76,155.99) and BI (-lnL=76,547.00 for run 1; -lnL=76,548.29 for run 2) arrived at an identical topology (Figure 1), which is our best hypothesis for frog phylogenetic relationships. Phylogenetic analyses based on the combined amino acid data set (data not shown) also recovered the same topology (Figure 1). Five major clades of frogs were recovered with high support: non-neobatrachian lineages branched off successively as (i) Amphicoela (Ascaphus + Leiopelma), (ii) Discoglossoidea, (iii) Pipoidea and (iv) Pelobatoidea, which was the sister group of (v) Neobatrachia, in agreement with recent molecular studies [36, 37, 39, 66]. Internal relationships within Discoglossoidea , Pipoidea [36, 37, 39, 66], and Pelobatoidea  fully agreed with recent morphological and molecular analyses.
Neobatrachia has traditionally been acknowledged as monophyletic [96, 97], a fact that has been corroborated by most morphological [34, 35] and molecular [36, 37, 39] studies, including the present work (Figure 1). The reconstructed tree recovered Heleophryne as sister group to all other neobatrachians, as reported by most recent molecular studies [25, 41, 66], and this placement received moderate support from ML and maximal from BI (Figure 1). The Australasian Lechriodus and the South American Calyptocephalella were recovered together in the same clade, and both as sister group to Nobleobatrachia (Figure 1). The two Seychellois Sooglossus species were grouped together as sister group of Ranoides (as early suggested by Savage ) (Figure 1). Although support for this relationship was maximal for BI, it received only moderate support for ML. Notably, other recent molecular studies failed to resolve the relative position of sooglossids with confidence [25, 39, 41, 66].
Estimation of divergence times
The two independent BEAST analyses gave very similar estimates of divergence times for each node (mean difference between runs was 0.638 ± 0.687 million years), and these estimates roughly agreed with other recent studies of divergence times among anurans [66, 67, 99]. The obtained estimates of divergence among amphibians were especially in agreement with those of a recent BEAST analysis , despite the differences in taxon sampling, choice of molecular markers, and calibration points. The origin of crown-group Anura was inferred in the Middle Triassic (about 230 mya) (Figure 2), and the initial diversification of non-neobatrachian frogs (successive branching of Amphicoela, Discoglossoidea, Pipoidea, and Pelobatoidea) in the Late Triassic – Early Jurassic (about 210–192 mya) (Figure 2), confirming other recent molecular dating studies [43, 99]. The split between Neobatrachia and Pelobatoidea was dated in the Late Triassic – Early Jurassic, before the initial break-up of Pangaea (mean 192 mya; 95% CI 219–166), and the initial neobatrachian diversification that led to Ranoides and Nobleobatrachia in the Late Jurassic – Early Cretaceous (160–130 mya) (Figure 2) in agreement with divergence time estimates of other recent studies [33, 67, 99–101].
Our estimates were younger than those inferred by earlier studies [67, 99, 100] that used MultiDivtime [102, 103], even though the 95% credibility intervals (CI) mostly overlapped. These discrepancies could be due to differences of both programs in methodological assumptions of rate change (auto-correlated in MultiDivtime, uncorrelated in BEAST), implementation of evolutionary models and prior calibrations, and techniques to calculate credibility intervals . Moreover, molecular dating estimates are much older than the first neobatrachian fossils (dated in the Early Cretaceous; ), indicating either that currently available fossils might be a poor indicator of this particular branching event  or that the molecular dating could be overestimated .
Congruence between mt- and nuclear-based phylogenies
Separate phylogenetic analyses of the mt (13,580 positions) and nuclear (7,083 positions) nucleotide data sets rendered two almost identical topologies (see Figure 3) with levels of support only slightly lower than those obtained in the phylogenetic analysis based on the combined nucleotide and amino acid data sets (not shown). The recovered phylogenetic tree based on the mt nucleotide data set (Figure 3a) placed sooglossids as the most basal neobatrachians (branching off before Heleophryne), which is likely an artifact due to the attraction of the extremely long branch leading to sooglossids towards that of Heleophryne and the stem branch of Neobatrachia (which is the longest branch in the tree) . Additionally, this phylogenetic tree (Figure 3a) favored with low statistical support alternative relationships for (Duttaphrynus + (Telmatobius + Hyla)) within Nobleobatrachia. Interestingly, the phylogenetic analysis of the mt amino acid data set differed from our best hypothesis (Figure 1) only in the basal position of sooglossids within Neobatrachia (not shown). The phylogenies reconstructed based on the nuclear nucleotide (Figure 3b) and nuclear amino acid (not shown) data sets recovered our best hypothesis for frog phylogenetic relationships (Figure 1).
In addition to assessing topological congruence, the separate analyses of the mt and nuclear data sets (both at nucleotide and amino acid levels) offer further information on the mode of evolution of these two different genetic systems. The comparison of mt- and nuclear-based trees (Figure 3) reveals that branch lengths in the mt tree (Figure 3a), which ultimately correspond to the underlying substitution rates, are on average 3.29 (0.42 – 12.15) and 2.82 (0.31-18.19) times longer than those in the nuclear tree (Figure 3b) at the nucleotide and amino acid levels, respectively. This confirms previous evidence of higher substitution rates in the mt DNA of metazoans [10, 12, 13]. More importantly, neobatrachians exhibit much longer branches in the mt trees compared to their non-neobatrachian relatives, with Heleophryne, Sooglossus, and natatanuran species having the longest branches. In contrast, neobatrachians do not display such conspicuous long branches in the nuclear tree, and a lineage-specific pattern of branch lengths is only subtle (Figure 3b).
Relative-rate tests and comparison of branch lengths
In order to test whether the mt substitution rate is significantly higher in neobatrachians, and to study putative lineage-specific rate changes in nuclear genes, we conducted RRTs, as well as a direct comparison of the ratios of the lengths of the same branch in the mt- and nuclear-based phylogenies (both at the nucleotide and amino acid levels). Overall, RRTs clearly showed that neobatrachians had significantly higher mt substitution rates compared to non-neobatrachians (K2 versus K1; Table 1). RRTs were applied to 16 individual mt gene alignments (13 protein-coding genes, 2 rRNA genes, and a single alignment combining all tRNA genes) at the nucleotide level, and neobatrachians had higher mean relative rates in all cases but the nad3 gene (Table 1). Differences were significant (p < 0.05) for 12 out of 16 mt genes, and highly significant (p < 0.001) for the concatenation of all mt genes (mt nucleotide data set; Table 1). RRTs of individual nuclear genes at the nucleotide level showed that six out of nine genes had higher mean substitution rates for neobatrachians (although significant only in the case of the genes rag2 and slc8a1). Genes h3a, pomc, and rho showed lower rates for Neobatrachia, but differences were non-significant (Table 1). Notably, the concatenation of all nuclear genes (nuclear nucleotide data set; Table 1) gave significantly higher rates for neobatrachians than for non-neobatrachians.
Results from relative-rate tests based on nucleotide single-gene alignments and mt and nuclear nucleotide data sets
Mean weighted substitution rates
Probability of relative-rate tests
p (1 vs. 2) p < 0.05
p (1 vs. 3) p < 0.0167
p (1 vs. 4) p < 0.0167
p (3 vs. 4) p < 0.0167
all mt genes
all nuc genes
K values are mean weighted substitution rates for (1) non-neobatrachians, (2) all neobatrachians, (3), species-poor neobatrachians (Heleophrynidae, Calyptocephalellidae, Limnodynastidae, Sooglossoidea), and (4) species-rich neobatrachians (Ranoides and Nobleobatrachia). Probabilities (p) of relative-rate tests correspond to the comparisons between the different groups (shown in parentheses). Statistically significant results (p < 0.05; or p < 0.05/3 = 0.0167 to correct for multiple testing) are in bold italics and marked with an asterisk. Abbreviations of mt genes follow ; see main text for full names of nuclear genes.
Additional RRTs based on the same nucleotide data sets mentioned above were performed to discriminate among three alternative hypotheses regarding the origin of the higher mt rates of neobatrachians: (i) mt substitution rates became accelerated in the origin of the clade and were maintained higher ever since; (ii) the higher substitution rates are specific to species-poor neobatrachian lineages; and (iii) the higher substitution rates are specific to species-rich neobatrachian lineages. In the second and third hypotheses, rate acceleration would be either punctual at the base of Neobatrachia (and thus, not present in species-rich lineages), or it would be a specific feature of the species-rich Ranoides and Nobleobatrachia (as suggested by ), respectively. The comparison of species-rich (Ranoides, Nobleobatrachia) and species-poor (Heleophryne, Calyptocephallela, Lechriodus, Sooglossus) neobatrachian groups against non-neobatrachian relatives (K4 versus K1, and K3 versus K1, respectively; Table 1) showed that both groups had consistently higher mt substitution rates, which were significant after Bonferroni correction for multiple comparisons (p < 0.05/ 3 = 0.0167). No significant differences in mt substitution rates were observed between species-rich and species-poor neobatrachians (K3 versus K4; Table 1). However, for the concatenation of all mt genes (mt nucleotide data set), species-poor neobatrachian lineages showed overall significantly higher rates compared to species-rich neobatrachians. In the case of the different nuclear genes, the separate comparison of species-rich and species-poor neobatrachians against non-neobatrachians gave similar results: the rag2 gene had consistently higher mean rates for both neobatrachian groups, but the higher neobatrachian rates found for the slc8a1 gene failed to be significant due to the lower significance threshold. No differences in nuclear substitution rates were observed among neobatrachians (K3 versus K4; Table 1). The concatenation of all nuclear genes (nuclear nucleotide data set) rendered non-significant comparisons (Table 1). All previously mentioned RRTs were also performed with protein-coding genes translated into amino acids, and produced consis-tent results with those obtained based on nucleotides (see Additional file 3).
In order to disentangle the contribution of the different codon positions to evolutionary rate acceleration in neobratrachians, we performed additional RRTs based on the mt (lacking first and third codon positions) and nuclear (lacking third codon positions) portions of the combined nucleotide data set, respectively (see Additional file 3). As expected, exclusion of most variable positions from the analyses produced lower mean weighted substitution rates in both mt and nuclear genes. This would in part explain the better performance of the combined nucleotide data set in resolving anuran phylogeny. Nevertheless, RRT results showed that conserved codon positions also contribute significantly to neobatrachian-specific rate acceleration, displaying similar trends to those experienced by more variable codon positions (see Additional file 3).
In order to compare the branch-specific rate bias of mt versus nuclear genes and quantify the acceleration of mt rates in Neobatrachia, we calculated the ratio of the lengths of each branch in the mt- and nuclear-based trees (both at the nucleotide and amino acid levels). For the ANOVAs, the ratios were log-transformed to meet the assumptions of normality (Shapiro-Wilk’s test on residuals p > 0.05) and homogeneity of variance (Levene’s test p > 0.05). A first set of ANOVAs supported a significant difference of the ratios of all neobatrachian versus all non-neobatrachian branches (p = 0.012 and 0.001 for nucleotide- and amino acid-based trees, respectively). The higher mt / nuclear ratios of neobatrachians (mean ± standard deviation for nucleotide- and amino acid-based trees, respectively: 3.90 ± 2.38 and 3.90 ± 3.67) compared to those of non-neobatrachians (2.66 ± 2.66 and 1.67 ± 1.27), showed that relative to the nuclear genome, the mt genome is approximately 46 and 130% more accelerated in Neobatrachia at the nucleotide and amino acid levels, respectively. The striking difference in percentage estimates between nucleotide and amino acid levels could be explained by saturation of more variable codon positions in the former. A second set of ANOVAs found significant differences in branch length ratios between species-poor neobatrachians, species-rich neobatrachians, and non-neobatrachians (p = 0.039 and 0.002 for nucleotide and amino acids, respectively). Further orthogonal contrasts found significant differences when non-neobatrachians were compared against species-poor (p = 0.027 and 0.001; nucleotide and amino acid levels, respectively) and species-rich (p = 0.049 and 0.011) neobatrachians. However, no significant differences were found between neobatrachian groups (p = 0.650 and 0.291).
Estimation of branch length ratios based on the mt (lacking first and third codon positions) and nuclear (lacking third codon positions) portions of the combined nucleotide data set rendered the same patterns of rate heterogeneity as derived from comparative analyses based on all codon positions and amino acids. Interestingly, however, when considering only conserved codon positions, we inferred that the mt genome is approximately 88% more accelerated in Neobatrachia relative to the nuclear genome. This percentage is closer to that obtained based on amino acid data, and further supports that the one based on all positions may be underestimated due to saturation.
Overall, RRTs support a statistically significant acceleration of the mt substitution rate at the origin of neobatrachians, which is shared by both species-poor and species-rich lineages. This result not only corroborates previous studies that suggested an unequal distribution of mt substitution rates between non-neobatrachian and neobatrachian frogs [25–28, 33], but also characterizes the evolutionary dynamics of the shift in rates within neobatrachians. The inferred evolutionary pattern is further corroborated by ANOVA results indicating that neobatrachian mt lineages are evolving 88 to 130% faster than non-neobatrachian mt lineages with respect to nuclear rates. Nuclear genes also exhibited an overall significant substitution rate acceleration, although much lower in absolute terms than the one experienced by the mt genome. In fact, the nuclear pattern was not as evident as that of the mt genes, since neobatrachian-specific higher substitution rates were found to be significant in only two out of nine nuclear genes. This might be the result of the different properties of mt and nuclear genomes, such as the recombination rate (very limited in mt DNA; ) or the different effective population size (smaller in mt DNA; ), which can influence the effectiveness of selection upon these two genetic systems . An alternative explanation is that although substitution rate acceleration is general for both mt and nuclear genomes, our results are biased by the use of particular nuclear genes, which cannot represent the complexity of the entire nuclear genome, with genes obviously subjected to very disparate selective regimes .
Changes in selective pressure and molecular synapomorphies
The assumption of a single selective coefficient (ω) for the frog tree (Figure 1) (null model) rendered ω values well below 1 for all different mt and nuclear genes (0.005-0.16) (Table 2), indicating the action of purifying selection to maintain gene function . To understand whether observed acceleration of the mt substitution rate in neobatrachians could be due to changes in selective pressure, and to compare the strength of selection acting on the mt and nuclear genomes, we tested for putative changes in the selective regime under four different scenarios for the Neobatrachia. All outcomes are available in the Additional file 3, and main results are here highlighted.
Results from branch models that assume branch-specific changes in the selection coefficient (ω) in Neobatrachia, based on nucleotide single-gene alignments and mt and nuclear nucleotide data sets
all mt genes
all nuc genes
The table shows ω values estimated for the whole frog tree (null model) and values estimated for specific (Branch) and remaining background (Backgr.) branches under the four alternative models tested. Bold italics highlight results that are significantly different to the null model (LRT p < 0.05). See Additional file 3 for full details.
For all mt genes, the independent ω values inferred for the stem branch of Neobatrachia were always higher than those estimated when a single ω was assumed for the frog tree (null model) (Table 2). These differences were only significant (LRT p < 0.05) for the genes cob, cox1, cox3, and nad1, as well as for the combination of all mt genes (Table 2). The independent ω estimates of mt genes under alternative models for “all Neobatrachia”, “all Ranoides” and “all Nobleobatrachia” branches were generally higher than those of the null model, but unlike the model of the stem branch Neobatrachia, ω was not higher for every mt gene, and fewer individual genes gave significant results (3, 1, and none for the three alternative models, respectively; see Table 2). These outcomes suggest that purifying selection acting on mt proteins could have been relaxed in neobatrachians.
In order to understand the relative support of the first four models tested (within Neobatrachia), and to further investigate the relevance of the obtained results, we compared all 10 models using the AIC. For the combination of all mt genes, the model of relaxed selection in the stem branch of Neobatrachia was clearly better than the remaining models (Table 3). All other models showed a difference of AIC values (ΔAIC) higher than 10: ΔAIC=48 for the second best model (independent ω for Pipoidea), ΔAIC=56 for the third (independent ω shared by all neobatrachian branches), ΔAIC=61 for the fourth (independent ω for Ranoides), etc. (Table 3). A notable exception to the above pattern was the gene cox1; despite evidence of relaxed selection in the stem branch of Neobatrachia, the comparison of all 10 models for cox1 strongly favored the relaxation along all branches of Neobatrachia (ΔAIC to the rest of models was always >10, and up to 44; Table 3).
Comparison of all 10 branch models that assume branch-specific changes in the selection coefficient (ω) (including the null model), based on nucleotide single-gene alignments and the mt and nuclear nucleotide data sets
“ All Discogl”
all mt genes
all nuc genes
The table shows the differences in AIC values among all 10 models.
For nuclear genes, most of the estimated independent ω values in all of the nine alternative models were lower than those of the null model, showing evidence of stronger purifying selection, although genes did not display a concordant pattern neither under particular models nor within specific genes (Table 2). However, there were two exceptions: (i) under the “all Neobatrachia” alternative model, relaxation of purifying selection on nuclear genes was frequently recovered (in six out of nine genes, even though it was statistically significant only for the genes rag1, rho, slc8a1, and the combination of all nuclear genes; LRT p<0.05); (ii) under the model of an independent ω for Amphicoela, for which relaxation of selection was also frequent (all nuclear genes except h3a, although the differences were only statistically significant for the genes pomc, slc8a1, slc8a3, and the combination of all nuclear genes) (Table 2). Using the concatenation of all nuclear genes, the comparison of the models assuming a second independent ω for (i) all Neobatrachia and (ii) Amphicoela favored the latter (ΔAIC=14; Table 3).
Overall, our results evidence the relaxation of purifying selection acting on mt DNA in Neobatrachia. Among all tested models, the assumption of relaxation at the stem branch leading to Neobatrachia outperformed the rest, although the very similar results obtained under the model of a general relaxation along the entire Neobatrachia indicates that this alternative hypothesis cannot be confidently rejected. In any case, these results suggest that overall relaxation in selection pressure could be, at least in part, responsible for the general acceleration of mt substitution rates at the origin of Neobatrachia, as found by RRTs and topological measures. Interpreting the results from nuclear genes was more complex because inferred relaxed selection along all neobatrachian branches could only explain the higher substitution rates found by RRTs for the gene slc8a1. Moreover, it is important to note that results derived from the comparison of different selection regimes should be taken with caution, because analyzed sequences are highly divergent and silent substitutions might be saturated, thus compromising the correct estimation of ω values .
Most identified amino acid synapomorphies corresponded to Neobatrachia in mt proteins (102 versus 49), whereas distribution of synapomorphies in nuclear proteins was only slightly higher in neobatrachians (24 versus 22). These differences were only significant (binomial test’s p < 0.05) for genes cox1, nad5 and rag2 (see Additional file 3). This is in agreement with the results of the RRTs, which revealed overall higher substitution rates in neobatrachians, and a more pronounced acceleration in mt genes. Most mt synapomorphies corresponded to leucine, serine, and alanine in both Neobatrachia and Pelobatoidea (18, 13, and 11 for Neobatrachia; and 10, 10, and 7 for Pelobatoidea, respectively). In the nuclear genes, the most frequent synapomorphies for Neobatrachia were serine (4), glutamic acid (3), and lysine (3), whereas for Pelobatoidea, they were leucine (6) and aspartic acid (3) (see Additional file 3). To further understand how proteins of neobatrachians could have accommoda-ted the corresponding mutations, we investigated whether synapomorphic amino acids showed any particular pattern of exposition to solvent, or whether they were associated with specific domains of trans-membrane proteins. However, distribution of neobatrachian synapomorphic changes were not related apparently to these functional traits, with mutations being distributed in a more or less uniform manner along mt proteins.
Could substitution rates be associated with life history traits, rates of diversification or mt gene rearrangements in frogs?
Our analyses recovered a fully resolved and robust phylogeny of frogs after new data on key lineages of Neobatrachia were added. These new lineages were essential to understand the origin of the observed higher substitution rates in neobatrachians. RRTs and branch length measures demonstrate the presence of a significant acceleration in both mt and nuclear substitution rates in Neobatrachia, which is shared by both species-rich and species-poor neobatrachian lineages.
The ultimate causes for among-lineage evolutionary rate variation are, in general, rather elusive. Several lines of evidence suggest that particular life-history traits may be responsible for rate variation , and at least, three main hypotheses have been put forward. (i) According to the generation time hypothesis, species with shorter generation times are expected to have higher substitution rates because their genomes are copied more often per time unit . However, estimating generation time is often difficult, and thus, the age at sexual maturity and the time to first reproduction are used as proxies [109, 110]. In frogs, generation time is generally short, and small- to medium-sized frogs typically reach sexual maturity in their first or second year of life . Sexual maturity one year after egg-laying is common in many anuran species, including both neobatrachian and non-neobatrachian frogs (e.g., Bombina, Discoglossus, Hymenochirus or Xenopus). On the other hand, several neobatrachians have longer generation times, such as Heleophryne (larval period of over 2 years; ), or Anaxyrus canorus (Bufonidae), whose females reach sexual maturity after 4–6 years . Data on sexual maturity available from the AnAge database (build 12; ) does not support the existence of significant differences in age of sexual maturity between neobatrachian and non-neobatrachian frogs, either for males or females (Student’s t p = 0.530 and p = 0.754, respectively). However, these results should be taken with caution, as the available number of age estimates of sexual maturity data in AnAge is still insufficient to draw robust conclusions (N between 32–34 and 2–3 for neobatrachians and non-neobatrachian males and females, respectively).
(ii) The longevity hypothesis proposes that long-lived or late reproducing species will have lower rates of molecular evolution , as they are expected to have more effective DNA repairing mechanisms . In most cases, maximum longevity data is derived from captive specimens, which have higher life expectancy that in the wild . Available data from the AnAge database suggests that short- and long-lived species are present both among neobatrachian and non-neobatrachian frogs. The comparison of longevity between neobatrachians and non-neobatrachians did not render significant differences (Student’s t p = 0.066; N = 10 and 75 for neobatrachians and non-neobatrachians, respectively). Inferences derived from this incomplete data set have to be taken with caution, though.
(iii) The metabolic rate hypothesis holds that rates of evolution are correlated with the production of free radicals during respiration . Metabolic rate is correlated with substitution rates in the frog family Dendrobatidae , but comparable data for other anuran groups is currently unavailable. A suitable proxy for metabolic rate might be genome size, which has often been shown to be inversely proportional to metabolic rate [118–120]. However, the comparison of genome sizes among anurans (data from the Animal genome size database; ) did not provide indications for consistently larger genome sizes in non-neobatrachians compared to neobatrachian frogs (Student’s t p = 0.770): C-values (mean ± standard deviation, minimum-maximum in parentheses) were 5.10 ± 3.65 pg (1.29-11.38 pg; N=9) and. 5.43 ± 2.59 pg (1.40-13.40 pg; N=27), respectively.
Some studies found a correlation of high substitution rates with more events of gene order rearrangements in the mt genome of some metazoans [122, 123]. However, frogs do not conform to this pattern, because the mt substitution rate became accelerated in the origin of Neobatrachia, and this acceleration is not exclusive of Natatanura, where most mt gene rearrangements are found in frogs [90, 124–127]. Furthermore, Kurabayashi et al.  found evidence for the absence of correlation between mt rates and number of gene rearrangements in one intensively studied lineage of neobatrachians (mantellid frogs from Madagascar).
Many other studies have found substitution rates to be correlated with species diversification [5–7, 128], and three main hypotheses have been proposed to explain this correlation . (i) Speciation is often associated with processes that can potentially increase substitution rates, such as adaptation to new environments or transient reductions in population sizes that reduce the efficiency of purifying natural selection [129, 130]. (ii) Higher substitution rates could produce higher net diversification, both by increasing speciation rate and/ or by reducing extinction rate . (iii) A third hypothesis rejects a causal relationship between substitution and diversification rates, and holds that this correlation is due to other factors that influence both rates simultaneously .
In frogs, it has been suggested that observed higher mt substitution rates of neobatrachians could be the product of faster recent speciation events in this clade (including more bottleneck events) . In addition, Dubois  hypothesized that direct-developing species (mostly within Neobatrachia) would tend to have higher substitution rates, and this would in turn promote speciation. Alternatively, it has also been suggested that higher substitution rates of neobatrachians could be responsible for higher diversification rates, due to shorter generation times and/or higher metabolic rates . Both species-rich and species-poor lineages of neobatrachians share higher mt substitution rates compared to non-neobatrachian relatives, but species diversity is highly unequally distributed among them, with most of the diversity corresponding to Ranoides and Nobleobatrachia . Moreover, neobatrachians do not fit the hypothesis of the different reproductive modes  because, most direct-developers belong to species-rich clades, whereas species-poor neobatrachian lineages are mostly indirect-developers [34, 132], but both groups share the presence of higher mt substitution rates. Therefore, the relationship between higher mt substitution and diversification rates in frogs remains elusive, and unless a rampant extinction of neobatrachian lineages external to currently species-rich clades (Ranoides and Nobleobatrachia) could account for the observed huge differences in diversity, it can be considered that substitution and diversification rates are decoupled in frogs. Unfortunately, the paucity of the current fossil record may hinder the answer to this question [105, 133, 134].
Using both complete mt genomes and partial sequences of nine nuclear loci, we inferred a robust phylogeny of frogs. RRTs and branch length measures found compelling evidence of higher substitution rates in the mt genome of neobatrachian frogs, and a subtle (but significant) trend in nuclear genes. Phylogenetic analyses suggest that the origin of this rate acceleration began at the stem branch leading to Neobatrachia, in the Early-Middle Jurassic period. Because substitution rates are determined, to a great extent, by the balance between selection and genetic drift , we studied the changes in selective pressure in frogs, and found that purifying selection acting on most mt and some nuclear proteins might have been relaxed in Neobatrachia. Therefore, we suggest that this relaxation of purifying selection could explain, at least in part, the general rate acceleration observed in this group.
We do not exclude the possibility that our results are slightly affected by some sort of phylogenetic artifact , but the neobatrachian-specific higher mt substitution rates are reinforced by compelling evidence of relaxed purifying selection on mt proteins. Furthermore, our results show that selection might have been relaxed also in nuclear genes, and thus justify the higher substitution rates found in the genes rag2 and slc8a1 in neobatrachians. Data from additional nuclear genes, which are likely to be gathered soon in the context of genome sequencing initiatives, hold the key to confirm or reject a putative general acceleration of evolutionary rates in neobatrachian frogs. Likewise, the clarification of the causes that relaxed purifying selection would need further, in-depth studies that investigate intrinsic and extrinsic factors that might have modified the fitness landscape of gene function .
With the exception of few particular linages (e.g., ), available data on life history traits for frogs is generally scarce and not representative for the main lineages within Anura. Available data suggests that no clear differences exist between neobatrachian and non-neobatrachian frogs with respect to generation time, longevity, or metabolic rate, but more data would be necessary to reliably test if substitution rates could be correlated with particular life history traits in frogs.
We are indebted to Renaud Boistel, Marius Burger, Alan Channing, Bill Love, and Stephen Richards for tissue samples or help in the field. We thank David Buckley and Silvia Perea for support with the BEAST software, and to Mario García-París for access to computational resources. We are grateful, to David Buckley, Mario García-París, Patrick S. Fitze, F. Borja Sanchíz, Luis M. San-Jose, Enrico Negrisolo, and four anonymous reviewers for insightful comments that improved earlier versions of the manuscript. This work was supported by the Spanish Ministry of Science and Innovation (CGL2010-18216 to RZ and RYC-2011-09321 to DSM), and the Consejo Superior de Investigaciones Científicas of Spain (CSIC) and the European Social Fund (ESF) (JAE-pre PhD fellowhip to II). We acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).
Department of Biodiversity and Evolutionary Biology, Museo Nacional de Ciencias Naturales, CSIC
Department of Animal Biology, University of Barcelona
Structural Computational Biology Group, Spanish National Cancer Research Center (CNIO)
UMR 72205 OSEB, Département de Systématique et Evolution, Muséum National d’Histoire Naturelle
Division of Evolutionary Biology, Zoological Institute, Technical University of Braunschweig
Simpson GG: Tempo and mode in evolution. New York: Columbia University Press; 1944.
Bromham L: Why do species vary in their rate of molecular evolution?Biol Lett 2009, 5:401–404.PubMedView Article
Smith SA, Donoghue MJ: Rates of molecular evolution are linked to life history in flowering plants.Science 2008, 322:86–89.PubMedView Article
Santos JC: Fast molecular evolution associated with high active metabolic rates in poison frogs.Mol Biol Evol 2012, 29:2011–2018.View Article
Eo SH, DeWoody JA: Evolutionary rates of mitochondrial genomes correspond to diversification rates and to contemporary species richness in birds and reptiles.Proc R Soc B 2010, 277:3587–3592.PubMedView Article
Lanfear R, Ho SYW, Love D, Bromham L: Mutation rate is linked to diversification in birds.Proc Natl Acad Sci USA 2010, 107:20423–20428.PubMedView Article
Barraclough TG, Savolainen V: Evolutionary rates and species diversity in flowering plants.Evolution 2001, 55:677–683.PubMedView Article
Bromham L: Putting the 'bio' into bioinformatics.Biol Lett 2009, 5:391–393.PubMedView Article
Mueller RL: Evolutionary rates, divergence dates, and the performance of mitochondrial genes in Bayesian phylogenetic analysis.Syst Biol 2006, 55:289–300.PubMedView Article
Moritz C, Dowling TE, Brown WM: Evolution of animal mitochondrial DNA: relevance for population biology and systematics.Ann Rev Ecol Syst 1987, 18:269–292.View Article
Cummings MP, Meyer A: Magic bullets and golden rules: data sampling in molecular phylogenetics.Zoology 2005, 108:329–336.PubMedView Article
Meyer A: Evolution of mitochondrial DNA in fishes. In Biochem Mol Biol Fishes Volume 2. Hochachka PW, Mommsen TP: Elsevier Science Publishers B. V.; 1993:1–38.
Brown WM, George MJ, Wilson AC: Rapid evolution of animal mitochondrial DNA.Proc Natl Acad Sci USA 1979, 76:1967–1971.PubMedView Article
Castellana S, Vicario S, Saccone C: Evolutionary patterns of the mitochondrial genome in Metazoa: exploring the role of mutation and selection in mitochondrial protein coding genes.Genome Biol Evol 2011, 3:1067–1079.View Article
Bogenhagen DF: Repair of mtDNA in vertebrates.Am J Hum Genet 1999, 64:1276–1281.PubMedView Article
Reyes A, Gissi C, Pesole G, Saccone C: Asymmetrical directional mutation pressure in the mitochondrial genome of mammals.Mol Biol Evol 1998, 15:957–966.PubMedView Article
Richter C, Park JW, Ames BN: Normal oxidative damage to mitochondrial and nuclear DNA is extensive.Proc Natl Acad Sci USA 1988, 85:6465–6467.PubMedView Article
Lynch M: The origins of genome architecture. Sunderland, Massachusetts: Sinauer Associates, Inc.; 2007.
Brown WM: Evolution of animal mitochondrial DNA. In Evolution of genes and proteins. Edited by: Nei M, Koehn RK. Sunderland, Massachusetss: Sinauer Associates; 1983:62–88.
Tsaousis AD, Martin DP, Ladoukakis ED, Posada D, Zouros E: Widespread recombination in published animal mtDNA sequences.Mol Biol Evol 2005, 22:925–933.PubMedView Article
Philippe H, Germot A: Phylogeny of eukaryotes based on ribosomal RNA: long-branch attraction and models of sequence evolution.Mol Biol Evol 2000, 17:830–834.PubMedView Article
Rodríguez-Ezpeleta N, Brinkmann H, Roure B, Lartillot N, Lang FB, Philippe H: Detecting and overcoming systematic errors in genome-scale phylogenies.Syst Biol 2007, 56:389–399.PubMedView Article
Felsenstein J: Cases in which parsimony or compatibility methods will be positively misleading.Syst Zool 1978, 27:401–410.View Article
Fuellen G, Wägele J-W, Giegerich R: Minimum conflict: a divide-and-conquer approach to phylogeny estimation.Bioinformatics 2001, 17:1168–1178.PubMedView Article
Hoegg S, Vences M, Brinkmann H, Meyer A: Phylogeny and comparative substitution rates of frogs inferred from sequences of three nuclear genes.Mol Biol Evol 2004, 21:1188–1200.PubMedView Article
Irisarri I, San Mauro D, Green DM, Zardoya R: The complete mitochondrial genome of the relict frog Leiopelma archeyi: Insights into the root of the frog Tree of Life.Mitochondrial DNA 2010, 21:173–182.PubMedView Article
Gissi C, San Mauro D, Pesole G, Zardoya R: Mitochondrial phylogeny of Anura (Amphibia): a case study of congruent phylogenetic reconstruction using amino acid and nucleotide characters.Gene 2006, 366:228–237.PubMedView Article
San Mauro D, García-París M, Zardoya R: Phylogenetic relationships of discoglossid frogs (Amphibia: Anura: Discoglossidae) based on complete mitochondrial genomes and nuclear genes.Gene 2004, 343:357–366.PubMedView Article
Hay JM, Ruvinsky I, Hedges SB, Maxson LR: Phylogenetic relationships of amphibian families inferred from DNA sequences of mitochondrial 12S and 16S ribosomal RNA genes.Mol Biol Evol 1995, 12:928–937.PubMed
Zardoya R, Meyer A: On the origin of and phylogenetic relationships among living amphibians.Proc Natl Acad Sci USA 2001, 98:7380–7383.PubMedView Article
Hedges SB, Maxson LR: A molecular perspective on lissamphibian phylogeny.Herpetol Monogr 1993, 7:27–42.View Article
Dutta SK, Vasudevan K, Chaitra MS, Shanker K, Aggrwal RK: Jurassic frogs and the evoution of amphibian endemism in the Western Ghats.Curr Sci 2004, 86:211–216.
Igawa T, Kurabayashi A, Usuki C, Fujii T, Sumida M: Complete mitochondrial genomes of three neobatrachian anurans: a case study of divergence time estimation using different data and calibration settings.Gene 2008, 407:116–129.PubMedView Article
Duellman WE, Trueb L: Biology of amphibians. New York: MacGraw-Hill; 1986.
Ford LS, Cannatella DC: The major clades of frogs.Herpetol Monogr 1993, 7:93–117.View Article
Roelants K, Bossuyt F: Archaeobatrachian paraphyly and Pangaean diversification of crown-group frogs.Syst Biol 2005, 54:111–126.PubMedView Article
Irisarri I, Vences M, San Mauro D, Glaw F, Zardoya R: Reversal to air-driven sound production revealed by a molecular phylogeny of tongueless frogs, family Pipidae.BMC Evol Biol 2011, 11:114.PubMedView Article
Pyron AR, Wiens JJ: A large-scale phylogeny of Amphibia including over 2800 species, and a revised classification of extant frogs, salamanders, and caecilians.Mol Phylogenet Evol 2011, 61:543–583.PubMedView Article
Frost DR, Grant T, Faivovich J, Bain RH, Haas A, Haddad CFB, de Sá RO, Channing A, Wilkinson M, Donnellan SC, Raxworthy CJ, Campbell JA, Blotto BL, Moler P, Drewes RC, Nussbaum R, Lynch JD, Green DM, Wheeler WC: The amphibian tree of life.Bull Amer Mus Nat Hist 2006, 297:1–370.View Article
Crottini A, Madsen O, Poux C, Strauß A, Vieites DR, Vences M: Vertebrate time-tree elucidates the biogeographic pattern of a major biotic change around the K-T boundary in Madagascar.Proc Natl Acad Sci USA 2012, 109:5358–5363.PubMedView Article
San Mauro D: A multilocus timescale for the origin of extant amphibians.Mol Phylogenet Evol 2010, 56:554–561.PubMedView Article
Burleigh JG, Mathews S: Assessing among-locus variation in the inference of seed plant phylogeny.Int J Plant Sci 2007, 168:111–124.View Article
Sambrook J, Fritsch EF, Maniatis T: Molecular cloning: a laboratory manual. Cold Spring Harbor, New York: Cold Spring Harbor Laboratory Press; 1989.
San Mauro D, Gower DJ, Oommen OV, Wilkinson M, Zardoya R: Phylogeny of caecilian amphibians (Gymnophiona) based on complete mitochondrial genomes and nuclear RAG1.Mol Phylogenet Evol 2004, 33:413–427.PubMedView Article
Boore JL: Animal mitochondrial genomes.Nuc Acids Res 1999, 27:1767–1780.View Article
Kurabayashi A, Sumida M: PCR Primers for the neobatrachian mitochondrial genome.Curr Herpetol 2009, 28:1–11.View Article
Nylander JAA, Wilgenbusch JC, Warren DL, Swofford DL: AWTY (are we there yet?): a system for graphical exploration of MCMC convergence in Bayesian phylogenetics.Bioinformatics 2008, 24:581–583.PubMedView Article
Felsenstein J: Confidence limits on phylogenies: an approach using the bootstrap.Evolution 1985, 39:783–791.View Article
Roelants K, Gower DJ, Wilkinson M, Loader S, Biju SD, Guillaume K, Moriau L, Bossuyt F: Global patterns of diversification in the history of modern amphibians.Proc Natl Acad Sci USA 2007, 104:887–892.PubMedView Article
San Mauro D, Vences M, Alcobendas M, Zardoya R, Meyer A: Initial diversification of living amphibians predated the breakup of Pangaea.Am Nat 2005, 165:590–599.PubMedView Article
Akaike H: Information theory and an extension of the maximum likelihood principle. In Second international symposium of information theory. Edited by: Petrov BN, Csaki F. Budapest: Akademiai Kiado; 1973.
Lanfear R, Calcott B, Ho SYW, Guindon S: PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses.Mol Biol Evol 2012, 29:1695–1701.PubMedView Article
Sarich VM, Wilson AC: Generation time and genomic evolution in primates.Science 1973, 179:1144–1147.PubMedView Article
Robinson-Rechavi M, Huchon D: RRTree: relative-rate tests between groups of sequences on a phylogenetic tree.Bioinformatics 2000, 16:296–297.PubMedView Article
Li P, Bousquet J: Relative rate test for nucleotide substitutions between two lineages.Mol Biol Evol 1992, 9:1185–1189.
Robinson M, Gouy M, Gautier C, Mouchiroud D: Sensitivity of the relative-rate test to taxonomic sampling.Mol Biol Evol 1998, 15:1091–1098.PubMedView Article
Kimura M: A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences.J Mol Evol 1980, 16:111–120.PubMedView Article
Yang Z: Computational molecular evolution. New York: Oxford University Press; 2006.View Article
Buschiazzo E, Ritland C, Bohlmann J, Ritland K: Slow but not low: genomic comparisons reveal slower evolutionary rate and higher dN/dS in conifers compared to angiosperms.BMC Evol Biol 2012, 12:8.PubMedView Article
Yang Z: PAML: a program package for phylogenetic analysis by maximum likelihood.Comput Appl Biosci 1997, 13:555–556.PubMed
Anisimova M, Bielawski JP, Yang Z: Accuracy and power of the likelihood ratio test in detecting adaptive molecular evolution.Mol Biol Evol 2001, 18:1585–1592.PubMedView Article
Hooft RWW, Sander C, Scharf M, Vriend G: The PDBFINDER database: a summary of PDB, DSSP and HSSP information with added value.Comput Appl Biosci 1996, 12:525–529.PubMed
Krogh A, Larsson B, von Heijne G, Sonnhammer ELL: Predicting transmembrane protein topology with a hidden Markov model: application to complete genomes.J Mol Biol 2011, 305:567–580.View Article
Sumida M, Kanamori Y, Kaneda H, Kato Y, Nishioka M, Hasegawa M, Yonekawa H: Complete nucleotide sequence and gene rearrangement of the mitochondrial genome of the Japanese pond frog Rana nigromaculata.Genes Genet Syst 2001, 76:311–325.PubMedView Article
Mindell DP, Sorenson MD, Dimcheff DE: Multiple independent origins of mitochondrial gene order in birds.Proc Natl Acad Sci USA 1998, 95:10693–10697.PubMedView Article
San Mauro D, Gower DJ, Zardoya R, Wilkinson M: A hotspot of gene order rearrangement by tandem duplication and random loss in the vertebrate mitochondrial genome.Mol Biol Evol 2006, 23:227–234.PubMedView Article
Moritz C, Brown WM: Tandem duplications of D-loop and ribosomal RNA sequences in lizard mitochondrial DNA.Science 1986, 233:1425–1427.PubMedView Article
Boore JL: The duplication/random loss model for gene rearrangement exemplified by mitochondrial genomes of deuterostome animals. In Comparative genomics. Edited by: Sankoff D, Nadeau JH. Kluwer Academic Publisher; 2000:133–147.View Article
Kurabayashi A, Sumida M, Yonekawa H, Glaw F, Vences M, Hasegawa M: Phylogeny, recombination, and mechanisms of stepwise mitochondrial genome reorganization in Mantellid frogs from Madagascar.Mol Biol Evol 2008, 25:874–891.PubMedView Article
Kurabayashi A, Usuki C, Mikami N, Fujii T, Yonekawa H, Sumida M, Hasegawa M: Complete nucleotide sequence of the mitochondrial genome of a Malagasy poison frog Mantella madagascariensis: evolutionary implications on mitochondrial genomes of higher anuran groups.Mol Phylogenet Evol 2006, 39:223–236.PubMedView Article
Mabuchi K, Miya M, Satoh TP, Westneat MW, Nishida M: Gene rearrangements and evolution of tRNA pseudogenes in the mitochondrial genome of the Parrotfish (Teleostei: Perciformes: Scaridae).J Mol Evol 2004, 59:287–297.PubMedView Article
Miya M, Takeshima H, Endo H, Ishiguro NB, Inoue JG, Mukai T, Satoh TP, Yamaguchi M, Akira K: Major patterns of higher teleostean phylogenies: a new perspective based on 100 complete mitochondrial DNA sequences.Mol Phylogenet Evol 2003, 26:121–138.PubMedView Article
Miya M, Kawaguchi A, Nishida M: Mitogenomic exploration of higher teleostean phylogenies: a case study for moderate-scale evolutionary genomics with 38 newly determined complete mitochondrial DNA sequences.Mol Biol Evol 2001, 18:1993–2009.PubMedView Article
García-París M, Buchholz DR, Parra-Olea G: Phylogenetic relationships of Pelobatoidea re-examined using mtDNA.Mol Phylogenet Evol 2003, 28:12–23.PubMedView Article
Duellman WE: On the classification of frogs.Occ Pap Mus Nat Hist Univ Kansas 1975, 42:1–14.
Reig OA: Proposiciones para una nueva macrosistemática de los anuros (nota preliminar).Physis 1958, 21:109–118.
Savage JM: The geographic distribution of frogs: Patterns and predictions. In Evolutionary biology of the anurans: Contemporary research on major problems. Edited by: Vial JL. Columbia: Columbia University Press; 1973:351–445.
Roelants K, Haas A, Bossuyt F: Anuran radiations and the evolution of tadpole morphospace.Proc Natl Acad Sci USA 2011, 108:8731–8736.PubMedView Article
Zhang P, Zhou H, Chen Y-Q, Liu Y-F, Qu L-H: Mitogenomic perspectives on the origin and phylogeny of living amphibians.Syst Biol 2005, 54:391–400.PubMedView Article
Vences M, Vieites DR, Glaw F, Brinkmann H, Kosuch J, Veith M, Meyer A: Multiple overseas dispersal in amphibians.Proc R Soc B 2003, 270:2435–2442.PubMedView Article
Thorne JL, Kishino H: Divergence time and evolutionary rate estimation with multilocus data.Syst Biol 2002, 51:689–702.PubMedView Article
Thorne JL, Kishino H, Painter IS: Estimating the rate of evolution of the rate of molecular evolution.Mol Biol Evol 1998, 15:1647–1657.PubMedView Article
San Mauro D, Agorreta A: Molecular systematics: a synthesis of the common methods and the state of knowledge.Cell Mol Biol Lett 2010, 15:311–341.PubMedView Article
Báez AM, Moura GJB, Gómez RO: Anurans from the lower cretaceous crato formation of northeastern Brazil: implications for the early divergence of neobatrachians.Cretaceous Res 2009, 30:829–846.View Article
Rodríguez-Trelles F, Tarrío R, Ayala FJ: A methodological bias toward overestimation of molecular evolutionary time scales.Proc Natl Acad Sci USA 2002, 99:8112–8115.PubMedView Article
Comeron JM, Williford A, Kliman RM: The Hill-Robertson effect: evolutionary consequences of weak selection and linkage.Heredity 2008, 100:19–31.PubMedView Article
Arbiza L, Dopazo J, Dopazo H: Positive selection, relaxation, and acceleration in the evolution of the human and chimp genome.PLoS Comput Biol 2006, 2:e38.PubMedView Article
Bromham L: The genome as a life-history character: why rate of molecular evolution varies between mammal species.Phil Trans R Soc B 2011, 366:2503–2513.PubMedView Article
Thomas JA, Welch JJ, Lanfear R, Bromham L: A generation time effect on the rate of molecular evolution in invertebrates.Mol Biol Evol 2010, 27:1173–1180.PubMedView Article
Wells KD: The Ecology and Behavior of Amphibians. Chicago: The University of Chicago Press; 2007.
Lannoo M (Ed): Amphibian declines. The conservation status of United States species. Berkeley: University of California Press; 2005.
De Magalhães JP, Costa J: A database of vertebrate longevity records and their relation to other life-history traits.J Evol Biol 2009, 22:1770–1774.PubMedView Article
Shao R, Dowton M, Murrell A, Barker SC: Rates of gene rearrangement and nucleotide substitution are correlated in the mitochondrial genomes of insects.Mol Biol Evol 2003, 20:1612–1619.PubMedView Article
Xu W, Jameson D, Tang B, Higgs P: The relationship between the rate of molecular evolution and the rate of genome rearrangement in animal mitochondrial genomes.J Mol Evol 2006, 63:375–392.PubMedView Article
Kurabayashi A, Yoshikawa N, Sato N, Hayashi Y, Oumi S, Fujii T, Sumida M: Complete mitochondrial DNA sequence of the endangered frog Odorrana ishikawae (family Ranidae) and unexpected diversity of mt gene arrangements in ranids.Mol Phylogenet Evol 2010, 56:543–553.PubMedView Article
Ren F, Tanaka H, Yang Z: A likelihood look at the supermatrix-supertree controversy.Gene 2009, 441:119–125.PubMedView Article
Sano N, Kurabayashi A, Fujii T, Yonekawa H, Sumida M: Complete nucleotide sequence of the mitochondrial genome of Schlegel's tree frog Rhacophorus schlegelii (family Rhacophoridae): duplicated control regions and gene rearrangements.Genes Genet Syst 2005, 80:213–224.PubMedView Article
Liu Z-Q, Wang Y-Q, Su B: The mitochondrial genome organization of the rice frog, Fejervarya limnocharis (Amphibia: Anura): a new gene order in the vertebrate mtDNA.Gene 2005, 346:145–151.PubMedView Article
Webster AJ, Payne RJH, Pagel M: Molecular phylogenies link rates of evolution and speciation.Science 2003, 301:478.PubMedView Article
Venditti C, Pagel M: Speciation as an active force in promoting genetic evolution.Trends Ecol Evol 2009, 25:14–20.PubMedView Article
Pagel M, Venditti C, Meade A: Large punctuational contribution of speciation to evolutionary divergence at the molecular level.Science 2006, 314:119–121.PubMedView Article
Davies TJ, Savolainen V, Chase MW, Moat J, Barraclough TG: Environmental energy and evolutionary rates in flowering plants.Proc R Soc B 2004, 271:2195–2200.PubMedView Article
Dubois A: Developmental pathway, speciation and supraspecific taxonomy in amphibians. 1. Why are there so many frog species in Sri Lanka?Alytes 2004, 22:19–37.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.