Efficient Sequencing of Anuran mtDNAs and a Mitogenomic Exploration of the Phylogeny and Evolution of Frogs

Abstract

Anura (frogs and toads) constitute over 88% of living amphibian diversity but many important questions about their phylogeny and evolution remain unresolved. For this study, we developed an efficient method for sequencing anuran mitochondrial DNAs (mtDNAs) by amplifying the mitochondrial genome in 12 overlapping fragments using frog-specific universal primer sets. Based on this method, we generated 47 nearly complete, new anuran mitochondrial genomes and discovered nine novel gene arrangements. By combining the new data and published anuran mitochondrial genomes, we assembled a large mitogenomic data set (11,007 nt) including 90 frog species, representing 39 of 53 recognized anuran families, to investigate their phylogenetic relationships and evolutionary history. The resulting tree strongly supported a paraphyletic arrangement of archaeobatrachian (=nonneobatrachian) frogs, with Leiopelmatoidea branching first, followed by Discoglossoidea, Pipoidea, and Pelobatoidea. Within Neobatrachia, the South African Heleophrynidae is the sister-taxon to all other neobatrachian frogs and the Seychelles-endemic Sooglossidae is recovered as the sister-taxon to Ranoidea. These phylogenetic relationships agree with many nuclear gene studies. The chronogram derived from two Bayesian relaxed clock methods (MultiDivTime and BEAST) suggests that modern frogs (Anura) originated in the early Triassic about 244 Ma and the appearance of Neobatrachia took place in the late Jurassic about 163 Ma. The initial diversifications of two species-rich superfamilies Hyloidea and Ranoidea commenced 110 and 133 Ma, respectively. These times are older than some other estimates by approximately 30–40 My. Compared with nuclear data, mtDNA produces compatible time estimates for deep nodes (>150 Ma), but apparently older estimates for more shallow nodes. Our study shows that, although it evolves relatively rapidly and behaves much as a single locus, mtDNA performs well for both phylogenetic and divergence time inferences and will provide important reference hypotheses for the phylogeny and evolution of frogs.

Introduction

Anura (commonly known as frogs) constitutes the vast majority (6,285 species out of a total of 7,125, or 88%; AmphibiaWeb 2013) of living amphibian diversity, and is also among the most diverse groups of vertebrates. Compared with other amphibians, frogs have a nearly cosmopolitan distribution and occupy a diversity of habitats such as forests, grasslands, rivers, waterfalls, mudflats, mountaintops, and deserts. Because frogs are often used as model organisms to address fundamental issues of morphological, developmental, and biogeographical evolution, a reliable phylogenetic hypothesis and divergence timescale are important. Traditional morphological characters that can be used for anuran systematics are still limited and alone are insufficient to resolve the anuran tree. Ongoing molecular phylogenetic studies now regularly update our understanding of the anuran tree (Roelants and Bossuyt 2005; San Mauro et al. 2005; Frost et al. 2006; Roelants et al. 2007; Pyron and Wiens 2011). However, certain parts of the tree remain unresolved. For example, although the clade Hyloidea is well supported, resolution among its major clades (Hylidae, Bufonidae, Dendrobatidae, etc.) is poorly supported in most analyses.

The mitochondrial genome (mtDNA) rarely undergoes recombination, and so the entire genome is inherited much like a single genetic locus. Although phylogenetic analyses of mtDNA have proven problematic in some cases, such as the position of snakes among squamates (Castoe et al. 2009) or the relationships of rodents among mammals (Reyes et al. 1998), analyses of mtDNA are often congruent with those derived from nuclear genes when appropriate sampling of taxa and phylogenetic analyses are used (Reyes et al. 2004; Kjer and Honeycutt 2007). In addition, not many nuclear genes are widely used in phylogenetic analyses of frogs; for examples, Pyron and Wiens (2011) used nine nuclear genes, but most of the taxa analyzed were incomplete for most of these genes. Thus, the moderate-size mitochondrial genomes are still an attractive data resource for anuran phylogenetics. Indeed, analyses of mtDNA have been shown to be effective for not only amphibian phylogenetics at different taxonomic levels but also for estimation of divergence times back to more than 300 Ma (Mueller et al. 2004; Zhang et al. 2005, 2006, 2008; Zhang and Wake DB 2009, Zhang and Wake MH 2009).

Compared with salamanders and caecilians, mitochondrial genomes of frogs are more difficult to sequence because of methodology usually used to amplify the mitogenome. The usual practice involves amplifying the mitogenome in two or three overlapping long-polymerase chain reaction (PCR) fragments, but frog mtDNA often contains long and GC-rich D-loop regions that are difficult to amplify. These noncoding D-loop regions are not usually useful for higher-level phylogenetic analyses. Therefore, if we omit these troublesome noncoding regions, the remaining portions of the mitogenomes are an accessible data resource for anuran phylogenetic studies.

Here, we describe an efficient sequencing method for anuran mitogenomes by introducing a set of frog-specific primers for conventional PCR amplifications. We then explore the phylogenetic utility of mitogenomes in reconstructing the anuran tree. Our major concerns are 1) whether a reliable anuran tree can be achieved when using rapidly evolving mitogenomes, and 2) how to effectively extract phylogenetic signal from highly heterogeneous mitogenomic data. Moreover, if the mitogenome can provide a robust phylogenetic hypothesis for frogs, its relatively long sequence length should be a good source of information about ages of anuran divergence. To this end, we sequenced 47 near-complete anuran mitogenomes, omitting the D-loop, including 26 anuran families not sampled previously. Prior to this, 43 complete mitochondrial genomes were present in Genbank, representing only 15 of the 53 anuran families. By combining the new sequences and published data, we analyzed the largest set of anuran mitogenomes that has been examined to date (90 frog species plus 4 outgroups). Then, based on the resulting phylogenies, we estimated divergence times for major splitting events during anuran evolutionary history.

Results and Discussion

Efficient Sequencing of Anuran Mitogenomes

