Reg. Charity No. 1168330
Anurans have a long fossil record, spanning from the Early Jurassic to the Recent. However, specimens are often severely flattened, limiting their inclusion in quantitative analyses of morphological evolution. We perform a two‐dimensional morphometric analysis of anuran skull outlines, incorporating 42 Early Cretaceous to Miocene species, as well as 93 extant species in 32 families. Outlines were traced in tpsDig2 and analysed with elliptical Fourier analysis. Fourier coefficients were used in MANOVAs, phylogenetic MANOVAs (as significant phylogenetic signal was found) and disparity analyses across multiple ecological and life history groupings. The Neotropical realm showed higher disparity than the Australian, Palaearctic and Oriental realms (p = 0.007, 0.013, 0.038, respectively) suggesting concordance of disparity and diversity. Developmental strategy had a weak effect on skull shape (R2 = 0.02, p = 0.039) and disparity was similar in metamorphosing and direct developing frogs. Ecological niche was a significant discriminator of skull shape (F = 1.43, p = 0.004) but not after phylogenetic correction. Evolutionary allometry had a small but significant influence on the cranial outlines of the combined extant and fossil dataset (R2 = 0.05, p = 0.004). Finally, morphospace occupation appears to have changed over time (F = 1.59, p = 5 × 10−10). However, as with ecological signal, this shift appears to be largely driven by phylogeny and was not significant after phylogenetic correction (R2 = 0.26, p = 0.22). This study thus suggests that frog skull evolution is shaped more by phylogenetic constraints than by ecology.
Anurans (frogs), both extant and extinct, are the most speciose and diverse clade of Lissamphibia. Currently numbering 6955 species (AmphibiaWeb 2018) frogs are found in a variety of ecological niches with a near worldwide distribution. Anurans have a long fossil history with impressive past diversity, ranging from the monstrous Beelzebufo ampinga (Evans et al. 2008) to the tiny Indobatrachus pusillus (Owen 1847). Salientia, stem and crown group anurans, are present in the fossil record from the Early Triassic (e.g. Triadobatrachus, c. 250 Ma; Piveteau 1936; Ascarrunz et al. 2016) but probably originated in the Permian (Ruta & Coates 2007). A greater abundance of fossils is seen from the Early Jurassic, including the stem‐anurans Prosalirus (Shubin & Jenkins 1995) and Viaerella (Estes & Reig 1973). The first crown group anuran is Eodiscoglossus oxoniensis dating from the Middle to Late Jurassic and may be the earliest member of the Discoglossoidea (Evans et al. 1990). Eodiscoglossus is represented by many Early Cretaceous fossils. It is likely that some other modern clades had evolved by the start of the Cretaceous period, including Pipidae (e.g. Estes et al. 1978). The majority of fossil occurrences are based on isolated or disarticulated bone fragments (see Sanchíz (1998) for a review) but rare cases of exceptional preservation including soft tissue are known, including Eleutherodactylus from the Eocene amber of the Dominican Republic (Poinar & Cannatella 1987) and a mummified specimen of Thaumastosaurus gezei from 30 Ma (Laloy et al. 2013). Despite the remarkable preservation of some fossil anurans, three‐dimensionally preserved specimens remain rare. The vast majority of complete specimens are severely flattened, which may be partially attributable to their thin, rod‐like cranial bones. Consequently, morphological evolution has not been well studied across fossil frogs, as morphological studies are limited primarily to descriptions of individual specimens (see Estes & Reig 1973) and not comparisons across taxa.
Preservational and sampling biases can strongly influence the completeness of a clade's fossil record. Temnospondyls, with their generally large size and robust skulls, have an extensive fossil record, with around 300 species identified (Schoch 2013) and a near global distribution. Temnospondyl diversity can be tracked from the Lower Carboniferous to the Lower Cretaceous (McHugh 2012) and includes many complete skulls. However, their descendants, in particular the anurans, are considerably more affected by preservational biases, being generally small organisms in ecosystems which do not preserve readily, such as tropical forests (e.g. Duellman 1999). In fact, despite a present day species ratio of 1:3, there is a mammal to amphibian fossil record ratio of >5:1 during the Paleocene (Benson et al. 2016). This is probably a combination of less intensive sampling efforts, less robust skeletons of amphibians, and many fossil mammals named based on isolated teeth, and highlights the difficulty of inferring past diversity from present. Throughout their history there is a clear bias towards amphibians preserved from freshwater systems compared with terrestrial systems (Schoch 2014, ch. 10) but both systems preserve considerably less readily than marine settings. There is also disproportionate representation of amphibian fossils from Lagerstätten, creating a biased picture of past amphibian diversity (Schoch 2014, ch. 10). The under‐represented nature of amphibians in the fossil record therefore stresses the importance of incorporating all suitable fossils into morphological studies.
The quantification of morphology allows us to explore how environmental, developmental and functional pressures can influence organismal shape. Recent decades have seen many advances in the study of morphology, in particular in the field of three‐dimensional shape analysis, with anatomy captured using a range of approaches from anatomical landmarks (e.g. Goswami 2006) to entire surfaces represented by sliding semi‐landmarks (described by Bookstein 1997; Gunz et al. 2005; and implemented by e.g. Fabre et al. 2014; Felice & Goswami 2018). However, these methods require high quality data if subsequent analyses are to be robust. Specimens must have well‐preserved and undeformed morphology, which are highly constraining requirements for most fossil specimens. Thus, many studies choose to exclude fossils rather than compromise the strength of the morphological analysis. Many macroevolutionary studies have explored morphological evolution across large clades, including Gymnophiona (Sherratt et al. 2014), Aves (Cooney et al. 2017; Felice & Goswami 2018) and the anuran family Myobatrachidae (Vidal‐García & Keogh 2017), but these clades (and many others) suffer from a fossil record largely inadequate for inclusion in comprehensive, quantitative morphological studies. Past diversity has been studied in clades such as temnospondyls (Stayton & Ruta 2006; Angielczyk & Ruta 2012) but these datasets represent clades with unusually abundant fossils available for study. Nonetheless, these fossil studies highlight the rich and varied past diversity of clades such as amphibians. This then leads to the inevitable decision of whether to prioritize the capture of morphological data (subsequently sacrificing unsuitable specimens) or use a more inclusive approach, reducing sampling of morphology but representing more of past and present diversity.
The cranium is a particularly informative osteological structure for inferring anuran ecology and behaviour (e.g. Jared et al. 2005; Senevirathne et al. 2016). Whilst some consider the anuran limb system to be the most unique feature of this clade, anuran skulls exhibit extreme morphological variation, from hypo‐ossification (e.g. Bombina) to heavily ossified and ornamented skulls (e.g. Ceratophrys), as well as the presence of novel cranial bones (e.g. the prenasal of the casque‐headed Triprion petasatus). Previous studies have explored three‐dimensional morphological variation in the anuran cranium (e.g. Vidal‐García et al. 2014; Vidal‐García & Keogh 2017) but fossils have not been included in these studies as three‐dimensional specimens are exceptionally rare. However, biological information can still be extracted from two‐dimensional anuran cranial morphology, especially as frog skulls are reasonably dorsoventrally compressed relative to the crania of other tetrapods (Trueb 1993, p. 311). Pipoid and non‐pipoid frogs have been compared using two‐dimensional cranial and lower mandible landmarks (Yeh 2002a), which revealed that pipoids differ ontogenetically from salamanders, but are more similar to them morphologically than non‐pipoids are. Two‐dimensional cranial landmarks were used by another study (Yeh 2002b) to compare skull shapes of miniaturized frogs and closely related larger species. Miniaturized frogs were found to have relatively larger braincases and a reduction of some bone elements. Despite these promising studies demonstrating that two‐dimensional data successfully discriminates skull shape in extant frogs, this approach has not yet been applied to fossil frogs.
Here, we quantify the two‐dimensional morphology of extant and fossil anuran skulls to take advantage of the many exceptionally preserved yet severely flattened fossil specimens. Because of the difficulty in consistently identifying homologous landmarks in many fossil frog specimens, we apply elliptical Fourier analysis to the outlines of anuran skulls. Using these data, we assess whether phylogeny, development or ecology is a greater influence on anuran skull morphology, and quantify how anuran skull morphology changes through time. First, we evaluate ecological and developmental signals in extant frog skull shape and compare disparities across groupings based on these factors. We then combine data from the extant and fossil datasets to measure phylogenetic signal and evaluate the influence of size on skull shape across all frogs. Finally, we examine shifts in frog skull shape through time and differences in disparity among temporal groupings of frogs in order to provide a comprehensive analysis of frog skull shape evolution, as quantified by two‐dimensional outlines.
Material and methodPhotos, diagrams, and reconstructions of 93 extant and 65 fossil frog skulls were obtained either from the published literature or from original specimens (Bardua et al. 2018, tables S1, S2; see Fig. 1 for an example of a fossil frog). Phylogenetic and ecological diversity is well sampled within the dataset, with 32 extant families represented, ranging from the aquatic Pipidae to the fossorial Pyxicephalidae, and from the tiny Eleutherodactylus pipilans (snout–vent length (SVL) reaching 29.4 mm; AmphibiaWeb 2018) to the world's largest frog, Conraua goliath (SVL up to 320 mm; AmphibiaWeb 2018). This study also has temporal coverage spanning most of anuran evolution, with 65 specimens of 42 fossil species ranging in age from the Middle‐Late Jurassic stem‐anuran Notobatrachus degiustoi (176–145 Ma, Báez & Nicoli 2004) to fossil members of living genera, such as early Middle Miocene Rana (Roček et al. 2011). Adult specimens were chosen to the best of our ability, but it can be difficult to discern between adult and sub‐adult fossil specimens. Fossil specimens were included if at least one bilateral half of the skull was traceable, and this traceable half was mirrored when required (Bardua et al. 2018, table S2). Initially, multiple specimens of each fossil species were included in the main dataset to assess the effect of deformation qualitatively, but the final analyses used only the least deformed specimen or best reconstruction for each species. Subsets of this dataset were used for different analyses, depending on the ecological, life history and body size traits available for each species, as detailed below. Taxa with unknown or ambiguous phylogenetic position were removed from phylogenetic analyses.
Figure 1 Open in figure viewerPowerPoint Aerugoamnis paulus FMNH PR2384 from the Early Eocene of Green River Formation, Wyoming, demonstrating the dorsoventral flattened preservation typical of specimens in this study (photo from Henrici et al. (2013) used with permission; FMNH, Field Museum of Natural History, Chicago, USA). Scale bar represents 3 mm. Colour online.Outlines were taken as the largest tracing around the skull, which in some cases meant that ventral attributes were being traced in dorsal aspect. However, in most cases, it was not possible to discern details beyond the largest outline. The program tpsUtil (Rohlf 2013) was used to build a tps file from images, and this file was loaded into tpsDig2 (Rohlf 2015). The outline of each image was traced using the closed curve function in tpsDig2 and this curve was then resampled to 100 points. The starting point of each curve was the anteromedial extreme of the skull, which was oriented to face upwards. Individual points were then slid around the outline, as points did not need even spacing because the analyses conducted on the outlines only looked at the outline not the individual points. Consequently, areas of high detail could be represented by more points than areas of conserved curvature (Fig. 2). Some fossils did not have scales included in the published images, but this did not affect their inclusion in analyses because scale, as well as position, are non‐shape aspects that were removed prior to analysis. Once all images had been outlined in the tps file, the curves were appended to landmarks using tpsUtil. These landmarks were then imported into R (R Core Team 2017) for analysis.
Figure 2 Open in figure viewerPowerPoint The outline of Adelotus brevis KU 56242, an extant specimen used in this study (see Bardua et al. 2018, table S1). A, traced using tpsDig2 and resampled to 100 curve points. B, generated using the Momocs package, using the 100 resampled curve points. (KU, University of Kansas, Museum of Natural History, USA). Scale bar represents 20 mm.The phylogenetic tree used for the extant frog dataset was adapted from the molecular phylogenetic analysis and divergence estimates of 2871 amphibians (Pyron & Wiens 2011). This was the most comprehensive phylogeny available for living amphibians. The R packages geiger (Harmon et al. 2008), ape (Paradis et al. 2004) and phytools (Revell 2012) were used to prune this tree and add 28 extant species by assigning them to the positions of their respective genera. Polytomies were randomly resolved in Mesquite (Maddison & Maddison 2016). Six species were excluded from phylogenetic analyses due to uncertain phylogenetic positions, including Crossodactylodes pintoi and Zachaenus parvulus which are both assigned to Hyloidea incertae sedis (Pyron & Wiens 2011).
Fossil anuran species in our study were manually added to the extant phylogeny for our combined analyses (Fig. 3). The topology for the fossil species was based on a recent composite tree analysis (Marjanović & Laurin 2014), which collated many fossil lissamphibian species and placed them within the framework of the extant amphibian phylogeny we used (Pyron & Wiens 2011). Incongruences exist between the Marjanović & Laurin (2014) tree and more focused phylogenetic studies; for example, Beelzebufo ambinga has been variably placed within Calyptocephalidae (Agnolín 2012; Marjanović & Laurin 2014) or Ceratophryidae (Laloy et al. 2013; Evans et al. 2014). However, as the Marjanović & Laurin (2014) tree has considered all major attempts to place fossil lissamphibians into trees and is congruent with the extant topology used in this study (Pyron & Wiens 2011) this composite tree topology was further modified for the combined analyses.
Figure 3 Open in figure viewerPowerPoint Composite phylogenetic tree of sampled extant and extinct anurans. Alternate shaded bands indicate periods of 50 million years, with the outer extreme at Recent.A composite tree was created by manually inserting extinct species into our extant phylogeny, using branch lengths that we first calculated by creating a topology of just the extinct species. To create the extinct species topology, the positions of 28 extinct species in our study were exactly replicated from the Marjanović & Laurin (2014) tree. Eight species (mainly Liaobatrachus and Rana specimens) were not included in that work, but congeneric species were represented, so they were placed within their respective genera, resolving polytomies at random. For taxa without representatives in the Marjanović & Laurin (2014) tree, placement was estimated based on other recent studies. These taxa include Thaumastosaurus genei (Laloy et al. 2013) and Aerugoamnis paulus (Henrici et al. 2013). Four specimens with unknown phylogenetic placement were removed from the phylogenetic analyses: Eupsophus sp., Lutetiobatrachus gracilis and two specimens of Anura indet. This topology was then dated using the R package Paleotree (Bapst 2012). The first and last occurrence dates of each fossil taxon were taken from the literature, and each was placed into 5 million year time bins. The sampling rate was set at a low value of 0.01, as much of the frog fossil record consists of isolated bones that could not be included in the study. Speciation and extinction rates were assumed to be equal (Foote et al. 1999). The cal3 dating method (Bapst 2013) created 1000 dated trees; the average branch lengths from these trees were then used to create the composite tree.
To combine our extant and fossil phylogenies, fossil taxa were manually added to the cropped extant phylogeny, using their calculated branch lengths as guidance. However, specific branch lengths could not easily and completely be superimposed due to some variation in topology, so some nodes required manual manipulation in Mesquite to ensure all branches and taxa were aged as accurately as possible with the current information on phylogenetic relationships and occurrences.
Ecological data were compiled from the literature and online databases (IUCN 2016; Myers et al. 2017; AmphibiaWeb 2018) on ecological niche, environment, zoogeographical realm occupation, IUCN Red List Status and developmental strategy (Bardua et al. 2018, table S1). Ecological niche (N = 86) was divided into terrestrial (ground dwelling), fast aquatic (rivers and streams), slow aquatic (lakes and ponds), semi‐aquatic, fossorial (underground), arboreal and semi‐arboreal, and species with broad niches or no available niche information were excluded. Environment (N = 66) was defined as: arid open, temperate forest, temperate open, tropical forest and tropical open. Frogs inhabiting a range of environments were excluded from this analysis. The zoogeographical realm (N = 85) was divided into the following: Afrotropical, Australo‐Papuan, Nearctic, Neotropical, Oriental, Palaearctic and Widespread (Duellman 1999). ‘Widespread’ was defined as species found in more than one realm. Extinction risk (N = 78) was classified by the Red List Index (IUCN 2016) with the following categories: 1 (Least Concern), 2 (Near Threatened), 3 (Vulnerable), 4 (Endangered), 5 (Critically Endangered). ‘Data deficient’ specimens were removed from this analysis. Extinction risk was treated as an ordered factor in analyses. Finally, developmental strategy (N = 78) distinguished between species that undergo metamorphosis (indirect) or don't (direct).
Snout–vent length (SVL) was taken as the measurement for body size as it was the most readily available measurement for fossil frogs and is the standard body measurement for extant frogs (see from Roček et al. 2010, fig. 17). Fossil frog SVL was taken from the literature as the measured or estimated SVL of the specimen, or the maximum SVL found for a specimen of that species. SVL for extant frogs was taken as the maximum female SVL, as females are generally the larger sex in frogs (Shine 1979), or the upper limit of any SVL range given. SVL data were available for 112 of the sampled species, 29 of which were fossil frog taxa and 83 were extant frog taxa.
Two‐dimensional projections of three‐dimensional shapes can be effective correlates of overall shape, potentially introducing only small error into the data (for a review see Cardini 2014). Two‐dimensional shape can be described by landmarks (e.g. Wund et al. 2012), linear measurements (e.g. Jojić et al. 2012) or by outline (e.g. Zhan & Wang 2012). It is generally accepted that outline analysis is preferable for capturing shape information when homologous landmarks are missing (see Temple 1992). Outline approaches include eigenshape analysis (Lohmann 1983), polar Fourier analysis (e.g. Kaesler & Waters 1972) and perimeter‐based Fourier analysis (Foote 1989) but the most common approach is elliptical Fourier analysis (see Rohlf & Archie 1984; Temple 1992). The Fourier series was first developed by Joseph Fourier (1768–1830), who demonstrated that a complex waveform can be expressed as a sum of trigonometric functions. This approach was later developed to apply this mathematics as a tool for quantifying shape in biological studies. Elliptical Fourier analysis (EFA) uses the Fourier series to describe the shape of a closed curve (Giardina & Kuhl 1977; Kuhl & Giardina 1982). Decomposing a complex, closed waveform into a series of simpler, sinusoidal waves of varying period and magnitude allows us to understand the primary aspects of shape, as well as reducing dataset dimensionality. EFA has been used to study shape in a variety of subjects, from assessing shape differences in forest fragment populations of leaves (Andrade et al. 2010) to comparing water basin shapes (Bonhomme et al. 2013). It has also been used in morphological studies, from comparing bird caudal skeletal morphology across the ecological spectrum (Felice & O'Connor 2014) to assessing artificial cranial deformation in human crania (Frieß & Baylac 2003). Thus, EFA has been shown to be a successful method for assessing two‐dimensional shape in a range of akinetic biological shapes.
Outline analysis was used to capture the shape of anuran skulls. Unlike other outline approaches, EFA does not require landmarks or a biological centre to the outline which is homologous throughout the sample (Crampton 1995). In addition, evenly spaced points are not a requirement, allowing increased sampling from areas of high curvature to better represent the full outline. Elliptical Fourier (EF) coefficients provide an efficient measure of the major aspects of shape to a desired degree of precision, reducing the dimensionality of the dataset by removing unnecessary harmonics (as determined by those cumulatively representing less than a set proportion of the shape).
Tps files were loaded into R and analysed using the Momocs package (Bonhomme et al. 2014). The data was firstly made into a ‘Coo’ object, which is the format required to analyse the outlines. These outlines were then centred and scaled to remove all non‐shape data from the analyses. Elliptical Fourier analysis was applied to the outlines using the efourier function in Momocs. Normalization was not performed in this step, as the outlines (being symmetrical and roughly circular) are at risk of bad alignment. Orientation effects had already been considered by orienting all images with the snout facing upwards. The Fourier coefficients were then extracted to use as high dimensional variables in analyses. Each harmonic n is represented by four coefficients, two from each of the x and y projections, so that the total number of coefficients is 4n. The number of harmonic coefficients required for the outline analysis was determined by choosing the number which corresponded to a cumulative sum harmonic power of 99%. We therefore used 14 EF coefficients for the extant dataset and 12 for the fossil dataset. Although the number of extra harmonics nearly doubled for 99% compared to 95% cumulative power, the extra harmonics were necessary to accurately represent all the highly detailed information in the posterior region of the skull. This reduces the dimensionality of the dataset down from 200 (100 two‐dimensional coordinates) to 56 and 48 dimensions for the extant and fossil datasets respectively.
To further reduce dimensionality, principal components analysis (PCA) was conducted in the R package Momocs on the EF coefficients which represented 99% cumulative power. Principal components analysis converts a dataset into a set of mathematically uncorrelated variables, which represent the principal axes of shape variation in a given dataset. This allows a large dataset to be represented by far fewer variables, making multivariate analyses tractable to interpret. Principal components should not be analysed separately as univariate variables, and a sufficient number of principal components should be retained, as using just the first few components can introduce bias into the analysis (Uyeda et al. 2015). Principal component (PC) scores that described 99% of the shape were retained and used in subsequent analyses.
We first conducted a series of sensitivity analyses to assess the effects of pooling photos, figures and reconstructions, and including specimens preserved in different aspects or with different degrees of completeness. Some diagrams of fossil specimens depicted the outline more clearly than a photo. To test whether diagrams from the literature are accurate representatives of the true skull shape, PC scores of skull outlines from diagrams and corresponding photos of 33 species were compared using the MANOVA function in Momocs. As well as diagrams, many authors reconstructed the skull in two‐dimensions, using one or multiple specimens as a reference (e.g. Báez & Nicoli 2004). As these reconstructions attempt to eliminate the deformation present in the original specimens, they may provide a better representation of the true skull shape. To assess the validity of using reconstructions, all specimens along with the reconstruction for each species were outlined and compared to see whether reconstructions plot in the same region of morphospace as the specimens. In addition, specimens of fossil frogs commonly exist either in dorsal or ventral aspect but not always as both part and counterpart. To test whether aspect is an important factor for skull outline shape variation, photos of extant frog skulls from 27 species were outlined in both dorsal and ventral view. As noted above, outlines were traced as the largest outline as it would look from a shadow projection and not by following bones from purely the dorsal or ventral side, as individual bones could not be determined in most fossil frog images.
Finally, deformation and crushing evident in fossil frog specimens introduces distortions in quantification of the true skull shape. The posterior outline of the fossil frog skull generally appears worst affected, probably due to the higher‐level detail here which is more easily affected by crushing. In particular, the squamosal and pterygoid regions were occasionally difficult to trace, and the occipital section was frequently difficult to decipher because of the overlapping of vertebrae. To assess the sensitivity of results to inaccuracy in outlining the posterior skull, a sample of 32 species was outlined three times, firstly completely (D1), secondly by deliberately removing information about the posteromedial section and replacing this with a straight‐line perpendicular to the anteroposterior plane (D2), and thirdly by similarly removing the posterolateral information (D3) (Fig. 4). The D2 outlines therefore exclude the occipital region, and the D3 outlines exclude the pterygoid region. PC scores from the three outlines for each specimen were then compared through pairwise MANOVAs. We can therefore assess the degree to which uncertainty in the outline tracings can affect the analysis. These extreme examples of incomplete outlines create realistic maximum values for the difference between true and observed outlines.
Figure 4 Open in figure viewerPowerPoint The outline of Adelotus brevis KU 56242 adjusted for sensitivity analyses assessing the effect of missing or degraded posterior morphology on results. A, outlined completely (D1). B, deliberately removing posteromedial information (D2). C, deliberately removing posterolateral information (D3). Scale bar represents 20 mm.Phylogenetic signal of the PC scores was evaluated using the Kmult statistic (Adams 2014), which is a multivariate generalization of the K statistic (Blomberg et al. 2003), implemented in the geomorph package (Adams & Otárola‐Castillo 2013) in R. Closely related species are likely to have similar traits due to shared ancestry, and this signal of phylogeny is measured relative to expectations under a model of Brownian motion. Non‐independence of data points means that once phylogenetic signal is detected, it is necessary to apply a phylogenetic correction to the shape data to investigate evolutionary differences.
Ecological and developmental data from extant species only were used in MANOVAs (MANOVA function in Momocs) and phylogenetic MANOVAs (herein pMANOVAs; procD.pgls function in geomorph), to assess multiple potential factors influencing skull shape. Whilst these two tests differ in their assessment of significance (the procD.pgls function is non‐parametric and the MANOVA function is parametric) results are comparable. The Pillai–Bartlett trace was used as it is considered to be the most robust test statistic (see Johnson & Field 1993). A correction was applied to the MANOVAs and pMANOVAs to account for an elevated false‐positive rate with increased number of tests (Benjamini & Hochberg 1995). Morphological disparity (as measured by Procrustes variance) was also compared across ecological groupings, using the morphol.disparity function in geomorph.
Allometric shape changes for the combined fossil and extant dataset were visualized by plotting the predicted Procrustes residuals at minimum and maximum size using the procD.allometry function in geomorph. Specimens were also plotted on a morphospace, colour coded by SVL, to display the distribution of cranial sizes along the main axes of shape variation. The influence of evolutionary allometry (shape correlates of size among species) was assessed using pMANOVAs (procD.pgls function in geomorph) for both the extant‐only and the combined extant and fossil dataset, with log‐transformed SVL set as a factor against PC scores. Log‐transformed size variables better fitted a normal distribution.
To assess differences in morphospace occupation based on species age, a MANOVA and pMANOVA were applied to PC scores using species age as a factor. For fossil frogs, this is the first occurrence date within the fossil record; for extant species this is the first occurrence date of that species if it has a fossil record, or ‘Present’ if no fossil record exists. The full dataset was binned into time bins of 5 myr.
ResultsElliptical Fourier analysis successfully captured the frog skull outlines in terms of a series of ellipses which reduced the dimensionality of the dataset down to 12 (fossil frogs) or 14 (extant frogs, and combined dataset) harmonics, to an accuracy of 99% cumulative power (Bardua et al. 2018, fig. S1). Seven harmonics would have been sufficient to describe 95% of the shape in both extant and fossil frogs but the last 4% of shape information was important for capturing appropriate detail in the posterior region of the skull.
Principal component scores were retained in each analysis so that 99% of the shape was described, resulting in 22–24 PC scores retained for extant frog analyses (depending on the subset of data) and 15–17 for the fossil subsets. For the combined dataset of extant and fossil frog outlines, 25 PC scores were retained. The first principal component of the combined dataset accounted for 37.7% of the variation, and PC2 accounted for a further 18.1%. Similar values were found for the extant (PC1 = 38.3%, PC2 = 18.3%) and fossil (PC1 = 40.5%, PC2 = 19.4%) only datasets. In the combined extant and fossil dataset, PC1 corresponds to the relative length ratio of the skull; one extreme shape is anteroposteriorly longer than it is wide, whereas the other extreme is wider than it is long (Fig. 5). The negative extreme of PC2 corresponds to a more triangular outline, whereas the positive extreme is more circular.
Figure 5 Open in figure viewerPowerPoint Morphospace occupation of 125 frog skull outlines through time. Points are coloured from purple (youngest) to red (oldest) in 5 myr bins based on species first appearances. Extreme outline shapes representing the positive and negative extremes along PC1 and PC2 are displayed.Sensitivity analyses supported the pooling of all data types, and allowed us to visualize the information lost from the posterior region of the skull. Dorsal and ventral outlines were shown to be not significantly different (MANOVA, approx. F = 0.42, p = 0.98). They were also shown to have similar morphological disparity as illustrated by comparing Procrustes variances (dorsal = 0.002, ventral = 0.002, p = 0.73) so that one aspect did not capture more morphological variation. The outlines of frog skulls taken from photos and their corresponding diagrams were also not significantly different from one another (MANOVA, F = 1.11, p = 0.37) and morphological disparity was similar (photo = 0.0002, diagram = 0.0002, p = 0.70). Outlines of fossil reconstructions were compared to outlines of original specimens. Morphological disparity within each species is small (range 0.003–0.022), except for Notobatrachus degiustoi which had a Procrustes variance significantly larger than all other species (disparity = 0.037), which may reflect greater deformation in these specimens. These results, and a visualization of morphospace (Bardua et al. 2018, fig. S2), suggest that deformation does significantly influence the outline shape, with some specimens of the same species showing large morphospace distributions. Lastly, in a pairwise MANOVA taking Outline (D1/2/3) as a factor, the partial outlines (D2 and D3) were both significantly different to the complete outline (D1–D2 approx. F = 19.19, p = 2.22 × 10−15; D1–D3 approx. F = 2.74, p = 0.003). This analysis reveals the importance of tracing the posterior region as accurately as possible, which justifies the use of EFA instead of radial or tangent Fourier analysis since a greater sampling of points is allowed in regions of greater complexity. Results from the various sensitivity analyses therefore support the inclusion of diagrams, reconstructions and both dorsal and ventral aspect outlines in further analyses. However, caution must be taken for fossils whose outlines are hard to trace; some fossils were consequently rejected from the study (Bardua et al. 2018, table S2) and reconstructions were prioritized over photos.
The observed phylogenetic signal of the shape data (using PC scores) was significant for both the extant frog dataset (N = 87, 24 PCs, Kmult = 0.49, p = 0.001) and the combined dataset (N = 125, 25 PCs, Kmult = 0.48, p = 0.001). These results suggest that closely related species resemble one another, but phenotypically less than expected under a Brownian motion model.
Ecological niche was a significant discriminator of skull shape (MANOVA, 24 PCs, F = 1.43, p = 0.004) even after multiple‐test correction (p = 0.02) but not after phylogenetic correction (pMANOVA, R2 = 0.09, p = 0.37) (Bardua et al. 2018, fig. S3). Procrustes variances for all ecological niches were similar and not significantly different, ranging from 0.011 to 0.039. Frogs inhabiting the five categories of environment were not found to have significantly different skull outlines (Bardua et al. 2018, fig. S4) before or after phylogenetic correction (MANOVA, 22 PCs, F = 0.94, p = 0.62; pMANOVA, R2 = 0.07, p = 0.87) nor did they display different morphological disparities (Procrustes variance range 0.029 to 0.045, no significant differences). There was no significant difference in the outline shapes of the different zoogeographical realms, in either a MANOVA (23 PCs, F = 0.87, p = 0.84) or a pMANOVA (R2 = 0.05, p = 0.53). Disparity in zoogeographical realms in decreasing order was: Widespread (0.044), Neotropical (0.041), Nearctic (0.035), Australo‐Papuan (0.022), Afrotropical (0.019), Palaearctic (0.018) and Oriental (0.009). The Neotropical realm showed significantly higher disparity than the Australian, Palaearctic and Oriental realms (p = 0.007, 0.013, 0.038, respectively). The other differences in Procrustes variances were non‐significant. There was no significant difference in outline shape between the different levels of extinction risk (MANOVA, 23 PCs, F = 0.84, p = 0.83; pMANOVA, R2 = 0.06, p = 0.61) as defined by the IUCN Red List status. Finally, developmental strategy was a significant discriminator of frog skull shape (MANOVA, 23 PCs, F = 2.04, p = 0.016) even after correcting for multiple‐tests (p = 0.04) and phylogeny (R2 = 0.02, p = 0.039) although not after correcting for both (p = 0.20) (Bardua et al. 2018, fig. S5). Disparity was similar in metamorphosing and direct developing frogs (direct = 0.033; indirect = 0.001; p = 0.85).
Evolutionary allometry had a small but significant effect on the extant frog skull outlines (R2 = 0.08, p = 0.0004). Pooling extant and fossil frog data showed a smaller but still significant allometric effect (R2 = 0.05, p = 0.004). Visualizing morphospace distribution with size set as a factor revealed an overlapping distribution with a general trend of larger species having shorter, wider skulls and smaller species having narrower skulls (Bardua et al. 2018, fig. S6).
A MANOVA test illustrated a significant difference between taxa based on their age of first appearance (25 PCs, F = 1.59, p = 5 × 10−10) binned into 5 million year time bins, but not once phylogeny was incorporated (R2 = 0.26, p = 0.22). A visualization of the morphospace occupation (Fig. 5; and see Bardua et al. (2018, fig. S7) for a phylomorphospace) along the first two principal components reveals a distribution seemingly independent of species age. Triprion petasatus (extant) and Beelzebufo ampinga (fossil) have the most extreme skull outlines along PC1 and Shelania pascuali (fossil) and Rana basaltica (fossil) represent the extremes along the PC2 axis.
DiscussionThe anuran skull is highly reduced relative to the ancestral condition in terms of the size and number of cranial bones. This reduced cranial morphology is already seen in stem anurans (e.g. Notobatrachus degiustoi) after which the skull appears relatively unchanged. This qualitative assessment was confirmed quantitatively in our study, which found that morphospace occupation remained similar from the Middle–Late Jurassic until the present day, as measured by fossil and extant frog skull outlines (Fig. 5). The small change we observed in morphospace occupation is phylogenetically structured and probably reveals a preservational bias rather than truly changing morphospace, from the absence of some taxa in the fossil record (e.g. Centrolenidae) and abundance of others (e.g. Pipidae).
With a widespread distribution and adaptations to many diverse environments, frogs successfully occupy many ecological niches in zoogeographical realms across the globe. It has long been known that the Neotropical realm is the most diverse region for flora and fauna, boasting the highest diversity of frogs (Duellman 1999) with 2585 recorded species, 97% of which are endemic to this region (Bolaños et al. 2008). This is thought to be greatly underestimated (Fouquet et al. 2007) with a suggested 170–460 frog species unrecognized in this realm. Its consistent moist and warm climates with little seasonal temperature change are conducive to diversity. Although frogs have adapted to extreme environments such as brackish water (Lepidobatrachus asper) and deserts (Craugastor augusti), their requirements of moisture and warmth have meant that diversity is highest in areas satisfying these needs. Results have demonstrated that, as measured by skull outline, the Neotropical realm is the most disparate realm, indicating possible concordance of disparity and diversity in frogs. Convergence of body shape has been found in anuran species with similar ecological niches, including similar toe morphology in arboreal species and an association of body rotundity with fossoriality (Duellman & Trueb 1986; Wells 2007). However, ecological and environmental signals in our dataset appeared weak in extant anuran cranial outlines (see Bardua et al. 2018, figs S3, S4). Frogs from different environments showed no significant differences in skull outline, and the influence of ecological niche on skull outline was only significant before phylogenetic correction. Simplification and categorization of ecosystem data may dilute any ecological signals, as environments and ecologies exist on a continuous scale and are hard to discretize, but it is also possible that skull outline does not capture ecological differences among anurans. For example, degree of ossification is thought to be related to certain feeding and life habits (Trueb 1993, p. 332), but cranial outline may not capture differences in ossification. Incorporation of skull depth may therefore provide further insights into ecological influences on cranial morphology.
Developmental strategy is easier to discretize, and was found to have a weak effect on cranial morphology, but no effect on cranial morphological disparity (see Bardua et al. 2018, fig. S5). Tadpole morphology has been postulated to impose constraints on the resulting adult morphology, for example the loss of a free‐living larval stage may have influenced the high morphological diversity of plethodontid salamanders (Wake & Hanken 1996). This ontogenetic constraint hypothesis can be contrasted with the idea of antagonistic selection creating pressure for different ontogenetic life stages to evolve independently, allowing them to respond to different environments and requirements. A larval stage does not appear to impose constraints on the body morphology of Australian anurans, with a recent study finding adult morphology was shaped strongly by phylogenetic signal and tadpole morphology influenced more by homoplasy (Sherratt et al. 2017). This study (Sherratt et al. 2017) illustrates how developmental strategy is only one aspect of development; tadpole morphology can evolve independently from adult morphology in at least Australian anuran species. Incorporation of tadpole specimens in studies of cranial morphology may therefore aid the understanding of how development influences adult cranial shape. Our analyses of extant anurans found that cranial disparity was not significantly different between direct and indirect developers, further suggesting that anuran morphology is not restricted by metamorphosis, despite the prerequisite of a tadpole morphology. Developmental strategy did significantly discriminate frogs based on skull shape, but accounted for only a small amount of the shape variation, suggesting that life history does not greatly influence the adult skull shape. Our study therefore finds the existence of a free‐living larval stage to neither greatly facilitate nor constrain adult morphology, in terms of cranial outline.
Instead, our study found morphology to be more driven by phylogenetic relationships than ecological or developmental influences (see Bardua et al. 2018, figs S3–5). Ecological and developmental signals, when significant, were often erased by correcting for phylogenetic signal. This result suggests that closely related species share similar niches and lifestyles, despite resembling each other phenotypically less than expected under a Brownian motion model. Phylogenetic conservatism has been found previously in Australian anurans, in both skull morphology (Vidal‐García & Keogh 2017) and body shape and size (Vidal‐García et al. 2014), suggesting the prevalence of a phylogenetic signal in anuran morphology. Furthermore, a study of toad skulls has revealed that evolutionary constraints can restrict morphological responses to past climate change (Simon et al. 2016). Fossil amphibian studies have had similar findings: evidence has been found for morphological convergence within groups, and divergence between groups, both for Palaeozoic temnospondyls (Angielczyk & Ruta 2012) and Mesozoic stereospondyls (Stayton & Ruta 2006). Our results are therefore concordant with previous studies on extant and fossil amphibians, suggesting the prevalence of phylogenetic constraints on anuran evolution through deep time.
Allometry can be studied across individuals of a given age (static), across different life stages (ontogenetic) and across a phylogeny (evolutionary). Craniofacial evolutionary allometry has been studied across placentals and marsupials, with a general trend of small species exhibiting shorter faces (e.g. Cardini & Polly 2013; Cardini et al. 2015). This trend has, however, been found to be variably present across tetrapods, evidenced in temnospondyls (Angielczyk & Ruta 2012) in which small size is associated with a short, broad snout and large orbits in some clades but not others. Our study found allometry to be a significant influence on anuran cranial outlines through deep time, although the effect was small. Visualization of the morphospace (Bardua et al. 2018, fig. S6) showed that larger species appear to have relatively wider skulls, whilst smaller species had narrower skulls. This result suggests that the pattern of craniofacial evolutionary allometry across mammals may not be true outside amniotes, or indeed Mammalia. Our findings therefore suggest that size‐dependent selection pressures or constraints on cranial morphology may differ between anurans and mammals.
Beyond anuran evolution and biology, our results have implications for future studies of fossil anurans, as well as other taxa. We described the complex outline of extant and fossil frog skulls by a series of simple curves, and showed that whilst extant frogs required 14 harmonics to describe 99% of the shape, fossil frogs only required 12 (although we used 14 for both in the combined analysis). The additional two harmonics required to explain the extant frog skull outlines suggest that high order detail may have been lost from the fossils, since the fossil skull outlines can be described using more simple curves. It has been shown experimentally by studying the decay of Amphioxus and ammocoete larvae (Sansom et al. 2010) that higher detail is more vulnerable to loss of information through decay. The most phylogenetically useful information is more susceptible to loss in the fossil record whilst the plesiomorphic characters are more robust. This results in fossils being interpreted as more basal than their true phylogenetic position. It has been suggested that this stemward slippage is an important bias often overlooked in the fossil record (Sansom et al. 2010). In a similar fashion, higher frequency harmonics are more readily lost to noise than the lower frequency harmonics. Comparing harmonic requirements in extant and fossil datasets using EFA may therefore help to illustrate the missing information in fossils in future studies.
The potential caveats of an outline study on fossil specimens (the inclusion of a variety of data types and different levels of deformation) were addressed in this study through sensitivity analyses, but require further investigation. Sensitivity analyses supported the inclusion of fossils in both dorsal and ventral aspects as well as the use of diagrams and reconstructions. This may have implications for the study of similarly flattened fossils, although the amount of shape information captured in different two‐dimensional aspects will be strongly dependent on the fossil. A method has been developed which allows the comparison of two and three‐dimensional shapes (Cardini 2014) providing that both of these data can be collected from the same dataset, which is a useful avenue for further work along these lines. In addition, this study begins to assess the impact of variable deformation within a dataset, and future studies should further explore these effects where possible. Our study did find significant differences between fossil specimens of the same species (Bardua et al. 2018, fig. S2) but the similarity observed in morphospace occupation between extant and fossil frog skulls suggests that deformation of fossil frog skull outlines did not introduce excessive noise into our results.
ConclusionThe flattened nature of most fossil anuran skulls renders a three‐dimensional analysis incorporating most fossil taxa currently unachievable. Whilst the incorporation of skull depth information would better illustrate the morphospace occupation of fossil and extant frogs, this requires three‐dimensional reconstructions of flattened fossils, or the discovery of better preserved fossil anurans. Elliptical Fourier analysis was therefore used to quantify two‐dimensional extant and extinct frog skull shape. This study compared extant and fossil frog morphologies to evaluate phylogenetic, ecological and developmental signals. We found that phylogeny overrode most ecological or developmental signals, and the observed change in morphospace occupation through time appears to have been phylogenetically driven as well. This study aimed to highlight the possibility of including less‐well preserved specimens in studies of morphology. This increased inclusivity does sacrifice the robustness of results to some extent, but it means greater temporal breadth and the inclusion of many fossil specimens which would otherwise remain excluded from morphological studies, limiting the sampling of past diversity.
Author contributionsCB collected and analysed data, and wrote the first draft of the manuscript. All authors contributed to the project design, provided ideas and discussion, contributed to the writing of the manuscript, and read and approved the final version.
We would like to thank Ryan N. Felice, Thomas Halliday and Marcela Randau for help with analyses; Katherine P. Belessiotis, Sean Wilcox, Marc Jones and Ryoko Matsumoto for the collection of extant frog images; and our two reviewers, Liping Dong and Andrew McIntosh, for helpful feedback and comments. This research was funded through ERC grant STG‐2014‐637171 (to AG).