Phylogenomics reveals rapid, simultaneous diversification of three major clades of Gondwanan frogs at the Cretaceous–Paleogene boundary

Data Characteristics.

We assembled a de novo 164-species dataset by using 95 nuclear genes (Table S1) and 88,302 nt from 156 frog species and 8 outgroups; this matrix is 89.6% complete. To increase coverage of anuran families, we added sequences of RAG1 and CXCR4 from GenBank of 145 additional anuran species. This 309-species dataset contains 88,386 nt and is 48.2% complete. The 164-species and 309-species matrices are available from the Dryad Digital Repository (doi:10.5061/dryad.12546).

Higher-Level Phylogenetic Relationships of Frogs.

Maximum-likelihood (ML) and Bayesian analyses of concatenated genes in the 164-species dataset produced identical trees except for two nodes with low support (Figs. S1 and S2). The ASTRAL species tree differed from the ML tree at eight poorly supported nodes (Fig. S3). Overall support is high: 94% of 155 nodes within frogs have a bootstrap value (BS) ≥70% (Fig. S1), and 97% have Bayesian posterior probabilities (BPPs) ≥0.95 (Fig. S2). The ASTRAL species-tree method produced BSs ≥70% at 84% of the nodes (Fig. S3). Although ML analysis of the 309-species dataset has weaker support (78% of nodes have BSs ≥70%; Fig. S4), the basic topology of the ML tree is similar. Therefore, we used the ML tree as our primary hypothesis (Fig. 1A) for estimating the chronogram.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Time-calibrated phylogenetic tree of frogs and the pattern of net diversification rate across time. (A) Evolutionary chronogram based on 95% nuclear genes and 20 fossil age constraints. Gray bars represent the 95% credibility interval of divergence time estimates. Divergence time estimates and corresponding 95% credibility intervals for all nodes are provided in Table S2. Note that the initial diversification of the three major frog clades: Hyloidea (blue), Microhylidae (purple), and Natatanura (green) took place simultaneously near the KPB (dashed red line). (B) Rate-through-time plot of extant frogs indicates an increase in diversification rate at the end of the Cretaceous.

Archaeobatrachian relationships are identical to those of most recent studies (refs. 2, 5, 6, 13; but see refs. 4, 10), including an analysis of mitogenomes of 90 anuran species (7). Relationships among the early branches of Neobatrachia are identical to the mitogenomic phylogeny (7). We corroborate the placement of Heleophrynidae, which is from extreme southern Africa, as the sister taxon to all other neobatrachian frogs (Fig. 1A), as found by some authors (4⇓⇓–7, 9), but not by others (10, 11). We find that Sooglossidae (known only from the Seychelles Islands) is the sister taxon of Ranoidea (BS = 100%, BPP = 1.0; Figs. S1–S3) in the 164-species analyses, and that Sooglossidae + Nasikabatrachidae (known only from the Western Ghats of peninsular India) form the Sooglossoidea, which is the sister group of Ranoidea (BS = 100%; Fig. S4) in the 309-species analyses. Sooglossoidea is placed as the sister taxon of all other neobatrachians (11), all neobatrachian frogs to the exclusion of Heleophrynidae (6), Hyloidea + Myobatrachidae + Calyptocephalellidae (4), and Ranoidea (3, 5, 7), but no previous studies recovered the placement of Sooglossoidea with strong statistical support, including a mitogenome phylogeny (7).

Taxa within the superfamily Ranoidea have been consistently grouped into three clades (4⇓⇓–7): Microhylidae, Afrobatrachia (i.e., epifamily Brevicipitoidae, which includes Brevicipitidae, Hemisotidae, Hyperoliidae, and Arthroleptidae), and Natatanura (Fig. 1A). We found that Microhylidae is the sister group of Afrobatrachia + Natatanura (BS = 72%, PP = 1.0; Figs. S1 and S2); in contrast, many previous studies placed Microhylidae as the sister group of Afrobatrachia (4⇓–6, 14) or of Natatanura (7). Relationships within the Afrobatrachia mirror those found in other studies (5⇓–7, 14).

