Terrestrial force production by the limbs of a semi-aquatic salamander provides insight into the evolution of terrestrial locomotor mechanics

The evolutionary invasion of land was a seminal event in vertebrate history that has received intense study, with integrative approaches filling major gaps in our understanding of how this transition transpired (Ashley-Ross et al., 2013; Molnar et al., 2021). Skeletal remains and trackways in the fossil record can provide some of the most direct evidence about morphological transformations and gait patterns, respectively, as vertebrates became terrestrial (Ahlberg, 2018), yet also present general challenges given the limits to which fossils preserve biomechanical parameters, such as muscle leverage (Molnar et al., 2021) and bone loading (Hohn-Schulte et al., 2013), that are often simulated or validated based on empirical data from living taxa. Consequently, extant amphibious fishes, amphibians and reptiles have been used as functional models to infer the biology of extinct tetrapodomorphs (tetrapods and tetrapod-like sarcopterygian fishes) (Nyakatura et al., 2014; Nyakatura et al., 2019; Pierce et al., 2013), with extant taxa representing alternative strategies for invading land and potentially serving as models for different functional stages as vertebrates became terrestrial. Investigations of extant taxa exhibiting morphological and/or behavioral traits that are consistent with those of extinct tetrapodomorphs offer particularly intriguing opportunities to evaluate how tetrapods might have left the water's edge (Molnar et al., 2021; Pierce et al., 2013).

In considerations of locomotor performance during the invasion of land, salamanders are often used as functional analogs for early tetrapods as they move between water and land (Karakasiliotis et al., 2012), and exhibit a relatively generalized tetrapod Bauplan that has not changed substantially for at least 150 million years (Gao and Shubin, 2001) and possibly as long as 230 million years (Schoch et al., 2020). Living salamanders have been used to study the in vivo biomechanics and muscle physiology of walking underwater (Ashley-Ross et al., 2009; Azizi and Horton, 2004; Deban and Schilling, 2009; Frolich and Biewener, 1992), walking on land (Ashley-Ross et al., 2009; Brand, 1996; Deban and Schilling, 2009; Delvolvé et al., 1997; Frolich and Biewener, 1992; Kawano and Blob, 2013; Pierce et al., 2020; Sheffield and Blob, 2011) and transitioning between water and land (Ashley-Ross and Bechtel, 2004). However, the extant salamanders that have traditionally been studied in this context [e.g. Ambystoma (Sheffield and Blob, 2011), Salamandra (Pierce et al., 2020) and Dicamptodon (Ashley-Ross, 1994)] have typically been highly terrestrial (following ecological classifications from Fabre et al., 2020). For instance, Ambystoma tigrinum are found in various terrestrial habitats, ranging from conifer forests to deserts, and adults rarely spend extended durations in water for reasons other than reproduction (Petranka, 1998). Consequently, existing literature on the locomotor biomechanics of salamanders may reflect a biased sample that does not consider the continuum from aquatic to terrestrial habitats.

Given the greater effect of gravitational loads on the musculoskeletal system in terrestrial environments (Young et al., 2017), one of the most fundamental requirements for moving on land is the ability to support body weight for maintaining posture and generating the physical forces for locomotion. Locomotor kinetics of individual limbs during terrestrial locomotion have been assessed through measurements of ground reaction forces (GRFs) in the tiger salamander, A. tigrinum, and demonstrated that the hindlimbs support a similar proportion of body weight to the forelimbs but have a greater contribution to acceleration (Kawano and Blob, 2013). Moreover, force plate analyses demonstrated that tiger salamanders and tuataras use inverted pendulum mechanics to recover mechanical energy during walking and running (Reilly et al., 2006), which would benefit terrestrial invaders by reducing the costs of locomotion on land. However, it is unclear whether these patterns are consistent across salamanders that spend less time on land as adults. Additional studies that compare salamanders with different terrestrial capabilities could provide insight into whether terrestrial species exhibit traits that would improve their locomotor biomechanics on land.

Because salamanders have a diverse range of habitat preferences and life histories (Wake, 2009), they provide an opportunity to compare the locomotor biomechanics of species that spend different proportions of time in aquatic versus terrestrial environments. The Spanish ribbed newt, Pleurodeles waltl Michahelles 1830, is one of the largest salamanders that spends a greater proportion of its life cycle living in water (Obst et al., 1988) compared with more commonly studied (terrestrial) salamanders, yet still normally makes terrestrial excursions (Fig. 1) (Karakasiliotis et al., 2012). Previous studies on the muscle activity (Delvolvé et al., 1997), bone microanatomy (Canoville and Laurin, 2009; Laurin et al., 2004) and kinematics and morphology (Karakasiliotis et al., 2012) of P. waltl provide a foundation for comparisons with data from more commonly studied terrestrial ambystomatid salamanders (Ashley-Ross and Barker, 2002; Bennett et al., 1989; Deban and Schilling, 2009; Kawano and Blob, 2013; Laurin et al., 2004; Stokely and Holle, 1954). Thus, GRF data from the semi-aquatic P. waltl and more terrestrial taxa (e.g. A. tigrinum) can be used to evaluate whether limb function differs between salamanders that have different propensities to move on land.