Here, we report an efficient and accurate method for sequencing anuran mtDNA by developing frog-specific primer sets. Our previous papers described similar methods for sequencing salamander and caecilian mtDNAs (Zhang et al. 2008; Zhang and Wake DB 2009, Zhang and Wake MH 2009). Evolutionary rate and genome structure are relatively conservative for mtDNAs of salamanders and caecilians; this trait facilitates the identification of conserved regions for designing primers. In contrast, frog mtDNA has considerably faster rates of sequence evolution and more types of gene arrangements, which made it difficult to find conserved regions for primer design. Primers reported here are based on investigations of 43 published anuran mitochondrial genomes and repeated optimizations. Most target fragments are approximately 1,000–1,600 nt, which is suitable for typical Sanger sequencing. However, to span some highly variable regions, we designed primers with somewhat longer intervals using primer walking. Our primers do not cover the mitochondrial D-loop region. We designed some conserved primers flanking this region but very few succeeded. Even using species-specific primers, we were unable to amplify most frog D-loop regions.

We amplified all 47 near-complete anuran mitochondrial genomes in 12 overlapping fragments with a suite of 32 universal primers (fig. 1, tables 1 and 2). Of the 564 target PCR fragments, we successfully amplified 545. Of the remaining 19 fragments, we amplified 7 in 3 longer fragments (F1–F3 in Leptolalax pelodytoides; F6–F7 in Callulina kreffti; F2–F3 in Petropedetes sp.) by different primer combinations. The PCR success rate of our frog-specific primer sets was more than 95% among 47 diverse anuran species. We were able to sequence all target PCR fragments with PCR primers, which greatly reduced the experimental workload. Sequencing a frog mitochondrial genome using this method typically took 3–5 days.

Gene Rearrangements in Anuran Mitogenomes

All but 3 of the 47 new anuran mitogenomes we obtained were complete (table 1) except for the region between the 3′ end of Cytb and the tRNA-Phe gene (mainly the D-loop region). The new sequences were deposited in GenBank under accession JX564852–JX564898.

Generally, we observed two common mitochondrial gene orders in nonneobatrachian and neobatrachian species (fig. 2). In all but one nonneobatrachian mitogenome sequenced here, we found the most frequent gene order to be that previously reported for nonneobatrachian species. In the case of Leptolalax pelodytoides (Megophryidae), we found a unique pattern not previously reported: the tRNA-Val and ND1 genes have become pseudogenes, shortened to 43 bp and 213 nt, respectively; the tRNA-Trp gene is missing from its original position (fig. 2). Because we amplified the corresponding region in a large fragment (from 16S to COII; fig. 1B), this unique pattern is unlikely to be an artifact resulting from sequencing nuclear copies. This is the second report of gene order change for nonneobatrachian frogs; the first is Leiopelma archeyi (Irisarri et al. 2010).

In the neobatrachian species, we found eight different gene orders not reported previously (fig. 2). Most of the novel mitochondrial gene arrangements occur in the WAN(O_L)CY region, consistent with the hypothesis that this region may be a hotspot of gene duplication by virtue of its association with the origin replication of L-strand (O_L) (San Mauro et al. 2006). Rare changes in mitochondrial gene order have attracted great interest because of their potential to provide homoplasy-free evidence of phylogenetic relationships (San Mauro et al. 2006). Given our taxon sample, none of the new gene arrangements diagnoses a clade (figs. 3 and 4). However, many family-level taxa are represented by only one species, so it would be worthwhile to sequence more anuran mitogenomes to explore the phylogenetic utility of these gene arrangements.

Phylogenetic Relationships within Anura

We analyzed two data sets using various partition schemes (table 3): a 94-species dataset and a 91-species data set. The 94-species data set combines two rRNAs, the concatenated tRNAs, and 12 protein-coding genes (all third codon positions are excluded); this alignment contains 11,007 characters. The 91-species data set is a subset of the 94-sp data set, and excludes three species (Hyperolius ocellatus, Arthroleptis poecilonotus, and Phrynobatrachus accraensis) for which we lack portions of the mitogenome. Bayesian analysis of the 94-species data set yielded a reasonably well-resolved phylogeny for frogs: of the 87 internal nodes within anurans, 82 (94%) were recovered with a posterior probability more than 0.95 (fig. 3). Maximum likelihood (ML) analyses on the 94-species data set produced overall similar topologies and comparable branch support to those from the Bayesian analysis (fig. 3). Of the 87 internal nodes recovered by Bayesian analyses, only two were not recovered in the ML analysis and 69 (79%) received strong support (BP > 0.7). The Bayesian and ML analyses of the 91-species data set produced very similar patterns of topology and branch support (fig. 4 and supplementary fig. S2, Supplementary Material online; some differences in fig. 3 are discussed later).

Five major clades of archaeobatrachians (nonneobatrachians) were recovered (figs. 3 and 4): Leiopelmatoidea (Leiopelmatidae + Ascaphidae), Discoglossoidea (Bombinatoridae + Alytidae), Pipoidea (Rhinophrynidae + Pipidae), Pelobatoidea (Scaphiopodidae + Pelodytidae + Pelobatidae + Megophryidae), and Neobatrachia, consistent with most recent studies (Roelants and Bossuyt 2005; Frost et al. 2006; Roelants et al. 2007; Wiens 2007; Pyron and Wiens 2011; Irisarri et al. 2011, 2012). Nonneobatrachian frogs were recovered as successively branching lineages, with Leiopelmatoidea branching first, followed by Discoglossoidea, Pipoidea, and Pelobatoidea (see also Roelants and Bossuyt 2005; Roelants et al. 2007; Pyron and Wiens 2011; Irisarri et al. 2011; but not Frost et al. 2006), different from the branching order of Leiopelmatoidea, Pipoidea, Discoglossoidea, and Pelobatoidea found by some studies (San Mauro et al. 2005; Frost et al. 2006; Wiens 2007). Notably, most previous mitogenomic studies with sparser taxon sampling have placed Discoglossoidea as the closest relatives of Pipoidea (Gissi et al. 2006; but see Irisarri et al. 2011). This grouping was not recovered in our Bayesian and ML analyses (figs. 3 and 4), suggesting that sparse taxon sampling may be responsible for those earlier results (Hillis 1996; Zwickl and Hillis 2002; Heath et al. 2008). Relationships among the four families of Pelobatoidea (figs. 3 and 4) were the same as recovered by the aforementioned studies with the exception of Frost et al. (2006).