Natatanura is a large clade of extant anurans (24% of species) and mainly found in the Old World. Our ML and Bayesian topologies of Natatanura are identical. All nodes in the Bayesian tree have a BPP of 1.0, and only three nodes in the ML tree have BSs <90%. The 309-species topology is identical, but with low support among the deeper branches, likely because of missing data. Notably, we found that endemic African continental lineages (Conrauidae, Odontobatrachidae, Petropedetidae, Phrynobatrachidae, Ptychadenidae, Pyxicephalidae) form a clade that is the sister group to the clade of the remaining North American, Eurasian, Melanesian, and Malagasy lineages (Ceratobatrachidae, Dicroglossidae, Mantellidae, Rhacophoridae, and Ranidae; Fig. 1A and Figs. S1–S4). This African clade has low bootstrap support (56%) but high Bayesian support (1.0). The clade of the remaining non-African families is strongly supported (BS = 100%), and the internal branches are strongly supported (BS = 100%, BPP = 1.0), although they are short. In other studies, this group of African lineages is not monophyletic (6, 7, 10, 14, 15). The phylogenetic position of the continental African lineages has important biogeographic significance (as detailed later).

Relationships among the subfamilies of Microhylidae (one of the largest anuran families, including 8.8% of all species), which have significant radiations on most continents and the large islands Madagascar and New Guinea, have proven difficult to resolve (4, 6, 7, 15⇓⇓⇓–19). In contrast, our Bayesian tree has strong support; all 10 of the deepest nodes have a BPP of 1.0 (Fig. S2). The ML topology is identical, but 3 of 10 nodes have a BS <90% (Fig. S1). ML analysis of our 309-species dataset (Fig. S4) recovered the same topology, but support is weaker. Nonetheless, our tree is more strongly supported than others except for a phylogenomic analysis (19) of 66 anchored loci for 48 taxa and 7 Sanger-sequenced loci for 142 taxa.

Significantly, our study resolves relationships among one of the most diverse clades in the anuran phylogeny: Hyloidea, which contains 54% of extant anuran species. ML and Bayesian analyses strongly support 53 of 55 of nodes in Hyloidea (BS > 80% and PP = 1.0; Figs. S1 and S2). This is particularly unexpected because even the most species- and character-rich studies of hyloids (5⇓–7) have recovered poorly resolved topologies. Our expanded 309-species topology (Fig. S4) is similar to our primary tree, even though many nodes have lower support, which is consistent with the higher degree of missing data for this dataset of more species.

Two arrangements within hyloids are noteworthy. The deepest divergences in Hyloidea are among southern South American taxa: Rhinodermatidae in the temperate beech forests of Chile, and Alsodidae and Batrachylidae in Patagonia (Figs. S1–S4). Similar relationships were previously reported (5, 20), albeit with weak support and a smaller sample of species. These relationships support a southern South American origin of Hyloidea. Moreover, most previous studies (4, 6, 7, 21) supported Terrarana, a large New World tropical clade (15% of extant anuran species), rather than southern South America groups, as the sister group of all other hyloids. This difference in the placement of Terrarana may result from long-branch attraction of mtDNA sequences, because Terrarana species apparently have higher rates of mtDNA evolution than other hyloid lineages (22).

In summary, our analyses corroborate many of the deeper neobatrachian nodes inferred by other studies, but with greater support (all BS values = 100%). Furthermore, we find strong support for many shallower nodes that are weakly supported in other studies; only 4 of 155 nodes in our Bayesian tree have a posterior probability <1.0 and only 9 of 155 nodes have bootstrap support <75%. Relationships among the deepest branches in the Hyloidea, Microhylidae, and Natatanura, for which previous studies have conflicting results, are now strongly supported in Hyloidea and Natatanura, and most nodes are well-supported in Microhylidae. Some clades that were not well-supported in the mitogenomic phylogeny (7) are rejected in this study. Our findings include the position of southern South American taxa as the earliest branches of hyloids, and the identification of a clade of endemic African taxa as the sister group of all other Natatanura.

A New Timescale for Extinction and Diversification.