Propulsion driven primarily by the hindlimbs (‘rear-wheel drive’) occurs during terrestrial locomotion in salamanders (Kawano and Blob, 2013), and may have originally evolved for aquatic locomotion when vertebrates were still fully aquatic (King et al., 2011). However, terrestrial propulsion in early stem tetrapods may have initially been dominated by the forelimb (‘front-wheel drive’), with the transition to hindlimb dominance (‘rear-wheel drive’) only occurring later in the water–land transition as the hindlimbs assumed a more important locomotor role (Boisvert, 2005; Molnar et al., 2020). The ‘front-wheel drive’ strategy proposed for the sarcopterygian fish Panderichthys (Boisvert, 2005) and the early tetrapod Ichthyostega (Pierce et al., 2012) might be appropriately modeled by extant amphibious fishes (i.e. Clarias catfish and mudskippers, respectively) whose pectoral fins are used for terrestrial excursions while the pelvic fins are not (Kawano and Blob, 2013; Pace and Gibb, 2014). In contrast, the ‘rear-wheel drive’ of more crownward early tetrapods might be better modeled using terrestrial salamanders (Pierce et al., 2013), such as A. tigrinum. Our past work showed distinct GRF patterns between these groups, with those produced by the pectoral fins of mudskippers directed more medially (∼17 deg) than those produced by the limbs of terrestrial salamanders, or almost any other terrestrial tetrapod (generally less than 11 deg) (Kawano and Blob, 2013). Such differences could have substantial implications for the stresses experienced by the limb bones by increasing moment arms of the GRFs that contribute to bending, which could make bones more vulnerable to fracture. However, it is unclear whether a salamander with lower terrestrial capabilities would have GRF patterns that are intermediate between those of the mudskipper fish and terrestrial salamanders, and in what ways the mechanisms that generate the GRF patterns differ across appendages.

In this study, we compared the kinetics of terrestrial locomotion in amphibious fish and salamanders to gain a better understanding of how locomotor biomechanics differ between fins and limbs. First, we compared GRFs of individual limbs produced by semi-aquatic salamanders (Spanish ribbed newts, P. waltl) with published data on terrestrial tiger salamanders, Ambystoma tigrinum Green 1825, and semi-aquatic African mudskippers, Periophthalmus barbarus (Linnaeus 1766) (Kawano and Blob, 2013). GRF data provide information about propulsion, maneuverability and stability, and represent the external forces that are loaded upon the bones as a result of making contact with the substrate. While peak force magnitudes have been commonly used to model peak stresses to the limb bones during terrestrial locomotion (Blob, 2001), rates of loading help trigger bone remodeling (Aiello et al., 2015). Thus, we also quantified the dynamics of external force development through calculation of ‘yank’, which represents the first time derivative of force (sensuLin et al., 2019). When evaluated at the organismal level, yank is measured as the rate of change in GRF. Positive and negative yank values would indicate that the rate of force production was increasing or decreasing, respectively. Yank has important implications for musculoskeletal function but was only recently defined formally, providing a promising avenue for future research on locomotor biomechanics.

Lin and colleagues (2019) presented case studies to highlight how yank measures can help explain dynamics in biological motor and sensorimotor feedback systems. Extending these analyses to a wider range of taxa (e.g. fishes and salamanders) may lend new perspectives on broadscale patterns of dynamic locomotor forces within and across taxa with diverse musculoskeletal designs. Yank was demonstrated to be an important measure for rapid movements such as jumping, where maximum take-off velocity is related to locomotor performance. Using a computer simulation, Lin and colleagues (2019) found that increasing yank decreased the time to take-off and increased jump velocity. Jumping animals can exhibit a trade-off between time-limited, submaximal force production (i.e. bullfrogs) versus non-time-limited, maximal muscle contraction (i.e. grasshoppers) (Rosario et al., 2016). Findings from the computational model led Lin and colleagues (2019) to propose yank as a mechanism to avoid the trade-off between reduced contact time versus faster take-off velocities. Yank can also be relevant to movements that are less time limited, such as sensory responses to environmental stimuli. For instance, Golgi tendon organs (GTOs) are distributed across the muscle–tendon junctions of vertebrates, respond to changes in muscle forces, and are primarily found in the skeletal muscles of tachymetabolic tetrapods (birds and mammals) versus in the tendons of bradymetabolic tetrapods (amphibians, lizards, turtles, crocodilians) (Granatosky et al., 2020). Differences in the functional morphology of GTOs allow tachymetabolic tetrapods to have more fine-scale resolution of internal force changes within muscles, compared with bradymetabolic tetrapods that would only detect coarse-scale force changes across the entire muscle, and more predictable loading regimes upon the hindlimbs that would reduce the metabolic costs of locomotion (Granatosky et al., 2020). Moreover, cursorial mammals, salamanders and tuataras exhibit comparable pendular mechanics but the last two taxa tend to exhibit larger changes in vertical displacement (i.e. potential energy) compared with mammals, which may be a plesiomorphic condition for tetrapod locomotion (Reilly et al., 2006). As a result, non-tachymetabolic taxa may exhibit larger yank magnitudes that reflect their variable locomotor mechanics and loading. Consequently, expanding yank analyses to salamanders and fishes may reveal how the dynamics of motor systems vary in taxa with more variable bone loading than mammals and birds, and shed light on the evolution of terrestrial locomotion in vertebrates.