Within Neobatrachia, we corroborate the placement of Heleophrynidae as the sister-taxon to all other neobatrachian frogs (Frost et al. 2006; Roelants et al. 2007; Wiens 2007; Pyron and Wiens 2011; Irisarri et al. 2012). This placement received strong support from both Bayesian and ML analyses (posterior probability [PP] = 1.0 and bootstrap proportion [BP] > 80%; figs. 3 and 4). The remaining neobatrachian frogs can be roughly divided into two large clades: one including Myobatrachidae + Hyloidea (=Procoela) and another including Sooglossidae + Ranoidea (=Diplasiocoela; figs. 3 and 4).

Sooglossidae is recovered as the sister-taxon to Ranoidea in our analyses with strong support in the Bayesian analyses (PP = 1.0) but not the ML analyses (BP = 62–69%) (see also San Mauro et al. 2005; Roelants et al. 2007; Wiens 2007; nuclear gene tree of Irisarri et al. 2012; but not Frost et al. 2006 or Biju and Bossuyt 2003). This placement has been widely accepted; that is, Sooglossidae (before Nasikabatrachus was discovered) has been considered the sister-taxon of ranoids, or more informally closer to Ranoidea than to Hyloidea (Lynch 1973). In contrast, based on mitogenomic data but sparse taxon sampling, Irisarri et al. (2012) found Sooglossidae to be the sister-group of all other neobatrachian frogs. Biju and Bossuyt (2003), who additionally included the Indian endemic family Nasikabatrachidae, found the clade Sooglossidae + Nasikabatrachidae (Sooglossoidea; Nasikabatrachidae not sampled by us) to be the sister-group of all other neobatrachians. Pyron and Wiens (2011) found a weakly supported relationship between Sooglossoidea and all neobatrachian frogs to the exclusion of Heleophrynidae (BP = 56%), whereas Frost et al. (2006) placed Sooglossoidea as the sister-group of Myobatrachidae + Calyptocephalellidae + Hyloidea. In summary, almost all possible arrangements have been found, and currently, neither mtDNA nor the available nuclear data resolve the placement of Sooglossoidea with decisive support although the arrangement ([Sooglossidae + Nasikabatrachidae], Ranoidea) is shared between our results and those of Roelants et al. (2007). Additional nuclear gene sequences of Nasikabatrachidae and other species of Sooglossidae may provide confirmatory evidence given the long branches that separate this taxon from Sooglossidae.

Within Hyloidea, 23 of 24 internal nodes (96%) are well supported (PP > 0.95; fig. 3); however, in half of the nodes bootstrap support is less than 80. Similarly, Pyron and Wiens (2011) found that only 11/27 of the nodes relating family group taxa had bootstrap values more than 50%. All hyloid families represented by two or more species were strongly supported except for Hylidae: Dendrobatidae, Hemiphractidae, Bufonidae, Centrolenidae, and Leptodactylidae.

In all analyses (figs. 3 and 4), we find strong support for placing a well-supported clade Eleutherodactylidae + Craugastoridae + Strabomantidae as the sister-group to all other hyloid frogs (see also Darst and Cannatella 2004; Frost et al. 2006; Pyron and Wiens 2011). However, two multilocus studies, Roelants et al. (2007; 4 nuclear genes and 1 mtDNA fragment) and Heinicke et al. (2009; 11 nuclear genes and 3 mtDNA fragments), placed Rhinodermatidae in the same position, but without strong support. Rhinodermatidae is the sister-group of Hylidae in our tree, but with only moderate support in Bayesian analyses (figs. 3 and 4), leaving its placement uncertain.

Our data recovered a well-supported clade including Odontophrynidae, Ceratophryidae, Hemiphractidae (represented by Gastrotheca and Cryptobatrachus, discussed later), Alsodidae, and Telmatobiidae (BP = 76%; PP = 1.0). The relatively close relationships among Odontophrynidae, Ceratophryidae, Alsodidae, and Telmatobiidae were also reported by Darst and Cannatella (2004) and Pyron and Wiens (2011), but these had weak support. The placement of Hemiphractidae (back-brooding treefrogs), broadly construed to include Gastrotheca, Stefania, Cryptobatrachus, Flectonotus, and Hemiphractus, has varied wildly, depending on the study. Frost et al. (2006) found the group to be nonmonophyletic, consisting of three independent lineages, but Wiens (2007) and Pyron and Wiens (2011) found it to be monophyletic. In contrast to our findings, most studies, including Heinicke et al. (2009) and Darst and Cannatella (2004), found at least some component of Hemiphractidae to be a deeply branching lineage among hyloids.

In summary, relationships among most of the family-level clades in Hyloidea have generally defied a consensus; nodes that are well supported in one study are often poorly supported in other studies. One general result is that the strongly supported clade Craugastoridae + Eleutherodactylidae + Strabomantidae (this node is a subset of the clade Terrarana [Hedges et al. 2008]; hence, we do not use that clade name) diverges very early within the Hyloidea. Because our taxon sampling for Hyloidea is limited and many families and subfamilies lack representative species, sampling more taxa for mitogenomic data will almost certainly alter our results, but it is not clear that it will result in a well-supported phylogenetic reconstruction for these groups.

Within Ranoidea (fig. 3), our phylogenetic estimate is generally similar to others (Darst and Cannatella 2004; Frost et al. 2006; Roelants et al. 2007; Pyron and Wiens 2011). A clade containing Brevicipitidae, Hemisotidae, Hyperoliidae, and Arthroleptidae is repeatedly recovered. van der Meijden et al. (2005) applied the epifamily name Arthroleptoidae (Dubois 1992) to this clade, although the name with priority is Brevicipitoidae Bonaparte, 1850. Frost et al. (2006) called this clade Afrobatrachia. When using all data (the 94-species data set), the Brevicipitoidae was weakly supported by the Bayesian analysis and not recovered by the ML analysis (fig. 3). When three species with missing data (Hyperolius ocellatus, Arthroleptis poecilonotus, and Phrynobatrachus accraensis) are removed (the 91-species data set), both Bayesian and ML analyses supported the monophyly of Brevicipitoidae (fig. 4). That missing data are causing the problem is clear because the topologies of the 91-species and 94-species analyses are otherwise identical.

