Do all frogs swim alike? The effect of ecological specialization on swimming kinematics in frogs

Frog locomotion has attracted wide scientific interest because of the unusual and highly derived morphology of these animals (Barclay, 1946; Estes and Reig, 1973; Zug, 1978; Frost et al., 2006). Frogs are characterized by a shortened trunk and tail, elongated ilia and elongated hind limbs. This morphology has been interpreted as being associated with a jumping life style and thus it has been suggested that jumping evolved early in the evolution of the lineage (Gans and Parsons, 1966; Shubin and Jenkins, 1995; Jenkins and Shubin, 1998) and many recent studies have attempted to infer locomotion in basal frogs (Prikryl et al., 2009; Essner et al., 2010; Reilly and Jorgensen, 2011; Sigurdsen et al., 2012; Venczel and Szentesi, 2012; Jorgensen and Reilly, 2013). However, kinematic and electromyographic studies indicate strong similarities between the mechanics of swimming and jumping in some frogs (Emerson and De Jongh, 1980; Peters et al., 1996; but see Nauwelaerts and Aerts, 2003), implying that morphological features associated with these two locomotor modes may not be that different. This may, in turn, complicate inferences of locomotor modes from anatomy as is often done for extinct animals. Despite their rather uniform morphology, frogs are an ecologically diverse and speciose group with over 5000 known species (Frost et al., 2006). Moreover, animals with different ecologies have evolved different morphologies and show different levels of locomotor performance (Moen et al., 2013), suggesting that locomotion differs in animals with different ecologies.

To date, most of our knowledge on frog locomotion is based on data for a limited set of derived frogs including ranoids [mostly ranids and bufonids (Calow and Alexander, 1973; Lutz and Rome, 1994; Kamel et al., 1996; Peters et al., 1996; Olson and Marsh, 1998; Gillis and Biewener, 2000; Nauwelaerts and Aerts, 2002; Nauwelaerts and Aerts, 2003; Nauwelaerts and Aerts, 2006; Nauwelaerts et al., 2001; Nauwelaerts et al., 2004; Nauwelaerts et al., 2005a; Nauwelaerts et al., 2005b; Johansson and Lauder, 2004; Stamhuis and Nauwelaerts, 2005)] and highly specialized aquatic pipids (Gal and Blake, 1988; Richards and Biewener, 2007; Richards, 2008; Clemente and Richards, 2013). A comparison of swimming kinematics between the highly specialized aquatic pipids and more generalized terrestrial species showed differences in joint kinematics, indicating differences in the underlying propulsive strategies of swimming across species (Richards, 2010). Although these data suggest that frogs with different ecologies differ in their limb kinematics, this remains to be tested using a broader sample of species with different ecologies and from different phylogenetic backgrounds. For example, the only study on swimming in primitive leiopelmatid frogs demonstrates an alternative swimming pattern consisting of an asymmetric swimming gait (Abourachid and Green, 1999) that may be related to the low locomotor speeds observed in these animals (Nauwelaerts and Aerts, 2002).

Here, we explore the diversity in hind limb kinematics during the propulsive phase of swimming by studying eight species of frogs from different families and with different ecologies (Table 1, Fig. 1). We include species with different ecologies (aquatic, terrestrial and semi-aquatic) and different phylogenetic affinities. Given the importance of pelvic girdle movements during locomotion in frogs (Emerson, 1976; Emerson, 1979; Videler and Jorna, 1985), we decided to use cineradiography rather than typical external high-speed video recordings to quantify swimming kinematics. Specifically, we test the hypothesis that species with different ecologies will differ in the kinematics of limb movement during swimming, with aquatic species showing greater velocities of movement and greater angular displacements at the distal-most joints associated with the rotation-powered swimming style observed in highly specialized swimmers (Richards, 2010).