Our extensive sample of nuclear loci produced not only a strongly supported phylogeny for frogs but also younger divergence time estimates with smaller CIs than previous studies. The estimated divergence times in the 164-species and the 309-species data sets are similar (Figs. S5 and S6); the average time deviation among comparable nodes between the two analyses is ∼5%. Following a thorough literature review, we selected 20 fossil calibration points, of which 13 are within crown Anura (Fig. S7). Jackknife removal of each fossil resulted in highly congruent sets of date estimates. The overall average time differences were less than 1% except for the two outlier fossils (the salamander Chunerpeton and the salientian Triadobatrachus), and even the two outliers only led to overall average time differences of 3–5.4% (Fig. S8).

For this discussion, we focus on the 164-species tree (Fig. 1A). Overall, our divergence times are notably younger than those found by most other studies (5, 10, 23, 24) (Fig. 1A and Table S2), and especially much younger than those in the mitogenome phylogeny (7). For example, we estimate the last common ancestor of crown-group Anura to be during the Upper Triassic at 210.0 Mya (95% CI, 199.1–220.4 Mya), and the age of crown Neobatrachia in the Lower Cretaceous at 142.1 Mya (95% CI, 132.8–149.8 Mya). In contrast, many other estimates place the age of crown Anura between the Permian and Middle Triassic and the origin of Neobatrachia between Upper Triassic and Upper Jurassic (Fig. S9).

A striking pattern is the synchronous origin of three species-rich neobatrachian clades—Hyloidea, Microhylidae, and Natatanura—at the Cretaceous–Paleogene (K–Pg) boundary (KPB;Fig. 1A, dashed red line). The diversification synchronicity of the three frog clades still existed when all calibration constraints in frogs were excluded (Fig. S10), suggesting that this pattern is unlikely the result of the choice of calibration points. The K–Pg boundary (KPB), dated precisely at 66 Mya (25), marks one of Earth's great extinctions, largely attributed to the impact of the Chicxulub bolide. However, other factors such as climate warming and the Indian Deccan volcanism, possibly related to the bolide impact, likely contributed to this extinction event (26). Although a near-KPB origin was separately reported for Hyloidea (5, 8, 10, 20, 27) and Microhylidae (10, 16) (Fig. S9), our study found that these three major clades originated at the same time, very near the KPB. The contemporaneous origin of the three large clades is highlighted in our results by the narrow 95% CIs for each and their overlap with the KPB. Nine of the 10 deepest nodes of Hyloidea have relatively narrow CIs that overlap the KPB (9.6–14.0 My; Fig. 1A and Table S2). Similarly, within Microhylidae, the CIs of the four deepest nodes encompass the KPB, and the CIs are similarly narrow (11.1–11.9 My). Finally, in Natatanura, the CIs of the six deepest nodes all overlap the KPB, also with similarly narrow CIs (10.4–10.7 My). The use of internal calibration points for these clades is unlikely to be driving this pattern because there was only one internal calibration point (within the Natatanura) used within these three clades.

There are several reasons why our study may have produced younger divergence times compared with previous studies. First, we included diverse outgroups, including lungfish, coelacanth, salamanders, caecilians, and amniotes. Second, we used a new set of calibration points rather than uncritically recycling calibration points from other studies; these calibration points have been carefully reexamined based on the fossil record. The choice of calibration points has significant impact on divergence time estimation for frogs and may be an important reason that older divergence times were obtained in previous studies (27). Third, our phylogenetic estimate is based on a much larger dataset (88 kb; 62% of sites are variable) derived solely from 95 nuclear loci (rather than primarily from mitochondrial loci) that exceeds other studies by at least sevenfold. We attribute the greater precision (i.e., smaller CIs) to the large amount of information within our dataset. A similar increase in precision was found for plethodontid salamanders by using the same nuclear marker toolkit (28).

Rate-through-time analysis indicates that a surge in net diversification rates of frogs occurred immediately following the KPB (Fig. 1B), which suggests a clear impact of the KPB extinction on frog diversification. The rapidity of the diversification is reflected in the short branches of the deepest nodes in each clade. In Hyloidea, all but one branch connecting the 10 deepest nodes are <4.6 My in duration. In Microhylidae, the branches connecting the four deepest nodes are <4.0 My long. In Natatanura, the branches connecting the six deepest nodes are <2.1 My. Even though these branches are short, they are strongly supported. Another line of evidence demonstrates three mass extinctions preceding explosive speciation. The branches subtending the Hyloidea, Natatanura, and Microhylidae are long: 47.3, 32.1, and 33.8 My, respectively (Fig. 1A and Table S2). The lack of other extant lineages originating from these three stem lineages corroborates our suggestion that an extinction event simultaneously decimated these lineages near the KPB.