The same missing data also influenced the placement of Microhylidae. In our analyses of the 94-species data set, which has missing data, Microhylidae is the sister-group of Ranoidae (Natatanura of Frost et al. [2006]) (fig. 3), an arrangement generally not supported by nuclear genes (but see van der Meijden et al. 2004). In our analyses of the 91-species (complete) data set, we recovered a sister-group relationship between Microhylidae and Brevicipitoidae in both Bayesian and ML analyses (fig. 4). This latter topology is reported by most studies using multiple loci (Frost et al. 2006; Roelants et al. 2007; Pyron and Wiens 2011) as well as 12S and 16S genes only (Darst and Cannatella 2004). Thus, contrary to several papers claiming that missing data are not problematic (many cited in Wiens and Morrill 2011), our analysis shows definitive effects of missing data.

The nine microhylid species (six subfamilies represented) formed a highly supported clade (figs. 3 and 4). Within this clade, relationships were not well resolved; similar lack of resolution has been reported in most other studies (van der Meijden et al. 2007; Roelants et al. 2007; Pyron and Wiens 2011), suggesting the possibility of rapid lineage divergence early in the evolution of this group.

The ML and Bayesian analyses placed Cophixalus sp. (Asterophryinae) as the sister-group to the remaining microhylids. However, given that this species has a longer branch than other microhylids and the internal branches are relatively short in this clade (figs. 3 and 4), this placement may have resulted from long-branch attraction (Felsenstein 1978). To explore this, we constructed a small data set including only the microhylids using Brevicipitoidae as the outgroup, and performed a CAT-model analysis, which is purported to be less susceptible to long-branch attraction artifacts. We found (fig. 5) that a clade containing Phrynomerinae and Gastrophryninae is the sister-group of other microhylids, whereas the rapidly evolving Cophixalus sp. (Asterophryinae) is the sister-group of Dyscophinae and Microhylinae. This result is very similar to that of van der Meijden et al. (2007), which was based on nuclear genes. However, the internal branches are very short in both analyses.

As in most previous studies (Frost et al. 2006; Roelants et al. 2007; Wiens 2007; Pyron and Wiens 2011), we recovered strong support for the clade Ranoidae (van der Meijden et al. 2005), also called Ranidae [sensu lato] in the traditional sense of Bossuyt et al. [2006], figs. 3 and 4). However, the family-level relationships within this clade are not well resolved, especially for the placement of the African families Phrynobatrachidae, Pyxicephalidae, Petropedetidae, and Ptychadenidae. We also found relatively strong support for Mantellidae and Rhacophoridae as sister-taxa and together as the sister-taxon of Ranidae (also reported by Pyron and Wiens 2011), and for the sister-group relationship between Pyxicephalidae and Petropedetidae (also reported by Frost et al. 2006; Roelants et al. 2007; Pyron and Wiens 2011). Given the very short internal branches and long terminal branches in this clade, resolving relationships within this group may require much denser taxon sampling (Hillis 1996).

Reliability of Single Mitochondrial Genes for Anuran Phylogenetics

Although sequencing a near-complete frog mitochondrial genome is neither difficult nor expensive using our method, many researchers typically sequence preferred mitochondrial genes and combine them with nuclear gene data. Which mitochondrial gene should be chosen first? To answer this, we need to better understand the contribution of each mitochondrial gene to the mitogenome phylogeny. We calculated Robinson–Foulds distance and K-scale factor (Soria-Carrasco et al. 2007) for each mitochondrial gene. In brief, Robinson–Foulds distance is the topological difference between the individual gene trees and the mitogenomic tree, our reference tree calculated from the concatenation of all genes. K-scale factor is the ratio between the global divergence (similar to average branch length) of the individual gene tree and the mitogenomic tree. We found important variability in the phylogenetic signal of the different genes (table 4). Six of the 14 genes were especially informative in topology resolution (lower Robinson–Foulds distance), in the order of 16S, ND2, ND5, ND4, 12S, and ND1. The widely used Cytb gene provides relatively little topological resolution, ranking third from the bottom (table 4). The resolution provided by COI, recommended by some workers for barcoding (e.g., Hawkins et al. 2007), was moderate (table 4). Therefore, 16S, recommended by Vences et al. (2005), could profitably be used together with COI for anuran DNA barcoding in which it is desirable to identify an unknown against a clade spanning 200 My. According to K-scale factors, trees based on 16S and 12S sequences were closest to the reference tree, followed by trees based on ND3, ND1, ND2, and ND5. The reasons for the different results indicated by the RF distances and K-scale factor measures with respect to the ND3 gene are unclear. Perhaps the short length of the ND3 fragment (226 nt) accounts for this difference.

Taking both criteria (RF distance and K-scale factors) into account by summing the ranks of genes in both categories, the five genes that provide the best phylogenetic resolution are 16S, ND2, 12S, ND5, ND1, and COI. In most previous molecular studies, 12S and 16S are often chosen as the two representative genes for mitochondrial genomes. However, our results suggest that if only two mitochondrial genes are to be selected for sequencing, 16S should be the first choice, but ND2 would be preferable to 12S. We note that the reference tree has divergences at all depths but only 16/88 divergences are less than 50 Ma. Cytb and COI remain reasonable choices for analysis of closely related species.

Origin and Evolution of Frogs

We used our mitogenome data set to estimate the divergence times for major lineages of anurans, using the methods implemented in MultiDivTime (Thorne and Kishino 2002) and BEAST (Drummond and Rambaut 2007). Overall, the two dating methods yielded similar results (supplementary table S1, Supplementary Material online); we specified the tree topology in both analyses. Some groups of nearby nodes showed a bias in which one method or the other gave an older estimate. It is not clear whether this effect might be related to the choice of calibration points (these were the same in both analyses) or of the computational methods in MultiDivTime and BEAST. Because we lacked evidence that one method provided consistently better estimates compared with the other, we averaged the time estimates of the two dating methods. We use these averages in our discussion and provide the divergence time for each analysis (supplementary table S1, Supplementary Material online).