Swimming in all species involved limb extension with significant movements at the hip, knee and ankle (Figs 2, 3 and 4; supplementary material Figs S1–S6). Whereas terrestrial and semi-aquatic species showed a clear proximo-distal extension sequence starting at the hip and ending at the ankle, this was not the case in specialized aquatic species where extension was initiated at the level of the knee, followed by the hip and the ankle. However, the greatest differences were observed in the movements at the distal-most segments (i.e. midfoot angles). Whereas in all species movements at the proximal foot were observed resulting in an extension of the foot fairly late in the kick, in specialized aquatics, the distal-most part of the foot (midfoot 2) was extended throughout the extension cycle. In the other species this angle showed little change and the distal foot remained extended throughout the extension cycle. Movements in the highly specialized aquatic species were also more stereotyped with lower variability, especially at the distal-most segments as suggested by the fact that they occupy only a small part of the kinematic space (Fig. 5).

A factor analysis performed on the mean kinematic variables per individual extracted four axes jointly explaining 79% of the overall variability in the data set (Table 2). Whereas the first axis (35.68%) was principally determined by extension of the limb, the velocity at the hip, knee and ankle as well as the total angular change at the hip and knee, the second axis (18.42%) was determined by the changes in midfoot angle as well as the minimal midfoot angles (Table 2). The third axis (14.2%) was determined by the pelvic shift and the change in pelvic angle (Table 2). The fourth axis (10.85%) was determined by the relative velocity of the animal (Table 2).

A multivariate analysis of variance performed on the raw kinematic variables that showed scores greater than 0.7 on the first two axes indicated a highly significant difference in swimming kinematics in animals with different ecologies, irrespective of the fact whether Bombina orientalis was classified as aquatic or semi-aquatic (B. orientalis aquatic: Wilks' lambda=0.095, F_20,26=2.92, P=0.006; B. orientalis semi-aquatic: Wilks' lambda=0.040, F_20,26=5.22, P<0.001). For the analysis with B. orientalis classified as aquatic, the subsequent univariate ANOVAs indicated that this difference was due to a significant effect on the delta knee angle, the delta hip angle and the delta and minima of the midfoot 1 and midfoot 2 angles (Table 3). Post hoc tests indicated that differences were significant between aquatic and terrestrial species in the delta hip angle with terrestrial species having a larger overall rotation at the hip. Moreover, differences were significant between the semi-aquatic species on the one hand and the aquatic and terrestrial species on the other hand for all midfoot angles, with semi-aquatic species having larger minimal angles yet smaller overall changes in angle. The only exception was for the minimal midfoot 2 angle, where aquatic and semi-aquatic species did not differ. For the analysis with B. orientalis classified as semi-aquatic, the univariate ANOVAs also indicated differences in the angular excursion at the hip and the midfoot (Table 3). Results of Bonferroni post hoc tests showed identical results to the analyses with B. orientalis classified as aquatic.

A plot of the phylogeny in the kinematic space constructed by using species means suggests that phylogeny is not driving the observed results (Fig. 6). For example, whereas the two terrestrial species Bufo calamita and Rhaebo guttatus are more closely related to Pelophylax esculenta they fall out with the terrestrial archeobatrachian Pelobates fuscus (Fig. 6). Thus the structuring in the kinematic space represents ecological affinities rather than representing phylogeny. The factor analysis performed on the species means shows a similar structuring to that observed using the individual data.

Our results show interesting differences in swimming behavior between species with different ecologies. Semi-aquatic species demonstrated a lack of changes in the midfoot angle during the extension phase, which was rather stably maintained. This is in contrast to specialized aquatic species such as Xenopus laevis and terrestrial species such as P. fuscus and R. guttatus where the midfoot actively contributes to generating propulsion. These results confirm previously published data on frog swimming (Richards, 2010), which demonstrated that the highly specialized aquatic X. laevis obtained nearly 100% of the total thrust during swimming through foot rotation involving tarso-metatarsal extension. Other species such as the semi-aquatic Rana pipiens or the terrestrial Bufo americanus had strong translational components to the kick. Interestingly, our analysis on species means suggested that terrestrial species have greater angular changes at the hip compared with aquatic and semi-aquatic species.