Our data from semi-aquatic P. waltl can give new insights into the functional differences between fins and limbs during terrestrial locomotion in taxa that exhibit a range of terrestrial capabilities. Simply by having limbs, P. waltl locomotor force production may be similar to that of terrestrial salamanders such as adult A. tigrinum. However, the habitual use of the limbs for aquatic locomotion by adult P. waltl might lead to kinetic similarities with other semi-aquatic taxa, such as mudskipper fish. These comparisons carry broader implications for generating hypotheses about the emergence of functional disparity between locomotor structures, such as whether changes in functional performance are coupled across major structural changes (e.g. transition from fin to limb) or through gradual steps related to loading regimes that are potentially decoupled from structural changes (e.g. the transition from aquatic to terrestrial habitats, regardless of which locomotor structure is used).

Five adult P. waltl (body mass: 16.60±0.40 g; snout–vent length: 0.083±0.001 m; total length: 0.186±0.003 m) were obtained from a commercial vendor. All values represent means±1 s.e.m. Animals were individually housed in glass aquaria that were aerated with sponge filters, kept on a 12 h:12 h light:dark cycle, and fed every 1–2 days on a diet of defrosted bloodworms and krill. Prior to testing, animals were starved for 2 days to reduce the effects of satiation on locomotor performance (Sass and Motta, 2002). Experimental and animal care procedures were approved by the Institutional Animal Care and Use Committee at Clemson University (protocols 2009-071, 2010-066).