Many studies have addressed the question of the origin of living frogs and dated the basal split during the Late Permian (∼265 Ma; San Mauro et al. 2005), the Early Triassic (225–250 Ma; Roelants and Bossuyt 2005; Roelants et al. 2007; Pyron 2010, 2011), or the Late Triassic (206 Ma; San Mauro 2010). Our chronogram (fig. 6; supplementary table S1, Supplementary Material online) places the early diversification of the modern frogs (Anura) in the Early Triassic (∼243.7 Ma; node 1, supplementary table S1, Supplementary Material online). This result is in line with previous molecular clock studies, in that the origin of crown-group frogs predated the Jurassic. Two well-characterized fossil species, Prosalirus bitis (Early Jurassic; Shubin and Jenkins 1995) and Vieraella herbstii (Early Jurassic; Reig 1961), belong to the stem-group Salientia rather than to the crown-group Anura. Citing these fossils, Marjanović and Laurin (2007) suggested that the crown-group frogs (Anura) originated much later, no earlier than 185 Ma, although they did not explain how they reached this opinion. The apparent conflict between fossil and molecule-based analyses might be due to the incomplete fossil record of crown-group frogs. Paleontological studies rely on fossils to directly infer minimal divergence times and molecule-based studies rely on fossils to calibrate chronograms. New fossil records, especially from basal crown-group frogs, will likely impact both types of studies. Until key fossil taxa are found, debates about the origin of living frogs will continue.

The radiation of neobatrachians (advanced frogs) is a major event during the anuran evolutionary history; this group comprises over 96% of all living anurans. Dating this event requires inclusion of five primary neobatrachian lineages: Heleophrynidae, Sooglossoidea, Myobatrachidae + Calyptocephalellidae, Ranoidea, and Hyloidea, but only a few studies have included representatives of all these groups. Based on the nuclear Rag1 gene, San Mauro et al. (2005) provided the first molecular estimate for the node Neobatrachia as 162 Ma. Using a five-gene data set of 120 anuran taxa, Roelants et al. (2007) found a similar estimate of 167 Ma. Using Rag1 sequences and treating fossils as terminal taxa, Pyron (2011) produced a younger time estimate of 125 Ma. Our estimate for the node Neobatrachia (163.1 Ma; MultiDivTime, 132.0, and BEAST, 194.2; node 5, supplementary table S1, Supplementary Material online) is close to the two older estimates. Currently, all molecular evidence indicates that the appearance of the neobatrachian clade took place either in the Late Jurassic or the Early Cretaceous.

The Hyloidea comprises approximately half of the extant frog species and previous molecular dating studies suggest that the initial diversification of hyloids commenced relatively late, near the K-T boundary ∼64 Ma (San Mauro et al. 2005; Roelants et al. 2007). Using nine genes, Heinicke et al. (2009) gave an older date (∼74 Ma) for this event. Our estimate for the origin of Hyloidea (node 28, supplementary table S1, Supplementary Material online) is 110.5 Ma (MultiDivTime, 96.1, and BEAST 124.9), older than previous estimates.

The frog fossil Beelzebufo ampinga was discovered from the Late Cretaceous (Maastrichtian, 70–65 Ma) of Madagascar (Evans et al. 2008; Ruane et al. 2011) and assigned to Hyloidea. Assuming that this fossil is within Hyloidea and not its sister-group, the appearance of hyloids, based on calibration using Beelzebufo, would have predated the K-T boundary (∼65 Ma), an estimate that is older than most based on molecular data. Hyloids occur primarily in South America and there are no living hyloids in Madagascar. Therefore, the origin of living Hyloidea may predate the time when Madagascar lost contact with Gondwana, or else Hyloidea went extinct on Madagascar after the island split from the rest of Gondwana. Although there is debate about the timing of splits in the Gondwanan landmass (reviewed by Krause et al. 2006), Madagascar has been isolated from other Gondwanan landmasses for approximately 120 Ma. Our estimates of the origins of Hyloidea (node 28, MultiDivTime, 96.1, and BEAST 124.9) are in line with this paleobiogeographic inference.

In a re-analysis, the data of Roelants et al. (2007) and Ruane et al. (2011) used 65–70 My to calibrate the node Ceratophryinae in concert with three other fossil calibrations. It yielded divergence dates for crown Ceratophryinae of 64.0–67.3 Ma, whereas its exclusion yielded divergence dates for Ceratophryinae of 11.6–13.6 Ma, a 4-fold, highly significant difference. The most plausible explanation is that the phylogenetic placement of Beelzebufo is incorrect. However, there seems to be little doubt that Beelzebufo is at least a stem-neobatrachian based on morphological character analysis (Evans et al. 2008).

Other Cretaceous putative neobatrachians (Báez et al. 2009, 2012) span the Albian-Aptian (108–112 Ma) to Maastrichtian (65–70 My). Our conservative use of the Maastrichtian period as the calibration for the divergence between Neobatrachia and Pelobatoidea (node 4) apparently did not bias the estimate toward a younger age because the estimated divergence of node 4 is the Upper Triassic. Although these several Cretaceous fossils are putatively “neobatrachian,” implying placement within crown Neobatrachia, in fact they cannot be excluded as stem neobatrachians that diverged along the branch between nodes 4 and 5 (fig. 6).

The current distribution of Ranoidea shows a prevalence of continental-scale endemism, suggesting that continental breakup (especially Gondwana) played a key role in the diversification of these frogs. Previous studies suggest the origin of Ranoidea took place 99 Ma (San Mauro et al. 2005), 119 Ma (Roelants et al. 2007), or 133 Ma (Van Bocxlaer et al. 2006). All these time estimates largely coincide with the breakup of Gondwana. Our time estimate for this event is 133.1 Ma (MultiDivTime, 108.0 and BEAST, 158.2, node 8, supplementary table S1, Supplementary Material online).