Even though the fossil record indicates mass extinction at the KPB for several taxa, including birds, squamates, and mammals, and subsequent radiation post-KPB (29⇓–31), molecular data often indicate that many clades had diversified, at least in the sense of lineage splitting, before the KPB event (32⇓⇓⇓–36). Unlike nonavian dinosaurs, whose demise is well documented in the rocks, the fossil record of frogs is so far largely uninformative about survival/extinction across the KPB (37⇓⇓–40). In the present study, the three parallel combinations of precise node ages overlapping the KPB, preceded by a long branch and followed by very short but strongly supported branches, indicate three extinctions followed by rapid divergence. These three radiations, which are coincident with the late Cretaceous extinction, account for 88% of extant frog species. This molecular phylogenetic perspective strongly suggests that frogs experienced a major extinction at the KPB.

Historical Biogeography of Frogs.

We performed a biogeographic analysis on the 309-species data set using BioGeoBEARS (41). The DEC+J model performs significantly better than the DEC model (Akaike information criterion; Table S3), indicating the importance of the J parameter, which models long-distance or “jump” dispersal. Our interpretation is that dispersal into a new area is accompanied by near-instantaneous speciation. The practical effect of adding the J parameter to the model is that ancestral ranges often comprise one area rather than several.

The most recent common ancestor (MRCA) of extant frogs was distributed in Eurasia + North America + Australia or North America + Australia (Fig. 2 and Fig. S11), a Pangaean origin. In contrast, neobatrachians (i.e., modern frogs) originated in Gondwanaland, most likely in Africa (relative probability, 93.6%; Fig. 2). The major neobatrachian lineages are also Gondwanan in origin: Hyloidea in South America and Ranoidea, Afrobatrachia, and Natatanura in Africa (relative probability, >90%; Fig. 2).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Ancestral-area estimates for 69 terminal taxa (families, subfamilies, and genera) of extant frogs using the DEC+J model in BioGeoBEARS. Circles on nodes represent the set of possible ancestral areas, and the color is associated with the area legends. The probabilities are given next to circles for the most probable ancestral area. Circles without values indicate that the probability of the ancestral area is >99%. Three important landmass breakup events are indicated: (A) the break-up of Pangea with division into Laurasia and Gondwana in the Middle Jurassic coincident with the origin of neobatrachian frogs; (B) the separation of Africa and South America in the Early Cretaceous coincident with the divergence of Ranoidea and Hyloidea as well as between the African and New World pipids; and (C) the separation of the Seychelles and India in the Late Cretaceous coincident with the divergence between the Sooglossidae and Nasikabatrachidae. Anuran taxa that contain at least some arboreal species are indicated in green and a tree icon (we do not imply that the last common ancestor of each of these families was arboreal). Note that all clades containing arboreal frogs originated after the KPB.

Breakup of land masses may be associated with three divergence events in the evolution of frogs. The first of these was the split of neobatrachians from ancestral groups (Fig. 2A), dated at 172–181 Mya, which agrees with the breakup of Pangaea into Laurasia and Gondwana in the Early Jurassic (∼180 Mya) (42). This breakup event, which isolated the MRCA of Neobatrachia in Gondwana, is associated with the initial diversification of the clade. The second breakup event caused the split between the two major lineages of neobatrachians, Procoela (Hyloidea and Myobatrachoidea, mostly South American) and Diplasiocoela (Ranoidea + Sooglossoidea, mostly African), as well as the split between the African and South American pipids (Fig. 2B). Spreading of the South Atlantic Ocean sea floor started in the Early Cretaceous (135 Mya), but the final physical separation between Africa and South America took place approximately 105 Mya (42). Our estimates for timing of these divergences, 125–130 Mya and ∼120 Mya, respectively, are highly congruent with this rifting process. The third breakup event took place when the Sooglossidae (Seychelles) and the Nasikabatrachidae (India) diverged (Fig. 2C). The Seychelles/India land mass continued to exist until 66 Mya, when new rifting severed the Seychelles from India (43). We estimate Sooglossidae split from Nasikabatrachidae at ∼66 Mya, which is considerably younger than the previous estimates (77–130 Mya) (10, 11, 23), in congruence with the geological event.