However, our results also show differences compared with previous studies. Notably, whereas Richards (Richards, 2010) found that foot rotation was greater in X. laevis compared with B. americanus, our results show that at least one of the bufonids (R. guttatus) shows greater foot rotation than X. laevis. The other species of bufonid included in our study (Bufo calamita), however, clustered with aquatic or semi-aquatic species depending on the classification of B. orientalis as aquatic or semi-aquatic (Fig. 5). Moreover, B. calamita also showed early knee extension, as has been observed for B. americanus, in contrast to the other terrestrial species in our data set (supplementary material Fig. S2). This suggests that differences in kinematic strategies exist within groups of closely related species with similar lifestyles. Further studies exploring swimming strategies in terrestrial bufonids would be especially insightful in this context.

In addition to confirming previous results (Richards, 2010), our results show significant differences between species with different ecologies. Indeed, our kinematic analysis showed that terrestrial species are significantly different from aquatic and semi-aquatic ones. Moreover, as our analysis included both primitive and derived species, this suggests that it is not a phylogenetic effect, but likely driven by the constraints of locomotion in different media. This is confirmed by the analysis on species means where a plot of the phylogeny in the kinematic space showed that structuring is largely according to ecological grouping rather than phylogeny (Fig. 6). Although differences between species and ecological groups were rather robust, the B. calamita included in our data set fell within the kinematic space of both highly specialized aquatic and semi-aquatic species, suggesting that interesting differences in locomotor strategies also exist within ecological groups. Moreover, our analysis on the individual means showed that one of the P. esculenta used in our analysis differed strongly from the other individual by showing much slower limb extension and a much lower contribution of the midfoot to overall propulsion. This result is hard to explain given the tight clustering around the species means of all other individuals used in the analyses. One possible explanation might be that this was a sub-adult individual, which could mean that that locomotor strategies vary throughout ontogeny; however, this remains to be tested.

Of the kinematic variables measured, those associated with the sliding of the pelvis did not contribute to the overall variation in kinematics of swimming. However, the highly specialized sliding pelvis of pipids has previously been suggested to play an important role during swimming by increasing the length of the power stroke (Palmer, 1960; Videler and Jorna, 1985). Despite the fact that two pipids were included in our data set, the average values of pelvic sliding were only slightly greater than those observed in other species that do not possess a highly specialized sacral joint which allows extensive sliding of the pelvis (supplementary material Fig. S6). Moreover, rather than lengthening, the distance between the tip of the ilium and the tip of the sacrum decreased, suggesting a forward sliding of the pelvis relative to the sacral joint during the extension phase of swimming. This suggests that the role of the pelvic joint needs to be re-evaluated and that its function may be related to escape behavior or even burrowing, as previously suggested (Whiting, 1961; Videler and Jorna, 1985).

Although previous studies found no trade-off between jumping and swimming kinematics or performance (Peters et al., 1996; Kamel et al., 1996; Nauwelaerts et al., 2007) our results suggest subtle but important differences in the kinematics of swimming that might be the result of specializations to different lifestyles. The principal differences observed are overall changes at the hip, which appear to characterize terrestrial species and differences in the kinematics of the distal limb elements, more specifically the foot, which appear to characterize semi-aquatic species. Whereas aquatic and terrestrial species appear to actively recruit the foot in generating propulsion, semi-aquatic species appear to have a relatively invariant foot angle throughout the limb extension cycle. This might be due to stiffer distal elements which could diminish the potential for the foot to contribute to the generation of propulsion, but this remains to be examined further. These results also suggest that locomotor inferences on extinct animals benefit from an examination of these distal elements rather than the often used proximal elements, such as the hip and proximal femur (e.g. Jorgensen and Reilly, 2013; Venczel and Szentesi, 2012). Unfortunately, such elements are rare in the fossil record, thus hampering our understanding of the evolution of locomotion near the base of the anuran tree.