Previous divergence time estimates for frogs were mainly based on nuclear gene data. In general, mitochondrial genes evolve much faster than most nuclear genes, and it has been argued that ancient divergence times cannot be correctly estimated by these fast-evolving genes because of mutational saturation (Bossuyt and Roelants 2009; Zheng et al. 2011). In a study of ray-finned fishes, Near et al. (2012) found a tendency for mitochondrial data to yield older divergence estimates than nuclear data. To check the overall deviation between times from our mitogenomic and previous nuclear studies, we plotted the divergence times for 27 nodes in this study against the ages of corresponding nodes of the phylogeny from Roelants et al. (2007), which are inferred from multiple genes (fig. 7). Our analysis yielded estimates that were close to those of Roelants et al. (2007) for deep nodes (>150 Ma) but were older for shallow nodes (<150 Ma, fig. 7). The congruence between the two studies in the deep nodes may be due to the use of similar calibrations at those deeper nodes. The main difference in time estimates between the mtDNA and nuclear data lies within Neobatrachia. Interestingly, Irisarri et al. (2012) found that within Neobatrachia the mtDNA substitution rate is faster than that for nuclear genes. Other studies (Bossuyt and Roelants 2009; Zheng et al. 2011; Near et al. 2012) have suggested that mtDNA data yield older time estimates because of uncorrectable saturation. Therefore, we infer that within Neobatrachia the older time estimates from mtDNA are perhaps due to the saturation of mtDNA, and also partly derived from different substitution rates between mitochondrial and nuclear genes. Clearly, additional work focusing on heterogeneous rates of nuclear and mitochondrial genes is needed (see Irisarri et al. 2012).

We noted earlier that relationships among the families of hyloids, families of Ranoidae, and subfamilies of Microhylidae were poorly to moderately supported. These three clusters of nodes (fig. 6) occur between 53.3 (node 37; MultiDivTime, 56.1; BEAST, 50.5) and 116.3 (node 50; MultiDivTime, 90.8; BEAST 141.7), spanning the Early Cretaceous to the early Paleogene. A similar pattern was noted by Roelants et al. (2007), who argued that the diversifications of Ranoidae and Microhylidae occurred before the K-T boundary (∼65 Ma), but that the diversification of Hyloidea followed the K-T boundary (see also Heinicke et al. 2009). In our analyses, almost all nodes within Hyloidea diverged before the K-T boundary (fig. 7). The 95% credible intervals of divergence times of each of the three clusters overlap considerably, and thus we argue for one worldwide diversification event of three primarily Gondwanan lineages. Most of the neobatrachian lineages ranked as families or subfamilies diverged during a Cretaceous “plague of frogs,” and a modern herpetologist who was transported in time to the Cretaceous would likely have little trouble identifying most frogs that s/he encountered to the family level.

Taxonomy

Our phylogenetic analyses of mitogenomes generally confirm the results of analyses of multilocus data sets using fewer nucleotides. Given that our analyses use only the mitogenome, we make no new taxonomic proposals, but some comments about how our results relate to previously proposed taxonomic arrangements of anurans, and particularly to those proposed by Frost et al. (2006), who provided new names for almost all nodes above the family-group level.

In cases in which Frost et al. (2006) applied names to previously undiscovered taxa, and our results support the monophyly of those groups, we have used those names (e.g., Phthanobatrachia). In cases in which a clade is associated with an earlier well-known name, we continue to use the earlier name. Two salient examples are Hyloidea and Ranoidea. Frost et al. (2006) proposed Hyloides and Ranoides to replace these. They used Ranoidea for a much more restricted node (only 2 families, rather than 14 as in AmphibiaWeb [2013]) and abandoned the use of Hyloidea.

We continue the use of Hyloidea because Darst and Cannatella defined the name in a phylogenetic context (sensu de Queiroz and Gauthier 1992) and it has been employed in that sense by several workers: for example, Lynch (1973), Cannatella and Hillis (2004), van der Meijden et al. (2005), Wiens (2007), Roelants et al. (2007), and Pyron and Wiens (2011). Similarly, Ranoidea was defined phylogenetically by Ford and Cannatella (1993), although the content was emended slightly following placement of Dendrobatidae as within Hyloidea (Hedges and Maxson 1993; Darst and Cannatella 2004). Ranoidea was also used with the same content by Lynch (1973), Cannatella and Hillis 2004, van der Meijden et al. (2005), Wiens (2007), and Pyron and Wiens (2011).

Similarly, we prefer existing and well-established names for other taxa: for examples, Pelobatoidea rather than Anomocoela; Leiopelmatoidea rather than Amphicoela; Pipoidea rather than Xenoanura; Procoela rather than Nobleobatrachia; and Diplasiocoela for Sooglossidae + Ranoidea (a clade not recovered by Frost et al. [2006]). A similar rationale applies to our use of Bombinanura, Pipanura, and Acosmanura (fig. 6). Thus, no part of our taxonomy is unprecedented. We also note that the International Code of Zoological Nomenclature does not apply to ranks above the family-group level (including superfamily, epifamily, etc.).

Summary

We describe a new method for rapid and efficient sequencing of mitogenomes. We used a large sample of new mitogenomes (90 species, of which 47 were new) to infer the phylogeny of anurans. The sample includes 39/53 families across a wide taxonomic coverage of anurans. Several gene order rearrangements were discovered, but when combined with published gene orders, most rearrangements are autapomorphies, suggesting that sampling at a finer scale is needed. Our phylogeny confirms strongly supported relationships reported by several studies using multiple nuclear genes. Similarly, the regions of the tree that were poorly resolved (e.g., within Hyloidea) are typically resolved by these other studies. This suggests that, despite criticisms of the use of mitogenomes in phylogenetic inference, these remain an important and accessible source of data (e.g., Powell et al. 2013).

Estimates of divergence times of deeper node in the mitogenome tree are generally similar to those of multilocus studies, whereas shallower nodes (Neobatrachia and younger) have older divergence times. However, the 95% credible intervals of these estimates overlap broadly in many cases. Comparison of our chronogram with that of Roelants et al. (2007) (the largest such study) indicates that the timing of poorly resolved rapid divergences within three clades (Hyloidea, Ranoidae, and Microhylidae) cannot be distinguished from each other. This finding calls into question the proposed divergence of Hyloidea after the K-T boundary.

Materials and Methods

Taxon Sampling

According to AmphibiaWeb (March 2013), Anura contains 53 recognized families. Before this study, complete mitochondrial genomes for 43 species were available in GenBank, representing only 15 families: Ascaphidae, Leiopelmatidae, Bombinatoridae, Alytidae, Rhinophrynidae, Pipidae, Pelobatidae, Hylidae, Bufonidae, Microhylidae, Ranidae, Dicroglossidae, Rhacophoridae, and Mantellidae. We sampled 47 frog species covering 39 anuran families, of which most families were not previously sampled. Our new data plus published frog mitogenomes total 90 frog species. At the time that our manuscript was submitted, Irisarri et al. (2012) reported six new anuran mitogenomes from five anuran families. All but one of these five families (Calyptocephalellidae) was included in our sample. Complete mitogenomes of two salamanders and two caecilians were retrieved from GenBank to serve as outgroup taxa. The details for newly determined sequences and published mitogenomes are given in table 1 and supplementary table S2, Supplementary Material online, respectively.