The MRCA of Hyloidea and of Myobatrachidae + Calyptocephalellidae occurred in South America (Fig. 2). The divergence between Myobatrachidae (Australia) and Calyptocephalellidae (South America) and the split between Phyllomedusinae (South America) and Pelodryadinae (Australia and New Guinea) took place ∼100 Mya and ∼50 Mya, respectively (Fig. 2). During this period, South America and Australia were distantly separated but connected intermittently via Antarctica (42). From the late Cretaceous to the early Tertiary, the Earth experienced a period of global warming (44). Climate data from plant fossils, sediments, and geochemical indicators show that the mean annual temperature of the Antarctic Peninsula region was 10–20 °C from 100 to 50 Mya (45), which is sufficiently warm to allow dispersal of frogs through this region. Thus, the origin of the Australian myobatrachids and pelodryadine hylids is most likely explained by dispersal from South America to Australia through Antarctica, and later extinction in Antarctica because of the formation of ice sheets. The role of Antarctica in dispersal of frogs among Gondwanan land masses has received little attention (but see ref. 46), but our time-calibrated phylogeny predicts that paleontological research in the Cretaceous and early Cenozoic of Antarctica could shed light on the early evolution and dispersal, especially of hyloids.

The nearly worldwide distribution of Microhylidae presents a longstanding and challenging biogeographic puzzle. Vicariant origin and long-distance oceanic dispersal have been proposed to explain the current distribution of this family (10, 15, 16, 18, 23). Which scenario dominates interpretation depends largely on the time and location of origin for the Microhylidae. Our analyses place the initial divergence of the Microhylidae at ∼66 Mya, with a relative probability of 79.5% for an African origin (Fig. 2). By this time, Gondwana was already highly fragmented and Africa was separated from South America and Madagascar by ocean. Accordingly, overseas dispersal of the major microhylid lineages on Gondwanan landmasses is required.

The origin and diversification of the second large clade, Natatanura, is controversial. Two hypotheses have been proposed. The Out-of-Africa hypothesis (47), argues that the origin of the Natatanura lies in Africa, with subsequent dispersal to other continents. The Out-of-India hypothesis (48) postulates that the ancestor of Natatanura originated on the Indian plate and was confined there until the plate collided with Asia approximately 55 Mya. We find that Natatanura originated in Africa (relative probability, 92.5%; Fig. 2). Considering that Natatanura is the sister group of Afrobatrachia (apparently of African origin) and the basal split within the Natatanura phylogeny separates all African continental lineages from other Asian, Indian, and Madagascar lineages (Fig. 2), an African origin seems more likely. The endemic Indian natatanurans (Nyctibatrachidae, Micrixalidae, and Ranixalidae) are nested within Asian lineages, and the divergences between them and their closest Asian relatives occurred between 55 and 60 Mya (Fig. 2), which is consistent with the timing of the India–Asia collision (∼55 Mya) (49). Intriguingly, the divergence between the Asian tree frogs (Rhacophoridae) and the Malagasy mantellids also happened during this period (Fig. 2), implying that the Indian plate served as a stepping stone for the long-distance dispersal from Asia to Madagascar.

The three long branches of similar duration leading to Hyloidea, Microhylidae, and Natatanura suggest that the K–Pg mass extinction event may have triggered explosive radiations by emptying ecological space. Notably, all terminal taxa (family and subfamily ranks) with arboreal species originate after the K–Pg boundary (Fig. 2). The rebounding of forests after the massive loss of vegetation (50) at the KPB may have provided new ecological opportunities for the subsequent radiation of largely arboreal groups such as Hylidae, Centrolenidae, and Rhacophoridae. The observation that no lineages of frogs originating before the KPB (archaeobatrachians, Heleophrynidae, Myobatrachoidea, and Sooglossoidea) have truly arboreal species, and that all origins of arboreality (e.g., within hyloids or natatanurans) follow the KPB, supports the hypothesis that the K–Pg mass extinction shaped the current diversity of frogs.