Locomotion and palaeoclimate explain the re-evolution of quadrupedal body form in Brachymeles lizards

(a) Animal collection and morphometric data

We conducted locomotor trials with 147 individuals belonging to 13 species of Brachymeles plus Lygosoma bowringii in the Philippines and Thailand (electronic supplementary material, table S1). We captured animals using pitfall traps and by hand while raking leaf litter and detritus to uncover the animals [38]. Capture locations were marked using a handheld GPS unit (Garmin Ltd., Olathe, KS, USA). We used a pocket penetrometer (LR-281, Forestry Suppliers, Jackson, MS, USA) to measure load-bearing capacity of the soil [39] five times (taking the average) within a 50 cm radius of the exact site of capture for each animal. We also measured percentage soil moisture using a Kelway soil tester (Kel Instruments Co., Wyckoff, NJ, USA) for 161 specimens belonging to nine species of Brachymeles. We completed trials within 3 days of capture at base camp. Animals were then preserved and housed in the collections of the Sam Noble Museum of Natural History, Oklahoma, the Kansas University Biodiversity Institute, or the Zoological Museum of Kasetsart University, Thailand (electronic supplementary material, table S2).

We collected basic morphometric data from each preserved specimen. We counted the number of fingers and toes and took measurements using a Mitutoyo (Kanagawa, Japan) digital caliper to the nearest 0.01 mm. Measurements followed [40] and included: head length, width and height; snout–vent (SVL), tail, front and hind limb lengths; and body width. All measurements were taken in triplicate, using the average for analysis. We also photographed the head of each specimen from dorsal and lateral views using a Nikon D90 (Tokyo, Japan) camera with a Nikkor 60 mm macro lens. We then used ImageJ 1.52p [41] to measure head width and internasal distance from dorsal view, and head height and nasal height from a lateral view, similar to [42], to quantify snout pointiness and slope. Pointiness was calculated as internasal distance ÷ head width from dorsal view and nasal height ÷ head height from a lateral view, with smaller ratios indicating pointier snouts. These are the inverse of ‘rostral lateral decline' and ‘rostral angulation’ defined by [42]. The slope was calculated as (head width − internasal distance) ÷ head length from dorsal view, and (head height − nasal height) ÷ head length from a lateral view. Smaller slopes represented more gradually tapering snouts.

(b) Locomotor trials

We painted dots on each animal using non-toxic white paint dorsally at the occiput, pectoral girdle, mid-body, pelvic girdle, cloaca, elbows and knees [43]. The mid-body was the mid-point between the pectoral and pelvic girdles. Limbless species had subtle depressions where limb buds once developed, allowing us to identify girdle positions externally.