MtDNA Amplification and Sequencing

The conventional strategy of sequencing an amphibian mitogenome is to amplify the genome in two or three large fragments by long PCR (Zhang et al. 2005). However, we were unable to obtain long PCR products for many frog species despite many optimization efforts, especially for fragments spanning the D-loop region. Therefore, we amplified the frog mitogenomes in contiguous and overlapping short fragments. We aligned available mitogenomes and designed a suite of 32 primers (table 2). A frog mitogenome (D-loop region not included) can be amplified in 12 fragments using these primers (fig. 1A, supplementary fig. S3, Supplementary Material online).

Total DNA was purified from frozen or ethanol-preserved tissues (liver or muscle) using the DNeasy Blood and Tissue Kit (Qiagen). PCRs were performed with Ex Taq DNA Polymerase (Takara) in total volumes of 25 μl, using the following cycling conditions: an initial denaturing step at 96 °C for 2 min; 35 cycles of denaturing at 94 °C for 15 s, annealing at 45–55 °C for 60 s, and extension at 72 °C for 2 min; and a final extension step of 72 °C for 10 min. PCR products were purified either directly via ExoSAP (USB) treatment or gel-cutting (1% TAE agarose) using the Gel Purification Kit (Qiagen). Sequencing was performed directly with the corresponding PCR primers using the BigDye Deoxy Terminator cycle-sequencing kit v3.1 (Applied Biosystems) on an automated DNA sequencer (ABI PRISM 3730) following the manufacturer's instructions. For some large PCR fragments, specific primers for primer walking were designed using newly obtained sequences. To ensure we did not amplify nuclear copies of mitochondrial fragments, we carefully examined the contig assemblies and found no incongruence in any overlapping regions and no stop codons in protein-coding genes.

Sequence Alignments and Data Partitions

The L-strand encoded NADH6 gene and the D-loop were not used because of their heterogeneous base composition and poor phylogenetic performance. Twelve protein-coding, 22 tRNA, and two rRNA gene sequences were separately aligned using the program PRANK (Löytynoja and Goldman 2008) at default settings for nucleotides (for RNA genes) or translated amino acids (for protein genes). PRANK is a probabilistic multiple alignment program that generates posterior probabilities for each aligned site. Therefore, for alignment refinements, sites with posterior probability less than 0.95 were excluded. A few tRNA and protein genes are incomplete or missing in some mitogenomes. These sites were treated as missing data. All 22 tRNA alignments were then combined into a concatenated alignment. Finally, all 15 refined alignments (2 rRNAs, 1 concatenated tRNA, 12 protein-coding genes) were concatenated for further analyses. Two data sets were analyzed: the first one included all 94 species and the second one excluded three species: Hyperolius ocellatus, Phrynobatrachus accraensis, Arthroleptis poecilonotus, because their mitogenomes are incomplete, hereafter “91-species” (although the mitogenome of the second Phrynobatrachus species, P. keniensis, is complete). To investigate the extent of substitution saturation, we plotted uncorrected pairwise distance against corrected pairwise distance (GTR + Γ + I model) for six major partitions of our data: 12S, 16S, tRNAs, all first codon positions, all second codon positions, and all third codon positions. We found slight saturation for the first five partitions but very strong saturation for the third codon positions (supplementary fig. S1, Supplementary Material online). Therefore, all third codon positions of protein-coding genes were excluded from further phylogenetic analyses.

To improve the fit of the substitution model to the data, we compared different data partitioning schemes according to the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) using the program PartitionFinder v1.0 (Lanfear et al. 2012). We divided the data sets into 1, 5, 15, and 27 partitions, referred to as 1p, 5p, 15p, and 27p. Except for the unpartitioned scheme (1p), each scheme included a separate partition for 12S, 16S, and the concatenated tRNAs. Strategies differed in the partitioning of protein-coding genes: 5p used two separate partitions for all first and second codon positions; 15p defined 12 separate partitions for each of the 12 protein-coding genes; 27p defined 24 separate partitions for each first and second codon position in each gene. Both AIC and BIC identified the 27p as the best partitioning scheme for the 94-species and 91-species alignments (table 3).

Phylogenetic Analyses

The two data sets (94-species and 91-species) were subjected to ML and Bayesian inference (BI) analyses under various partitioning strategies. Partitioned ML analyses were performed using RAxML v7.2.6 (Stamatakis 2006) with the “-q” option. Separate GTR + Γ + I substitution models were defined for each partition. An analysis combining 100 replicate searches was applied to find the optimal tree (-f d option), and branch support for each node was evaluated with 500 rapid bootstrap replicates (-f a option). Partitioned BI analyses were conducted in MrBayes 3.2 (Ronquist et al. 2012). The best-fitting model for each partition was selected by using BIC implemented in PartitionFinder v1.0 (Lanfear et al. 2012). Parameters of the substitution models and among-site rate variation were unlinked across partitions, and partition-specific rate-multipliers were used to account for variation in evolutionary rates across partitions. Three independent Markov chain Monte Carlo (MCMC) runs were performed with one cold and three heated chains (temperature set to 0.1) for 20 million generations and sampled every 1,000 generations, yielding a sample of 20,000 generations. We judged the chain to be converged when the effective sample size (ESS) values of all parameters (estimated in Tracer version 1.5; Drummond and Rambaut 2007) all exceeded 100 after the first 25–50% generations were discarded. Topologies and posterior probabilities were estimated from the remaining generations of three runs and compared for congruence.