Two B. orientalis, one B. calamita, two R. guttatus, four Discoglossus pictus, three P. fuscus, three Pipa pipa, two P. esculenta and seven X. laevis of undetermined sex, yet from phylogenetically different backgrounds (Fig. 1), were used in the recordings. Animals were housed individually in a temperature-controlled room and provided with food consisting of crickets, earthworms and waxworms twice weekly. For each individual, the snout–vent length, the length of the femur, the tibiofibula and the tarso-metatarsus were measured on X-ray images of anesthetized frogs (MS222) by digitizing the proximal and distal ends of each limb segment (Table 1). All experiments were approved by the ethics committee at the University of Antwerp, Belgium.

Animals were recorded in dorso-ventral view while swimming using a Phillips Optimus X-ray unit with a 14 inch image intensifier and coupled to a Redlake Imaging MotionPro 2000 high-resolution digital video camera set at a recording frequency ranging from 250 frames s⁻¹. Swimming was recorded in an experimental tank of 120×25×50 cm with 10 cm of water restricting swimming to a single horizontal plane parallel to the image intensifier. Test temperature varied between 20 and 24°C for all swimming trials. Swimming was elicited by tapping the animal at the base of the urostyle with a long, thin metal rod. In all cases the stimulus was provided by the same person and such that the frog was unaware of the rod before it touched the animal. For smaller species or species that showed a poor degree of ossification (D. pictus, B. calamita), small radio-opaque markers were implanted at the different limb joints of interest to facilitate analysis of the kinematic data. Markers were implanted percutaneously using hypodermic needles under full anesthesia with MS222.

Five swimming sequences were recorded for each individual and those where the frog stayed in the plane parallel to the image intensifier were retained for analysis. This resulted in nine sequences for two individuals of B. orientalis, 17 sequences for four individuals of D. pictus, 26 sequences for seven individuals of X. laevis, eight sequences for three P. pipa, 16 sequences for four P. fuscus, five sequences for one B. calamita, eight sequences for two R. guttatus and nine sequences for two P. esculenta for a total of 98 analyzed sequences.

On each frame, 21 landmarks were digitized for the limb extension cycle using Didge (version 2.2.0.; A. Cullum) (Fig. 7) and the X- and Y-coordinates for each point were exported to a spreadsheet. Landmarks used were (numbers indicated for one side only; see Fig. 7): the tip of the snout (1), the center of the sacrum (2), the distal end of the ischium (3), the left and right iliosacral joints (4), the left and right proximal head of the femur (5), the left and right distal end of the femur (6), the left and right proximal end of the tibiofibula (7), the left and right distal end of the tibiofibula (8), the left and right proximal end of the proximal tarsals (9), the left and right distal end of the tarsal bones (10), the left and right distal end of the longest metatarsal (11) and the left and right distal end of the terminal phalanx of the longest toe (12). Next, coordinates were re-calculated to a frame of reference moving with the frog and with the x-axis parallel to the midline of the frog and the y-axis going through the sacrum, thus making landmark 2 the origin of our new reference frame.

Based on the X- and Y-coordinates of these landmarks, the following kinematic variables were calculated: the pelvic angle, being the angle subtended by the lines interconnecting landmarks 1 and 3, and 2 and 3, respectively; the hip angle, being the angle subtended by the lines interconnecting landmarks 1 and 3, and 5 and 6, respectively; the knee angle, being the angle subtended by the lines interconnecting landmarks 5 and 6, and 7 and 8, respectively; the ankle angle being the angle subtended by the lines interconnecting landmarks 7 and 8, and 9 and 10, respectively; the midfoot 1 angle being the angle subtended by the lines interconnecting landmarks 9 and 10, and 10 and 11, respectively; the midfoot 2 angle being the angle subtended by the lines interconnecting landmarks 10 and 11, and 11 and 12, respectively. The hip angle, knee angle, ankle angle and both midfoot angles were calculated for both limb pairs. Additionally, the amount of pelvic sliding was calculated as the difference in the X-coordinates between markers 2 and 4. Finally, limb extension was calculated as the difference in the X-coordinates between marker 2 and a virtual tibio-tarsal joint marker calculated as the average between the X-coordinates of markers 8 and 9, respectively.