We conducted at least three running and three burying trials per individual on a fine and coarse substrate, resulting in over 2000 high-speed videos from the 147 individuals. Substrates were obtained by sieving (Hubbard #548, Forestry Suppliers Inc., Jackson, MS, USA) soil from the animals' habitat that we had first dried in the sun to ensure the same moisture content. The fine and coarse substrates were composed of particles 0.25–0.50 mm and 3–5 mm in diameter, respectively.

For running trials, we used a 15 cm wide by 1 m long field-portable racetrack with the substrate to a depth of 0.5 cm. For burying trials, we used a 31 × 18 cm plastic tub filled to a depth of 8 cm with the substrate. We recorded trials from dorsal view with a field of view of approximately 30 cm long, at a resolution of 512 × 384 pixels at 240 Hz and a shutter speed of 1/640 s using Casio Exilim cameras (Tokyo, Japan, models EX-ZR850 and EX-ZR1700). We coaxed animals to run or bury by gently tapping their hindquarters. A trial was successful when the animal either ran across the field of view without touching the walls of the track or buried continuously until the pelvic dot disappeared below the substrate. We conducted trials only when the ambient temperature was ≥27°C and measured body temperature of each animal immediately after each successful trial using an infrared thermometer (Raytek MT6, Santa Cruz, CA, USA) [28].

We also measured the penetration force exerted by each animal as it burrowed into each substrate at least three times. We used a 16 cm long plastic block with a 20 mm diameter tunnel drilled into it. One end of the tunnel was blocked with an aluminium plate mounted on a single element piezoelectric Kistler type 9203 force transducer connected to a Type 5995 charge amplifier (Kistler Instrument Corp., Amherst, NY, USA). The tunnel was filled with 5 cm of the substrate and the animal coaxed to burrow horizontally into this substrate, measuring the compressive force it exerted while burrowing at least 2 cm into the substrate [28,44].

(c) Quantifying running and burying performance and kinematics

We used DLTdv5 (Hedrick, 2008) in MatLab 2017a (MathWorks, Natick, MA, USA) to digitize the painted dots, providing us with xy coordinates for each point in each frame. The number of pixels per metre was calculated in ImageJ using a Philippine one Peso coin that was recorded in each video. We then calculated time from frame number and cumulative distance moved by the pelvic dot from the beginning of the video to each frame [43]. The pelvic dot approximated the centre of mass (CoM). We then used the curve-fitting toolbox in MatLab to fit a sixth-order smoothing spline to the cumulative distance and time data [39]. We adjusted the smoothing of the spline until secondary oscillations in the second derivative were removed and its maximum stabilized [43,45]. The first derivative provided us with frame-by-frame velocity estimates, from which we calculated maximum and average velocity while the animal was moving. We selected the trial with the highest average velocity for each individual per a mode of locomotion per substrate to avoid pseudoreplication and calculated kinematic variables for these trials.

We focused on axial kinematic variables because some species were limbless. For running, we calculated the average frequency, wavelength and amplitude of lateral undulations for the pectoral, mid-body and pelvic dots. We did this in R, using the pectoral point to model a linear path that the animal travelled and calculating the undulatory variables based on that path, as in [46]. Animals do not undulate regularly during burying and disappear from view as it progresses, so for this behaviour, we tallied the proportion of trials that each individual used its front and hind limbs just before each set of limbs disappeared beneath the substrate [47].

(d) Statistical analysis

We conducted all analyses in R v. 3.5.1 [48]. All analyses were phylogenetically informed and we used the most recent dated phylogeny of Brachymeles, which was based on a multi-locus gene dataset, and with most clades being supported with posterior probabilities of one [33]. A phylogenetic principal component analysis (pPCA) was used to study patterns of morphometric variation using the ‘phytools' package [49,50]. The pPCA was done using a correlation matrix, simultaneously estimating λ [51], a measure of phylogenetic signal in the data, including the number of digits (equal in front and hind limbs in our species), ln-transformed head length, and relative SVL, body width, and front and hind limb lengths. We calculated relative lengths by dividing each by head length, following other studies of body elongation and limb reduction [14,29,52]. We chose not to standardize by body mass because we found that in Brachymeles, body mass correlated highly with body form, as represented by pPC-1 (R² = 0.738, slope = −0.049 ± 0.008, t = −5.81, p < 0.001) and made the pPCA difficult to interpret. Head length correlated strongly with body mass (R² = 0.959, slope = 3.048 ± 0.182, t = 16.71, p < 0.001), indicating that it was an appropriate measure of body size. To understand how intraspecific variation mapped on our species PC space, we projected individuals onto the phylomorphospace [53]. There was relatively little intraspecific variation along the pPC-1 axis (electronic supplementary material, figure S1). Our pPCA results closely matched those obtained using a covariance matrix, and one where we decomposed the hind limbs into their constituent segments, supporting the robustness of the analysis (electronic supplementary material, table S3). Hence, we used pPC-1 as an index of body form for subsequent analyses.

We also modelled the evolution of size, represented by head length, for 40 species of Brachymeles and our phylogeny [33,40]. We did this by fitting a series of continuous trait evolution models in the ‘OUwie' package [54] and comparing them using AIC corrected for small sample size, interpretting a Δ_AIC > 2 as support for one model over another [55]. To test whether the pentadactyl clade evolved towards a different size optimum than the limbless and limb-reduced species, we fitted a single rate Brownian motion model, a single optimum OU model and two OU models with multiple optima. The first had an optimum for the pentadactyl clade, including its most recent common ancestor and an optimum for all limbless and limb-reduced species. The second had three optima: one for the pentadactyl clade, one for limb-reduced species and one for limbless species. We then compared whether the optimum size for the pentadactyl clade differed from the optimum for the other species using a t-test.

We tested for relationships between variables using phylogenetic generalized least-squares regression (PGLS) while estimating λ of the residuals [56], using the ‘caper' package [57]. For analyses of locomotor variables, we conducted separate regressions for coarse and fine substrates. We tested for relationships among pPC-1 factor scores, head shape indices, locomotor and soil variables using PGLS. Locomotor velocity, and the amplitude and wavelength of body undulations could be affected by body size, so we also repeated analyses with these variables divided by HL to adjust for size. Although our regressions were planned a priori, there were many comparisons, so we corrected for multiple comparisons while accounting for false discovery rate [58,59]. Because this is not implemented in ‘caper', we did these corrections post hoc, so some p-values that we present are <0.05 yet are interpreted as not significant.

(a) Evolution of body form in Brachymeles

Our pPCA characterized body form variation in Brachymeles, with the first component explaining 77% of this variation. pPC-1 provided an index of body form ranging from lizard-like to snake-like. Species with high pPC-1 scores were small and had relatively long, narrow bodies and relatively short limbs with few digits, so were more snake-like than species with lower scores (electronic supplementary material, table S3). Species with low pPC-1 scores were large, had short, thick bodies and were pentadactyl. Projection of individual data into the species PC space revealed that intraspecific shape variation was small along the pPC-1 axis, further supporting its use as an index of body form (electronic supplementary material, figure S1).

These body form variables are related similarly in other squamate clades [14,24,40,60], and the elongation of the body and reduction of appendages may be the primary axis of body form variation in vertebrates [61,62]. However, body size and body form relate differently across exemplar clades. In Lerista, both limbless snake-like and quadrupedal pentadactyl species are diminutive while intermediate species are large [29]. By contrast, Brachymeles head length, the proxy for size in elongate squamates, was correlated tightly with pPC-1 (PGLS regression: R² = 0.631, slope = −0.137 ± 0.030, t = −4.53, p < 0.001), indicating that more lizard-like species were larger than more snake-like ones.

In investigating the role of body size in the re-evolution of a quadrupedal, pentadactyl body form, the Ornstein–Uhlenbeck (OU) model with two selective optima best fit the data, followed by the OU model with three optima (electronic supplementary material, table S4), both indicating that the pentadactyl clade, including its ancestor, evolved towards a significantly higher size optimum (HL = 25.6 ± 5.7 mm) than species with reduced or absent limbs (HL = 5.1 ± 0.9 mm, t = 3.58, p < 0.001). Thus, body size appears to be an important part of the story behind the re-evolution of a quadrupedal, pentadactyl form in Brachymeles.

(b) Locomotion and the re-evolution of quadrupedal form

To study the coevolution of locomotion and body form in Brachymeles, we first tested whether there was a correlation between burying and surface locomotion performance. We found support for the facilitation hypothesis because there was a significant positive relationship between average burying and surface velocity on the coarse substrate, with a similar, marginally non-significant trend on the fine substrate (figure 2a; electronic supplementary material, table S5A). This was driven largely by increased maximum and average surface velocity in more lizard-like species (figure 2c,d; electronic supplementary material, table S5B), because there was only a non-significant trend between average burying velocity and body form (figure 2b; electronic supplementary material, table S5C). These results were unchanged when we analysed size-corrected locomotion data (electronic supplementary material, table S5D). Opposite to what we show here for Brachymeles, in Lerista, body form is strongly related to burying performance but not surface locomotion [28].

• Download figure
• Open in new tab
• Download PowerPoint

Animals increase surface locomotor speed by increases in the stride or undulatory frequency of the locomotor cycle, or the distance moved during that cycle [63,64]. To understand how more lizard-like species moved faster on the surface, we related body form to the characteristics of lateral undulation of the animals. We found that more lizard-like species adopted higher-frequency undulations, but with no change in undulatory amplitude or wavelength, resulting in their higher velocity compared with more snake-like species (figure 3; electronic supplementary material, table S5B). These patterns were also unchanged when considering size-adjusted amplitude and wavelength (electronic supplementary material, table S5D).

• Download figure
• Open in new tab
• Download PowerPoint

Locomotor data for Lygosoma bowringii also supported the facilitation hypothesis. This species had the lowest average burying velocity of any pentadactyl species on the coarse substrate, a mid-range burying velocity on the fine substrate, yet very high average surface locomotor velocity and undulatory frequency on both substrates (figure 2a, 3a). This makes sense because L. bowringii evolved from a cursorial, non-burrowing ancestor [65,66], while ancestors of pentadactyl Brachymeles have been estimated as snake-like burrowers using phylogenetic evidence [14,33].

Our facilitation hypothesis posits that re-evolved limbs might aid in fossorial locomotion in tougher soils. Wet substrates, like those inhabited by Brachymeles, have at least a 400% higher resistive force than comparable dry substrates [67]. Supporting our hypothesis, we found that species with limbs used them in the vast majority of trials during burying. This was especially true for the front limbs (figure 4a; electronic supplementary material, table S5C). Although L. bowringii has ancestral limb morphology, it lives in wet substrates similar to Brachymeles [37], and usually used its front limbs (>90% of trials, figure 4a), but rarely its hind limbs (<20% of trials, electronic supplementary material, figure S2A) while burying. The phalanges of limbed Brachymeles re-evolved to be more robust, possibly allowing them to exert more force on the substrate, something not seen in Lygosoma [33]. Although some lizard species fold their limbs against their bodies while burrowing [63], Chalcides occelatus, a skink that often inhabits damp sand, also uses its front limbs while burying [67]. All of these lines of evidence suggest that wet substrates may impose different selective pressures on fossorial animals than dry substrates.

• Download figure
• Open in new tab
• Download PowerPoint

We found that more snake-like species of Brachymeles were able to initiate burying while exerting less force, particularly on the fine substrate (figure 4b). Narrower objects require less force to penetrate a substrate [68,69], and snake-like Brachymeles have narrower heads (electronic supplementary material, figure S2B). Other aspects of head shape may also affect burying performance. For example, the sandfish skink (Scincus scincus) and South American gymnophthalmid species that burrow in looser soil have pointy, shovel-like snouts [42,70,71]. We tested this idea and found that more snake-like species had more gently sloping snouts (electronic supplementary material, figure S2C), but that pointiness was unrelated to body form (electronic supplementary material. figure S2D, table S5E). Despite the gradation in the slope of the snout, the head shape was unrelated to burying velocity or penetration force (electronic supplementary material, table S5F).

One reason for the lack of relationship between head shape and burying performance may be that all species inhabit equally firm substrates. Indeed, we found no relationship between either of the load-bearing capacity or the moisture of the substrate and body form (pPC-1), snout pointiness or slope, supporting this hypothesis. This is unlike the situation in gymnophthalmids [42,72], but Brachymeles lives in wetter soils with more plant debris.

(c) Palaeoclimate change and the re-evolution of quadrupedal form

Our locomotion data suggest a scenario in which Brachymeles evolved to be snake-like when the climate and resulting substrate properties were more amenable to fossoriality, and subsequently re-evolved a quadrupedal, pentadactyl form when the substrate became more wet and packed. Many other clades of snake-like lizards live in arid areas with dry substrates. For example, snake-like species of gymnophthalmids live in more arid regions with drier soils than quadrupedal species [72]. The same appears to be the case with South African Scelotes scincid lizards [73].

To investigate the relationship between body form evolution and climate change, we compared published palaeoclimatological and geographical data with our latest dated phylogeny of Brachymeles (figure 1) [33]. Our phylogeny estimates that Brachymeles evolved 62 ± 15 Ma, which was before the Philippines formed (figure 1), indicating that its origin was on mainland Southeast Asia when it extended to present-day Borneo [74]. Palaeoclimatological data are sparse for this time period, but suggest a drier and more seasonal climate than today [75]. From 50 MA to present, oxygen isotope, temperature, pollen and fossil algae data are available, resulting in more reliable reconstructions [75]. These data show that between 50 and 30 Ma the climate tended to be subhumid and seasonal with coniferous trees expanding during cooler periods (figure 1) [75,76]. Hence, Brachymeles probably evolved in habitats that were much drier, with looser substrates that were easier to penetrate than today.

Brachymeles’s origin in the Philippines was 31 ± 6 Ma, shortly after many Philippine islands separated from the Sunda Shelf [74] (figure 1). At this time, the climate remained hot and dry [75,77]. However, by 25 Ma, the climate had shifted to be much wetter, with the monsoons becoming established [75–77]. This change in climate would have resulted in wet, packed soils with more vegetation that are more difficult to penetrate [67], similar to today's environments. The pentadactyl clade evolved around 21 ± 5 Ma (figure 1), coinciding with this transition to a monsoonal ‘ever-wet' climate, when robust limbs and larger body size would facilitate burrowing in such soil, as our locomotion data show. These large-scale palaeoclimatological data cannot take into account microhabitat selection of drier or wetter soils as the climate changed. However, our field soil moisture and load-bearing capacity data showed no relationship with body form, despite a possible expectation of more snake-like species inhabiting less compacted or drier soils.

(d) On the persistence of snake-like forms

We showed that limbless species of Brachymeles are poorer surface and fossorial locomotors than quadrupedal species and live in soil that is not well suited to their locomotion. These species are often syntopic with large pentadactyl species [34], and so an outstanding question is how they persist. Our observations in the field suggest that the diminutive, snake-like species are more secretive and slower moving than large pentadactyl species. The snake-like species exclusively feed on arthropods that live in the leaf litter, while the pentadactyl species also prey on vertebrates, including lizards and snakes [78,79]. We also primarily captured pentadactyl species in pitfall traps, suggesting that they are more surface dwelling than snake-like species. Hence, high levels of performance may be less important to limbless Brachymeles than pentadactyl species on account of the limbless species being less frequently be exposed to predators due to their fossorial and secretive nature [80,81].

Size and dietary differences among species differing in body form may have also allowed for the persistence of intermediate forms in Australian Lerista skinks, but in that clade, the intermediate forms that are largest [29]. Quadrupedal species of Lerista tend to have diets of more surface-dwelling insects than more snake-like species [82] and tend to eat smaller prey [83], so in both clades, there appears to be a greater reliance on surface locomotion in quadrupedal species, further supporting a link between snake-like forms and fossoriality.

Divergence in size and activity patterns may have led to a divergence in a dietary niche, where high locomotor performance is necessary to capture small vertebrates. Many vertebrate radiations initially diverge in microhabitat followed by trophic morphology [84], but snake-like and quadrupedal, pentadactyl species of Brachymeles can be found in the same rotten log, yet may have diverged in size, diet and associated locomotor behaviour. Further work is needed to test this hypothesis, but similar patterns of diversification and coexistence have been documented in cichlid fishes [85].