Experimental procedures from a previous study on the GRFs of tiger salamanders and mudskipper fish (Kawano and Blob, 2013) were replicated in the present study (see ‘Criteria for selecting trials for GRF measurements’ and ‘Experimental procedures’ in Supplementary Materials and Methods) to obtain forelimb (n=50) and hindlimb (n=49) GRFs from P. waltl (Fig. 1). Custom code to replicate these analyses is available from GitHub (https://github.com/MorphoFun/kraken). The focal taxa examined in the present study represent models for distinct functional stages during the evolution of terrestrial locomotion: front-wheel drive in a semi-aquatic vertebrate (semi-aquatic mudskipper P. barbarus fish), a semi-aquatic vertebrate with a generalized tetrapod Bauplan that has yet to be confirmed as being front-wheel or rear-wheel driven (semi-aquatic P. waltl salamander), and rear-wheel drive in a terrestrial vertebrate with a generalized tetrapod Bauplan (terrestrial A. tigrinum salamander). Baseline data collected from the force plate were zeroed at the beginning of each trial, when no structures had recently contacted the force plate. GRFs in the vertical, mediolateral and anteroposterior directions were digitally filtered with a custom low-pass, zero-phase second-order Butterworth filter, and then interpolated to 101 points (0–100% of stance at 1% increments) using a cubic spline with the signal R package (https://CRAN.R-project.org/package=signal). The stance phase was defined to begin at the first video frame when the entire foot was flat against the ground, and to end at the frame penultimate to the beginning of the swing phase.

Interspecific and intraspecific comparisons were made between fins and limbs. First, GRFs from the forelimbs and hindlimbs of P. waltl were compared with published data from the pectoral and pelvic appendages of mudskipper fish and tiger salamanders (Kawano and Blob, 2013) to assess whether the kinetics of semi-aquatic salamander limbs were more similar to those of mudskipper pectoral fins or terrestrial salamander limbs. Second, GRFs were compared between the forelimbs and hindlimbs for P. waltl to determine whether this species used forelimb-driven or hindlimb-driven propulsion. Third, GRF data from these three species were used to conduct new analyses of yank (sensuLin et al., 2019) to quantify the rate of change in the GRFs produced by fins versus limbs. Linear mixed effects models (LMMs) were used to compare magnitudes and angles of orientation of GRFs when forces were maximal (‘peak net GRF’), and the range of yank values between the species and appendages (see ‘Statistical analyses’, below).

LMMs with random intercepts via the lme4 package (Bates et al., 2015) were used to compare the biomechanical variables of interest while accounting for the non-independence of multiple trials sampled per individual, and fitted with restricted maximum likelihood (REML) to produce unbiased estimates. All analyses were conducted in R version 3.6.1 (http://www.R-project.org/) and RStudio version 1.2.1335 (http://www.rstudio.com/). LMMs were conducted to account for the non-independence caused by trials being nested within individuals and the unequal sample sizes across individuals (for additional information, see ‘Details on statistical analyses’ in Supplementary Materials and Methods). LMMs also have the advantage that they do not make assumptions about balanced designs, homogeneity of variance, etc., so the model assumptions only include random sampling and normality of the residuals (reviewed in Smith, 2017). Evaluation of quantile–quantile (Q–Q) plots indicated that there were no major deviations from the reference line (see Figs S1–S3), suggesting that all variables reasonably met the assumption of normality. To assess whether any outliers affected our results, we re-ran our LMMs with the outliers removed and found that the point estimates were not drastically affected (Tables S1 and S2), so we present the results from our full dataset. However, violations of the model assumptions would not have affected our analyses given our focus on fixed effects, which tend to be unbiased even if assumptions are violated (LeBeau et al., 2018). ‘Individual’ was treated as a random effect, and ‘appendage type’ (pectoral versus pelvic) or ‘species’ was treated as a fixed effect for intraspecific and interspecific comparisons, respectively.

Large sample sizes are typically required to produce accurate estimates of a random slopes model (reviewed in Harrison et al., 2018), so we assumed that the slopes were similar across the individuals within a species based on our smaller sample size. Assuming common slopes in a LMM can increase the potential for false positives and false negatives when calculating P-values for null hypothesis testing (reviewed in Harrison et al., 2018); however, the present study focused on effect sizes (i.e. the magnitude of the effects or variables) rather than the P-values, so these outcomes do not detract from our analyses. Given that one of the objectives of this study was to quantify the mean differences (i.e. point estimates) and range of variation (i.e. interval estimates) in GRF production and yank between and within salamanders, our emphasis on calculating effect sizes is appropriate (Cohen, 1990) because of the limited information that often results from traditional methods of null hypothesis testing (Cumming, 2014).

To calculate the point estimates (i.e. mean values) and interval estimates (i.e. standard errors) of the fixed effects, we used the emmeans::emmeans() function (https://CRAN.R-project.org/package=emmeans). The coefficient of determination (percentage of the variation explained by each model) was calculated as Nagakawa's marginal R² (R²_LMM(m); variation due to the fixed effects only) and conditional R² (R²_LMM(c); variation due to the fixed and random effects) using the performance::r2_nagakawa() function (Lüdecke et al., 2019). Dynamic changes in GRF and yank were plotted by pooling individual trials for a given species, but results from the LMMs are reported in the tables to account for repeated trials within individuals.

We then ran additional statistical analyses to determine whether our kinetic variables were influenced by stance duration. Stance duration was relatively similar between the two salamanders for the forelimbs (P. waltl: 0.545±0.027 s; A. tigrinum: 0.526±0.026 s) and hindlimbs (P. waltl: 0.607±0.036 s; A. tigrinum: 0.596±0.036 s) but was lower for P. barbarus pectoral fins (0.388±0.026 s), so we created additional LMMs with the interaction between species and stance duration as a fixed effect for the comparisons between pectoral appendages. To determine whether including stance as an interaction term improved our models for the pectoral appendages, we used likelihood ratio tests (LRTs) between the reduced model (LMM with species as a fixed effect and individual as a random effect) and the full model (LMM with a species×stance duration interaction term as a fixed effect and individual as a random effect) using the stats::anova() function in R. Stance duration did not improve the model for the net GRF (LRT: chi-squared: 7.979, d.f.=4, P=0.092) and the vertical component of the GRF (LRT: chi-squared: 2.417, d.f.=3, P=0.491). Including stance as an interaction with species did improve the model for the mediolateral component of the GRF (LRT: chi-squared=20.323, d.f.=4, P=0.0004), the anteroposterior component of the GRF (LRT: chi-squared=41.086, d.f.=4, P<0.0001), the mediolateral angle of the GRF (LRT: chi-squared=19.818, d.f.=4, P=0.0005) and the anteroposterior angle of the GRF (LRT: chi-squared=48.640, d.f.=4, P<0.0001). However additional analyses (see Table S3 and ‘Comparison of stance durations’ of the Supplementary Materials and Methods) indicated that stance duration had minimal effects on the point estimates from our LMMs and that these differences had no effect on the main patterns of our results. As such, we report the results from our reduced models to provide more straightforward interpretations of the statistical analyses.

The forelimbs of P. waltl exhibited kinetics that were distinct from (and sometimes intermediate between) those of the pectoral fins of P. barbarus fish and the forelimbs of A. tigrinum salamanders (Fig. 2). Force production during terrestrial locomotion generally approximated a bell-shaped curve that was skewed to the right in P. barbarus and A. tigrinum, but relatively broad and shallow in P. waltl for the net GRF and the vertical component of the GRF (Fig. 2A,B). The mediolateral component of the GRF was directed medially during stance for all three species, yet had a wider range of values in P. barbarus (Fig. 2C). For the anteroposterior component of the GRF, forelimb kinetics in P. waltl were more similar to those of terrestrial A. tigrinum salamanders than to semi-aquatic P. barbarus fish. The pectoral fin is the primary propulsor during crutching for P. barbarus and its GRF was oriented anteriorly for all of stance (indicating an acceleratory role), whereas the forelimbs of A. tigrinum and P. waltl showed primarily posteriorly directed GRFs (deceleratory role), with only one or two periods in which the GRFs became slightly anterior (Fig. 2D,F). The GRF was angled medially in the horizontal plane throughout stance for P. barbarus and A. tigrinum, and for almost all of stance for P. waltl except for a slight lateral orientation starting around 90% of stance (Fig. 2E). These findings illustrate how different aspects of force production by the forelimbs of P. waltl share similarities with those of semi-aquatic fins and terrestrial limbs, yet also show distinct differences (e.g. mediolateral GRFs).

GRF parameters at the peak net GRF also differed between the pectoral appendages for the three groups examined. The timing of the peak net GRF occurred at 60.4±2.86% of stance for A. tigrinum, 56.0±2.86% for P. barbarus and 46.9±2.93% for P. waltl (mean±s.e.m.; R²_LMM(m)=0.226, R²_LMM(c)=0.462). The magnitude of the peak net GRF was comparable between P. barbarus and A. tigrinum at approximately 44% and 46% of body weight, respectively, but was roughly 40% of body weight and had a wider range of values in P. waltl (Table 1, Fig. 3A). Similarly, the vertical component at the peak net GRF was comparable for A. tigrinum and P. barbarus yet lower and broader for P. waltl (Fig. 3B). The mediolateral component at the peak net GRF was similar between the two salamanders (Fig. 3C), with values for P. barbarus being 1.5–2 times more medial than those of the salamanders. The anteroposterior component was also more similar between the salamanders at the peak net GRF, with values being anterior (acceleratory) for P. barbarus but posterior (deceleratory) for both salamanders (Fig. 3D). The angle of the GRF was medial for all pectoral appendages at peak net GRF, but mean values for P. waltl were intermediate (−14 deg) between those of P. barbarus (−17 deg) and A. tigrinum (−9 deg) (Fig. 3E). These comparisons further demonstrate that kinetics at the time of peak net GRF for P. waltl share similarities with those of semi-aquatic fins and terrestrial limbs, with some features intermediate between the other two species.

Kinetic profiles were similar in shape between the hindlimbs of semi-aquatic and terrestrial salamanders (Fig. 4), but there were differences in the magnitudes of the curves. Semi-aquatic and terrestrial hindlimbs had a peak net GRF occurring around 30% of stance; however, the net GRF and the vertical component at peak GRF for P. waltl hindlimbs [0.349 and 0.314 body weight (BW), respectively] were more similar to those of P. waltl forelimbs (at about 30% of body weight) than to those of the hindlimbs of A. tigrinum (∼50% of body weight) (Table 1, Fig. 5A,B). The GRFs were directed medially for all limbs (−9 to −17 deg) but were 1.5 and 2 times more medial in P. waltl than in A. tigrinum for the forelimbs and hindlimbs, respectively (Table 1, Fig. 5E). The forelimbs of both salamander species had a net deceleratory role (Table 1), whereas the hindlimbs had a net acceleratory role (Table 1), with the hindlimbs of A. tigrinum showing a mean value ∼1.5 times higher than those of P. waltl at the time of peak net GRF (Fig. 5D). However, the angle of the GRF in the anteroposterior direction was comparable between salamander hindlimbs, inclined 16–17 deg anteriorly (Fig. 5F). These findings demonstrate that hindlimb kinetics are not necessarily equivalent between salamander species with differing degrees of terrestriality, and that locomotor biomechanics differ between the forelimbs and hindlimbs during terrestrial walking in both species.

Average yank values for the pectoral appendages tended to be similar between the three species, but GRFs produced by the forelimbs of P. waltl generally reached their peak more slowly (as indicated by smaller absolute magnitudes of yank) and had a narrower range compared with those of the semi-aquatic fish and terrestrial salamanders (Table 2). The shapes of the curves for yank differed between P. waltl and the other two species. For P. barbarus and A. tigrinum, yank for the net GRF and vertical component of the GRF (Fig. 6A,B) approximated sine waves that started close to zero at the beginning of stance, increased to positive yank values before decreasing to an inflection point around 60% of stance, and then decreased to negative yank values before returning to the zero baseline. In contrast, the net GRF and vertical component of the GRF for P. waltl reached an inflection point at ∼50% of stance, and then approached zero again at ∼80% of stance, at a time when A. tigrinum and P. barbarus were both reaching their minimum yank values (i.e. peak in negative yank). The anteroposterior component was comparable between the two salamander forelimbs (∼0.006–0.008 BW for each percentage of stance), which increased faster than that for the pectoral fin of P. barbarus (0.004±0.001 BW for each percentage of stance). These findings demonstrate that rates of GRF production in the pectoral appendages of the semi-aquatic salamander, P. waltl, were generally comparable in magnitude but shifted at different times during stance and encompassed a narrower range of values than for the other two species.

Similarities between A. tigrinum and P. waltl were more evident in the shapes of the yank curves for the hindlimbs (Fig. 7). The curves for the net GRF (Fig. 7A) and the vertical component of the GRF (Fig. 7B) approximated sine waves that reached an inflection point at around 25% of stance. The peak for positive yank appears almost 2 times higher in A. tigrinum than in P. waltl (Fig. 7, Table 2), but the confidence intervals overlapped for the fixed effects. The minima for negative yank in the vertical component and net GRF were shallower in P. waltl than in A. tigrinum. As a result, rates of change in GRF production for the hindlimb were slower for the semi-aquatic P. waltl salamander compared with those in the terrestrial A. tigrinum salamander for the vertical component of the GRF and net GRF but profiles for yank in the mediolateral and anteroposterior directions were similar between the two taxa.

Pleurodeles waltl is a semi-aquatic salamander that spends most of its life history underwater. Similar to terrestrial salamanders, the predominant acceleratory forces in this semi-aquatic salamander are produced by the hindlimb, signifying rear-wheel drive. The orientation of the GRF in the hindlimb was about 17 deg in the anterior direction for both salamanders, but more medial in P. waltl (−17 deg) than in Ambystoma tigrinum (−9 deg). Similarly, the orientation of the GRF for the pectoral appendages was more medial in the semi-aquatic P. waltl and P. barbarus than in the primarily terrestrial A. tigrinum; however, the anteroposterior GRF of both salamanders indicated a net deceleratory force compared with the acceleratory role in P. barbarus (Table 1). These differences in orientation may affect how loads are applied to the appendages, potentially resulting in the lower vertical component and net magnitude of the GRF in P. waltl forelimbs and hindlimbs.

Although the forelimbs of both salamanders had a net deceleratory role as evidenced by the negative value for the anteroposterior component of the GRF, this value was larger in magnitude throughout almost all of stance in P. waltl, whereas A. tigrinum exhibited two periods of slightly net positive values (Fig. 2, Table 1). Comparisons between the limbs of P. waltl demonstrate that the hindlimbs are the primary propulsors, whereas the forelimbs may be operating as ‘struts’ that act as pivot points about which the body can rotate, based on observations that the extended limb angles produce a straightened limb and that the net role of the forelimbs is deceleratory (interpreted from the negative values of the anteroposterior component of the GRF). Semi-aquatic and terrestrial salamanders both exhibit functional disparities between their forelimbs and hindlimbs but the extent of these differences is more distinct in semi-aquatic salamanders, potentially suggesting greater functional decoupling of the locomotor roles (e.g. acceleratory versus deceleratory) in forelimbs and hindlimbs in more aquatic taxa.

Biomechanical differences between forelimbs and hindlimbs can be affected by limb proportions, kinematics and muscle anatomy, but these are unlikely to be the primary factors driving the patterns in our results. Differential limb function in running quadrupedal lizards was caused by interspecific differences in limb morphology, with a more medial orientation of the GRFs produced by the hindlimb in lizard species that had disproportionately longer hindlimbs than forelimbs (McElroy et al., 2014). In our study, patterns of GRF production differed between the forelimbs and hindlimbs within each salamander species, but the ratio between forelimb length and hindlimb length was ∼1 in P. waltl (Trochet et al., 2014) and the ratio between the diaphysis of the humerus to the femur was ∼1.19 and ∼1.03 in A. tigrinum and P. waltl, respectively (Molnar, 2021). Moreover, the muscular anatomy and kinematics of terrestrial walking were suggested to be relatively conserved within salamanders, but did exhibit some interspecific differences that could be due to the terrestriality of a given species, or other factors related to the physiological challenges that salamanders experience in their respective environments (Pierce et al., 2020). Differences in ecological factors between our salamander species (e.g. habitual environments) could, thus, play a larger role in shaping GRF patterns than differences in anatomy or kinematics.

The functional differences between (i) pectoral and pelvic appendages and (ii) aquatic and terrestrial species could relate to the different demands imposed by the primary environments in which the appendages of these taxa function. For example, point estimates for the medial orientation of the peak GRF in semi-aquatic P. waltl limbs (−14 to −16 deg, forelimb–hindlimb) fall between those of P. barbarus pectoral fins (−17 deg) and most previously evaluated tetrapod limbs (≤−11 deg), including primarily terrestrial A. tigrinum salamanders (Fig. 8, Tables 1 and 2). A shift to a GRF directed less medially could reduce joint moments and, thus, the stresses experienced by the appendicular bones during terrestrial locomotion (Kawano and Blob, 2013). However, the greater medial orientation of the GRF in semi-aquatic P. waltl likely relates to the greater lateral spread of their distal limb segments compared with those of terrestrial salamanders, so that the feet are placed lateral to the elbow or knee joint during stance (Fig. 1B), rather than directly below these joints (as in terrestrial salamanders: Fig. 1C). A recent comparison of terrestrial locomotion in salamanders found that the maximum extensions of the wrist (160 deg), ankle (159 deg) and knee (155 deg) were higher in P. waltl than in A. tigrinum (150, 150 and 145 deg, respectively), and that the maximum abduction of the upper arm was larger in P. waltl (128 deg) than in A. tigrinum (110 deg) (Pierce et al., 2020). These kinematic differences could collectively contribute towards placing the distal elements of the limb more laterally in P. waltl. In addition, aquatic salamanders generally have thicker and longer articular cartilage on the long bones of their limbs (despite the forelimbs and hindlimbs being relatively subequal within taxa) and the thickness of the articular cartilage as a percentage of bone length was almost 2 times greater in P. waltl (52.3–43.7%, humerus–femur) than in A. tigrinum (24.1–24.8%, humerus–femur) (Molnar, 2021). Increasing the amount of cartilage on articular surfaces increases muscle length by moving the insertion point more distally, thereby increasing muscle leverage for limb extension (Molnar, 2021). Taricha torosa, another semi-aquatic salamander, had greater lateral spread of the limbs when walking in water than on land, which was caused by more extended angles about the limb joints (Ashley-Ross and Bechtel, 2004). Specifically, T. torosa had greater extension of the elbow and greater excursion of the elbow than A. tigrinum and P. waltl, while having comparable maximum extension of the knee and maximum abduction of the thigh (Pierce et al., 2020). Thus, certain aspects related to differences in the extensions in the proximal limb joints and abduction of the proximal limb elements (e.g. Pierce et al., 2020), articular cartilage (e.g. Molnar, 2021) and locomotor forces (this study) may collectively provide a stronger foundation for distinguishing semi-aquatic versus terrestrial salamanders. Given that this more pronounced lateral position of the distal limb segment is also found in the mudskipper fish, a more medial orientation of the GRF also might be found in other semi-aquatic taxa.

The observation that P. waltl generally has more extended joint angles than terrestrial salamanders on land may mean that some salamanders exhibit stereotypy (sensuWainwright et al., 2008) in their locomotor biomechanics. The wider spread of the limbs that would result from such lateral foot placement might convey additional stability against currents or other flows in aquatic habitats (Martinez et al., 1998) by reducing pitching and rolling (Chen et al., 2006). However, when on land, species that tend to be more aquatic may not be able to adjust to the more upright orientations of distal limb segments that are seen in terrestrial taxa (Kawano and Blob, 2013). Producing more acute angles about the more proximal joints could facilitate elevating the body off the ground, and shift the bone loading regime to reduce bending and increase compression (Ashley-Ross and Bechtel, 2004; Kawano and Blob, 2013), but the latter has yet to be confirmed empirically. Lateral spread of the distal appendage may also contribute to the high medial orientation of the GRF in mudskippers (Figs 1 and 2), and could prevent rolling of the body axis during terrestrial crutching, given the lack of extended posterior appendages in this taxon. Thus, shifts to less sprawling limb posture by flexing the distal limb segments closer to the body could have major biomechanical consequences that could facilitate terrestrial locomotion.

Similarities in the overall shapes of GRF and yank profiles produced by individual hindlimbs suggest that rates of GRF production are consistent in the timing of their changes across hindlimbs of semi-aquatic and terrestrial salamander species, but the different peaks observed in these curves indicate that the maxima differ between them. The similarities in force production between semi-aquatic and terrestrial hindlimbs could indicate that the use of the hindlimbs as a primary propulsor imposes strong constraints on kinetics for acceleration, but that other aspects, such as maneuverability and body support, can be altered by changing the GRFs and yank. For example, the smaller absolute value in the vertical component of the GRF and net GRF for the hindlimbs in P. waltl, compared with A. tigrinum, suggest that there may be steadier (i.e. less variable) loads in the former. Loading the appendages in more predictable ways can have numerous advantages (Granatosky et al., 2020). For instance, higher safety factors (margins of protection against failure) may be found in structures that experience more variable loading (Alexander, 1997; Lowell, 1985). Consequently, steadier loading could reduce the safety factors needed in the limb bones and, therefore, reduce the energetic costs to produce the higher safety factors. If load magnitudes were elevated as a result of the highly sprawled posture of P. waltl (Fig. 8), then it would be advantageous to have less variation in bone loading to reduce the safety factors needed. The limb bones of A. tigrinum have been found to possess generally high safety factors against failure (Blob et al., 2014; Kawano et al., 2016; Sheffield and Blob, 2011) and studies are ongoing to determine whether a similar pattern is found in semi-aquatic salamanders.

The rate of load application can also affect the load magnitudes experienced by bones. Strain magnitude is positively correlated with strain rate for both cyclic and discrete loading regimes in vertebrates (Aiello et al., 2015), so smaller yank values in P. waltl may be associated with slower strain rates and lower peak strain magnitudes in their limb bones. Lower peak strains would keep bones further from failure, whereas slower strain rates could affect bone mass (Rubin and Lanyon, 1982). Thus, peak strains could have a direct impact on the safety factors of the limb bones, and strain rates could affect the remodeling of bones in response to different loading patterns. Strain rate has been found to be a main factor driving bone remodeling and was positively correlated with higher Young's modulus values (Aiello et al., 2015), which would confer greater stiffness and greater resistance to bending in bones. Such properties would be advantageous for terrestrial taxa that must habitually withstand the downward effects of gravity on land. Semi-aquatic taxa that spend most of their time underwater might not need higher stiffness in their limb bones while their bodies are primarily being supported by buoyancy; however, such traits could be beneficial during the terrestrial eft phase of many newt species (including P. waltl). Preliminary analyses on the cross-sectional geometry of the limb bones in 93 salamander species indicate that the second moment of area at the mid-shaft is greater in terrestrial salamanders, suggesting a higher resistance to bending, compared with aquatic salamanders (Huie et al., in press). Changes in the mechanical properties and loading regimes of limb bones across ontogenetic stages within salamanders and between species with different life histories (e.g. direct versus indirect developers) are the subject of subsequent studies by our labs and would provide further insight into the factors that shape limb function.

Comparative studies on the locomotor biomechanics of salamanders living in different environments lend valuable information about the realm of locomotor capabilities that are possible within a general tetrapod Bauplan and point towards new avenues for examining the functional differences between and within limbs in the context of water–land transitions (Pierce et al., 2013; Standen et al., 2014). Broadening empirical analyses of extant taxa to include variation in ecology and habitat use is particularly crucial to identify which characteristics may confer greater terrestrial capabilities. Although limb kinematics and neuromuscular control during terrestrial locomotion appear to be relatively conserved across extant salamanders (Pierce et al., 2020), other aspects have been found to vary between aquatic and terrestrial salamanders. Terrestrial salamanders tend to have greater resistance to bending (Huie et al., 2022), shorter and thinner articular cartilage (Molnar, 2021), and decreased compactness via larger medullary cavities (Canoville and Laurin, 2009; Laurin et al., 2011) in their limb bones compared with aquatic salamanders. Amphibious salamanders are a particular challenge given that they sometimes cannot be differentiated from their terrestrial counterparts (Laurin et al., 2011) or have features that overlap with those of both aquatic and terrestrial salamanders (Molnar, 2021). Our work adds to this area of research by demonstrating how differences in the locomotor functions of forelimbs versus hindlimbs within a salamander can be less pronounced in a primarily terrestrial taxon, prompting future studies to examine whether the locomotor modules within salamanders exhibit different levels of integration in aquatic versus terrestrial species and environments.

Taxa with intermediate morphologies and locomotor biomechanics may be particularly relevant for improving understanding of stem tetrapods. Morphological comparisons across the humeri of 40 tetrapodomorphs spanning the water–land transition demonstrated that stem tetrapods had intermediate morphologies that placed them in a valley between two adaptive peaks (one for aquatic fishes and a second for terrestrial crown tetrapods) (Dickson et al., 2021). This was a likely transitional state before natural selection was strong enough to favor the movement of stem tetrapods to the adaptive peak representing the terrestrial characteristics now associated with crown tetrapods (Dickson et al., 2021). Moreover, Pederpes is an early tetrapod from the Early Carboniferous (348–344 million years ago) and had some terrestrial adaptations (e.g. forward-pointing foot) (Clack, 2009) but recent analyses indicated that there were also other traits that were indistinguishable from the Devonian early tetrapod (∼365 million years ago) Acanthostega, which was likely fully aquatic (Molnar et al., 2021), suggesting that early tetrapods in the early Carboniferous possessed a mosaic of tetrapod-like and fish-like traits that may have limited extended excursions on land. In addition, Greererpeton burkemorani is another early tetrapod from the Early Carboniferous (about 335–331 million years ago) that was presumed to be benthic and primarily aquatic (Godfrey, 1989). Our preliminary analyses on the locomotor capabilities of G. burkemorani indicated that estimated peak stresses in the femur decreased when modeling it with peak GRF magnitudes that more closely resembled those of the semi-aquatic P. waltl and as the angle between the GRF and the femur resembled more sprawling limb postures (Kawano and Blob, 2020). These studies highlight how incorporating a range of terrestrial capabilities can help broaden perspectives on the biology of stem tetrapods that likely still possessed a mosaic of tetrapod-like and fish-like features that are more indicative of a semi-aquatic taxon, rather than a fully aquatic or fully terrestrial taxon.

Our results provide novel insights into the functional changes associated with tetrapods becoming terrestrial, with broader implications for whether functional capabilities could be closely associated with major structural changes (e.g. transition from fin to limb) or could occur through gradual changes that included intermediate stages where taxa possess a mix of fish-like and tetrapod-like features. The kinetics of terrestrial walking in semi-aquatic P. waltl salamanders differed between the limbs, with the hindlimbs sharing numerous similarities with the hindlimbs of primarily terrestrial salamanders but the forelimbs exhibiting kinetics that were distinct from those of the pectoral appendages of both semi-aquatic fish and primarily terrestrial salamanders. Differences in limb kinetics between P. waltl and A. tigrinum could relate to the primary environments they inhabit as adults (i.e. aquatic and terrestrial, respectively), which complements published work demonstrating that limb morphology relates to a salamander's degree of terrestriality. Some kinetic patterns in P. waltl that were intermediate between those of semi-aquatic P. barbarus fins and terrestrial A. tigrinum limbs (e.g. mediolateral orientation of the GRF), which could be related to P. waltl possessing other traits (e.g. extended joint angles) that are intermediate between those of fishes and terrestrial salamanders. These data provide additional context to account for differential limb function in models of musculoskeletal function, including those that use salamanders as a blueprint for estimating the locomotor function of early tetrapods.

We thank Zach Marion, Jim Fordyce and Yu Xia for helpful discussions on statistical analyses that were considered for this study. We are also grateful to Margaret Ptacek, Miriam Ashley-Ross and the anonymous reviewers for feedback on earlier drafts of the manuscript. An earlier draft of the manuscript was submitted by S.M.K. in partial fulfilment of a doctoral dissertation at Clemson University.