The displacements of all limb segments were plotted against time and smoothed using a zero phase shift fourth-order low-pass Butterworth filter with user defined cut-off frequency that was set iteratively to obtain smooth acceleration profiles without losing information in the displacement and velocity profiles (Winter, 2004). Next, the limb extension cycle was interpolated over 50 time points, allowing us to compare cycles across individuals and species. After interpolation, the velocity and acceleration of displacements and angular changes were calculated based on numerical differentiation of the displacement profiles. For statistical analysis, the peak snout velocity, peak snout acceleration, average velocity, the amount of pelvic sliding, the total limb extension, the peak limb extension and retraction velocity, the delta pelvic angle, the delta hip angle, the minimal hip angle, the minimal and maximal angular velocity at the hip (i.e. associated with hip extension and flexion respectively), the delta knee angle, the minimal knee angle, the minimal and maximal angular velocity at the knee, the delta ankle angle, the minimal ankle angle, the minimal and maximal angular velocity at the ankle, the delta and minimal midfoot 1 angles and the delta and minimal midfoot 2 angles were extracted (Table 4). As limb movements are not always perfectly symmetrical, the largest angular displacement and velocity of the right and left side was retained for further analysis.

Species were classified into three broad ecological groups on the basis of literature data. We considered as aquatic species, species that spend most of their time in water outside of the breeding season. As such, X. laevis, P. pipa and B. orientalis (Kaplan, 1992; du Preez and Carruthers, 2009; Ouboter and Jairam, 2012) were all classified as aquatic. As terrestrial species, we considered species that spend most of their time away from water outside the breeding season. These species thus cannot be found in the immediate vicinity of water outside the breeding season and include B. calamita, P. fuscus and R. guttatus (Arnold and Ovenden, 1978; Ouboter and Jairam, 2012). Finally, we classified as semi-aquatic, species that live near water outside the breeding season, yet typically jump into the water as an escape response. These species included D. pictus and P. esculenta (Arnold and Ovenden, 1978). However, given conflicting statements in the literature concerning B. orientalis, we ran all our analysis with this species classified both as aquatic and as semi-aquatic.

Next, all raw kinematic variables were averaged per individual. All variables were log₁₀-transformed and used as input for regression analysis with SVL as the independent variable. Where significant, residuals were extracted and saved as variables. Next, kinematic data (residual for those variables dependent on overall size) were used as input for a factor analysis with varimax rotation. Factors with eigenvalues over 1 were extracted and factor scores were saved. Factor scores were used to explore how species were distributed in kinematic space and to select kinematic variables for subsequent analysis. We selected all variables with loadings higher than 0.7 on the first two axes as input for a multivariate analysis of variance coupled to subsequent univariate ANOVAs. Bonferroni post hoc tests were then used to examine which groups differed from one another for each variable that showed significant effects.

Because species cannot be considered as independent data points or disconnected from their evolutionary history, comparative analysis has been advocated to take into account shared ancestry in explaining patterns of phenotypic or functional diversity. However, such approaches typically require a minimum number of species for the analysis to be robust. Given the time-consuming nature of kinematic analyses, our data set remains restricted. Thus, rather than doing explicit comparative analyses we decided to map the phylogeny onto the functional space, allowing us to evaluate whether or not structuring is driven by phylogeny. We did so using the phylomorphospace function in R (R Development Core Team, 2011) implemented in the ‘phytools’ library (Revell, 2012). We use two alternative phylogenies that differ in the placement of the basal most taxa (Pipoidea versus Bombinatoroidea) based on the phylogenies provided by Frost et al. (Frost et al., 2006) and Zhang et al. (Zhang et al., 2013) and pared down to include only the taxa in our analyses (Fig. 1). Moreover, we classified B. orientalis both as aquatic and as semi-aquatic. Branch lengths were computed using the Grafen method (Grafen, 1989) with the ‘compute.brlen’ function of the ‘Ape’ library (Paradis et al., 2004) in R (R Development Core Team, 2011).

We would like to thank Jan Scholliers for taking care of the animals in the lab.