The CAT model (Lartillot and Philippe 2004) implemented in PhyloBayes has been shown in some contexts to be less prone to systematic errors (e.g., long-branch attraction artifacts) than other models (Lartillot et al. 2007). In cases where long-branch attraction was suspected of having influenced phylogenetic inference, the CAT model Bayesian analyses were performed using PhyloBayes v3.2f (Lartillot et al. 2009) with Poisson exchange rates. The distribution of rates across sites is emulated by a discrete gamma distribution with four categories. Two separate runs were performed for 10,000 cycles. Chain stationarity between two runs was evaluated by using a threshold of 0.1 for the “maxdiff” parameter (calculated by the “bpcomp” program) as suggested by the authors, after discarding the first 50% of cycles as burn-in.

We characterized node support (either bootstrap proportions or posterior probabilities) in all analyses as strong (≥95.0), moderate (70.0–94.9), or weak (60.0–69.9).

Single Gene Reliability

To explore the contribution of each individual gene to the concatenated data set, 14 single-gene phylogenies (2 rRNA and 12 protein-coding genes) from the 91-species data set were reconstructed with GTR + Γ + I model using MrBayes (10,000,000 generations). We scored each single-gene tree based on Robinson–Foulds distances and the K-scale factor as calculated by the program Ktreedist (Soria-Carrasco et al. 2007) using the 91-species Bayesian topology as the reference tree.

Estimating Divergence Times

Two relaxed clock methods, MultiDivTime (Thorne and Kishino 2002) and BEAST 1.7.0 (Drummond and Rambaut 2007), were used for molecular dating analyses based on the 91-species data set. The data set was divided by genes into 15 partitions. We used Caudata as the first outgroup and Gymnophiona as the second outgroup. Some studies (e.g., Fong et al. 2012) have found limited support for ([Gymnophiona + Caudata], Anura) rather than ([Gymnophiona], [Caudata + Anura], suggesting that our outgroup should be [Gymnophiona + Caudata]). Even if this alternative tree were true, our results are not affected because we rooted the tree a priori between salamanders and frogs and calibrated that divergence using Triadobatrachus; that is, we did not use the analysis do determine the root of Amphibia. Given our a priori root, the position of caecilians as either the sister-group of Anura + Caudata or as the sister-group of Caudata only would have no effect.

The prior for the root age is 275 ± 50 Ma, which is sufficiently broad to include both recent molecular and fossil estimates (Marjanović and Laurin 2007; Anderson et al. 2008; Hugall et al. 2007; Zhang and Wake DB 2009; Zhang and Wake MH 2009; Pyron 2011). We imposed five time constraints to calibrate the molecular clock: 1) the salamander–frog split was constrained to be at least 250 Ma based on the fossil Triadobatrachus massinoti (Rage and Rocek 1989); 2) the split between Cryptobranchidae and Hynobiidae (salamanders) was constrained to be at least 145 Ma based on the fossil Chunerpeton tianyiense (Gao and Shubin 2003), which is a more conservative minimum age for this problematic fossil than the original assumption of a Middle Jurassic age by Gao and Shubin (2003); 3) a minimum age of 164 Ma for the Discoglossoidea-Pipanura split, based on the fossil Eodiscoglossus oxoniensis, of Bathonian age (Evans et al. 1990); 4) a minimum age of 151 Ma for the Rhinophrynidae-Pipidae split, based on the fossil Rhadinosteus parvus, of Kimmeridgian age (Henrici 1998); 5) a minimum age of 65 Ma for the Pelobatoidea-Neobatrachia split, based on the fossil Beelzebufo ampinga from the Late Cretaceous (Evans et al. 2008).

The calibration for the Neobatrachia-Pelobatoidea split (fig. 6, node 4) requires explanation. The Cretaceous Malagasy neobatrachian frog Beelzebufo was dated at 65–70 Ma by Evans et al. (2008). Their phylogeny placed Beelzebufo within the clade Ceratophryinae, as the sister-group of Ceratophrys. This indicates the presence of ceratophryines in Madagascar at 65–70 Ma. Given that some estimated ages of crown Hyloidea are 55–65 Ma (e.g., San Mauro et al. 2005), Evans et al. (2008) noted that it was difficult to reconcile the age of Beelzebufo with its putative inclusion within Ceratophryinae.

Analyses of several Cretaceous putative neobatrachians (Báez et al. 2012) indicate a diverse Gondwanan assemblage spanning the Albian-Aptian (108–112) to Maastrichtian (65–70 Ma). We chose the uppermost Cretaceous (65–70) as the calibration for the divergence between Neobatrachia and Pelobatoidea because the putative neobatrachians may actually be stem-neobatrachians rather than crown-neobatrachians.

For the MultiDivTime analysis, the priors for the mean and standard deviation (SD) of the ingroup root age were set as rttm = 2.75 and rttmsd = 0.25, respectively. The prior mean and SD for the gamma distribution describing the rate at the root node (rtrate and rtratesd) were both set to 0.2422. These values were based on the median of the substitution path lengths between the root and each terminal, divided by rttm (as suggested by the authors). The autocorrelation parameter prior (brownmean) and its SD (brownsd) were set to 0.727, such that brownmean multiplied by the rttm prior (2.75) equals 2.0. After a burn-in period of 250,000 generations, MCMC chains were sampled every 100 generations until 10,000 samples were acquired. Two separate runs were performed and similar results were observed.

In the BEAST analysis, the uncorrelated lognormal model was used to describe the relaxed clock and the GTR + Γ + I was used to describe the substitution model for the 15 partitions of the data set. The Yule process was used as the speciation model. A normal distribution prior (mean = 275, SD = 25) was used for the root. The five time constraints were used as the lower bounds with uniform distributions (the upper bounds of the uniform distributions were set to 10,000 Ma). Two separate runs were performed for 20 million generations, and sampled every 1,000 generations. Results of the two runs were compared for convergence after a burn-in of approximately 25%.

Acknowledgments

The authors thank Xing Xing Shen for general laboratory support. Two anonymous reviewers provided insightful comments on an early version of this manuscript. For loans of tissue samples, they thank Travis LaDuc, Texas Natural Science Center, University of Texas, Austin; Ronald Nussbaum, University of Michigan; and Jim McGuire and Carol Spencer, Museum of Vertebrate Zoology, University of California, Berkeley. This work was supported by National Natural Science Foundation of China grant 31172075 and 30900136 to P.Z. and by the National Science Foundation grant EF-0334952 to D.C.C. and D.M.H. and EF-0334939 to D.B.W.

References

Author notes