Consequences of evolving limbless, burrowing forms for the behavior and population density of tropical lizards.

: We quantified functional traits (escape strategy, sprint speed and predatory performance) and population density across 10 lizard species representing the acquisition of burrowing snakelike morphs (BSLM), from Brazil. We used phylogenetic mixed models to test if: a) morph and substrate affects flight strategy and speed, b) BSLM species get more access to different potential prey types than lacertoid ones, when in syntopy, and d) morph is correlated with population abundance and habitat use in a way expected from the output of the previous experiments. BSLM rigidly relied on burrowing as flight strategy, while syntopic lacertoid species changed their strategy according to the substrate. Also, sand had opposite effects on sprint speed depending on morph, making lacertoids run slower and BSLM faster. Despite BSLM were overall slower than lacertoids, they were as good hunters of challengingly fast prey, and better hunters of underground prey. In their shared habitats, prey is most abundant in the leaf litter layer, although a large fraction is stock underground, too, under bushes. Experimental results support the observed higher importance of sand for BSLM’s density and the higher importance of vegetation for lacertoids’ density. Finally, although BSLM species reached the highest populational densities among the studied species, a systematic effect of morphological evolution on the abundance of limbless lizards remains elusive. prey for effects of vegetation and microhabitat on both traits using a mixed model with a Poisson link function, in which sampling point entered as grouping variable.


Introduction
For more than half a century, researchers have been striving to understand how organismal function relate to the abundance and habitat use of populations [1] [2]. Classical morphofunctional theory (AKA, the organismal performance paradigm [3] [4] states that evolutionary changes in individuals' morphology, behavior, or physiology interact with the environment to influence their performance (i.e. speed) during tasks such like escaping predators, capturing prey, etc. Within this paradigm, natural selection filters individuals depending on their performance during important tasks for survival (Ex. prey speed during predatory attacks, [3]. Literature shows many examples of selection favoring higher sprint speed [5] and how different morphologies induce different levels of sprint speed, or correlate with individuals´ microhabitat preferences (Ex. in lizards [6] [7][8] [9]. Nonetheless, slow morphotypes have also evolved in nature (Ex. the slow worm, or the sloth), raising the question of how this process relates to changes in morphology, behavior and population abundance. Likely, evolving a slower morphotype requires ways of compensating for the low speed during tasks such like escaping or hunting. That compensation might come from protective morphological traits (eg. A tortoises' armor) or compensatory escaping or hunting strategies, as shown by many vertebrates, like fishes [10] or birds [12]. For instance, fast lizard species have been observed avoiding "slower" microhabitats (Ex. thin branches), where fast locomotion is difficult [7]. Still, how behavior compensates the evolution of slower morphotypes remains largely unexplored.
In addition, understanding how the evolution of phenotypes relates to population level traits (eg. abundance), remains a key gap for linking patterns across different levels of biological organization [2]. Existing evidence remains heavily fragmented across different animal groups. For example, antipredator behavior has been shown to explain habitat use [12], but their abundance has been more related to access to food [13] [14]. In addition, a match between species' phenotypes and microhabitat use is often observed [15]. Thus, as phenotypes evolve, subsequent changes in individuals' performance during predatory interactions might have consequences for ecological patterns at the population level.
Even individuals' performance during experimental tasks may not be necessarily relevant for the traits of their populations in the wild. For example, while changes in morphology have experimentally shown to affect prey handling in limb reduced lizards [16], a posterior study did not find the expected restrictions of evolving limb reduced morphotypes on diet [17]. Therefore, to better understand relationships between morphology, performance, and population abundance, it is recommended to combine experimental observations with field observations (Ex. predatory performance with natural prey availability [18]). Yet this type of integrative study is more rarely done, due to the large effort and time required.
The evolution of burrowing snake-like morphs (BSLM) in squamates provides an interesting example to observe links between functional evolution and ecological traits at the population level. Across the world, BSLM have evolved multiple times from typically lacertoid species, developing elongated bodies, reduced limbs and a burrowing lifestyle [19] [20][21] [22]. This morphology may make BSLM relatively slow on the surface [23] but able to burrow faster and deeper than many lacertoids [24] [25]. The evolution of their predatory performance is not well understood, with some studies considering them as diet specialists [26]; [15] and others not [27] [16]. Many BSLM lineages have evolved in open sandy environments [28][29] [30] to which BSLM's distribution seems to be constrained [21] [27]. Still, some studies suggest they might be locally abundant [31]; [32]. Taken together, these studies suggest that morphofunctional evolution in BSLM could be linked to their population level traits, such as population density. However, these links have never been evaluated.
Herein, we report how the evolution of BSLM may influence their performance in different situations, and related their performance with their variation in population density. For that, we first tested whether morph and substrate can explain individuals' speed and flight strategy. Then, we tested if different lizard morphotypes had different predatory performance over four prey types (fast and slow crickets, and buried and surfaceactive termites). Finally, we tested if the availability of sand and vegetation affects the population density of each morphotype, and whether limbless or lacertoid forms can be more abundant at sites where they are syntopic.

Species accounts
Our study species comprehend two monophyletic sister groups of Gymnopthalmidae lizards. One is composed by BSLMs (Calyptommatus, Nothobachia, and Scriptosaura) and the other is compound of surface active lacertoid ecomorphs (Vanzosaura, Micrablepharus, Procellosaurinus, Psilophthalmus; [28], but capable of burrowing[. They represent the unique known acquisition of BSLM within the lizard tribe Gymnophthalmini (Fig. 1a, b). As happens in other BSLM lineages [29], our study lizards of both morphs may live syntopically. In this case, at a few sandy spots within the open biome of Brazilian Caatingas. All of the studied species are endemic from sandy areas in such region [28], except for Vanzosaura multiscutata that occurs on other soil types and Micrablepharus maximilliani that extend widely across Brazil. We acknowledge that a single evolutionary acquisition of such morphs cannot demonstrate causation, and thus we simply aim to show how traits correlate among each other in these species [3]. In addition, the natural "common garden" in which these endemic species have diverged makes our study system especially suitable for documenting the interactions between phenotypic evolution, environment and population level traits. During field trips, specimens were kept individually or in small groups, when found under the same bush, in transparent plastic bags with substrate from the exact place where they were collected. A lamp maintained thermal and light gradients within the bags. Animals were fed with termites and sprayed each two or three days. In the lab, animals were kept in 15x60x40 cm plastic terraria, separated by species but maintaining the groups found together in nature under single bushes during fieldwork. All terraria received UV (L12:D12) and heating (L10:D14) light, allowing lizards to thermoregulate normally. Termites were provided and water was sprayed in the terraria three times a week. Experimental and collection procedures were approved by the official ethics committee at the Instituto de Biociências da Universidade de São Paulo. Due to permit limitations, all experimental procedures were done over different, partially overlapped subsamples of 15 individuals, always including animals of both sexes and juveniles. Any specimen noticeably unhealthy (skinny or slow) was excerpted from the experiments. Specimens run for each species and test are shown at the corresponding results table.

Effects of morph on flight strategy
Flight strategy of all specimens was observed within 48h of their respective collection and at air temperatures between 30 to 35 C (within preferred ranges of all species, see ). ACG made lizards escape within an opaque plastic box (measures 13x35x29 cm), with its bottom covered by 3 cm of sand, and a half of the surface also covered by leaf litter from the collection sites. Each lizard performed only two consecutive escape trials in which it was deposited over one of two treatments: nude sand or sand plus the leaf litter layer, in random order. The immediate following action of the animal was coded as "burrowing", if the animal burrowed completely under the sand; "hiding", if the animal left a part of the trunk exposed or just sheltered under the leaf litter layer; or "running", if the animal kept running after the first 2 or 3 steps/body undulations. The flight strategy for each species within each context was coded as the average of the strategies observed among all individuals of that species.
We fitted phylogenetic logistic mixed models to test for the relative effects of morph and environment on the odds of using a given flight strategy. Since we had three escape strategies, we fitted three separate models comparing the three escape options in one pair each time. For that, we used the R function glmer, implemented in the package phyr [33]. We used the phylogenetic relationships proposed for these species [34]. Individuals entered as the grouping variable, to account for the two measurements made for each individual.

Effects of morph on sprint speed over different substrates
We estimated the maximum sprint speed of each species in two substrates (i.e. sand paper and loose soil). The sand paper substrate was constructed with dune sand glued to a wood card whereas the loose soil context was generated by adding a 2-3mm layer of the same sand to the sand paper. That depth provided enough lateral support for body ondulations while preventing lizards from burrowing under it. Lizards were maintained for a week in terraria prior to initiate speed measurements. Species' different thermal optima and loss of motivation were accounted for by making each lizard run in a random sequence of five different temperatures (22ºC, 27ºC, 32ºC, 37ºC, and 42ºC). Daily, each lizard ran twice per substrate at one temperature, during its typical activity period. Each trial was recorded using a high speed video camera (Redlake MotionXtra HG-SE, 250 frames per second) and illuminated with fluorescent lamps. This procedure generated 2054 videos over which we calculated individuals' linear speeds using the program Tracker 3.0 (Cabrillo edu®). We considered linear speed as the distance run between 3 steps/body undulations of acceptably rectilinear advance divided by the time to do it. For each temperature, the fastest run of each lizard at each substrate was selected to represent its speed.
We again used phylogenetic mixed models, this time using a gaussian link function, to predict speed in terms of morph and substrate, and controlling for phylogenetic relationships. We again used individual identity as the grouping factor, so as to compare the effects of treatment among identical individuals across the different temperatures. Modelling the effects of temperature would unnecessarily increase the complexity of the model and these are out of the scope of this study. Thus, we did not test for their effects.

Effects of morph on predatory performance over different prey types
We evaluated the predatory performance of all lizard species in four prey types. Slow versus fast prey, and surface-active prey versus buried prey. In all experiments, we stimulated lizard's appetite by only providing water for 4 days before trials. To change prey speed, we used captive raised cricket hatchlings. These crickets often present injured legs and thus allow comparing similar prey with different speeds. In this case, we added 5 crickets with no jumping legs and 5 healthy crickets at the same time to 1 l plastic boxes, each one containing 2cm layer of sand and a lizard. We kept each lizard and its prey for 24h in the plastic box and then counted the existing leg injured and healthy crickets.
To evaluate whether any morph was constrained to predate on superficial or buried prey, we placed individual lizards in plastic boxes containing 5 termites buried under two centimeters of clean sand from their habitats, which is reachable by all the studied species ). An overnight observation confirmed that termites were unable to neither resurge from this deep nor dig into the sand or leave the box. After 24h, we counted how many termites each lizard ate. On a different day, we repeated this procedure but with 5 termites left above the sand. This experiment was done in two steps to help noticing if lizards would surface a buried termite with their own movements.
For each species, predatory performance was measured as the numbers of prey items ingested, of each type. We tested for the effect of morph and prey type on the number of prey items ingested using phylogenetic mixed models, separately for each predation experiment (i. e. prey speed and buried prey). This time we used a Poisson's link function, adequate for count data [35]. In this model, the number of prey items ingested is predicted by morph and prey type, and the interaction between them.

Estimating the availability of feeding resources
We measured the abundance and diversity of potential prey in two natural microhabitats where these species live: superficial leaf litter and underground soil, and two vegetation cover types: shaded by vegetation and exposed soil. We collected 1 substrate sample of 1 L volume at each combination of microhabitat and habitat in 33 sampling points, across four of our lizard sampling sites, totaling 132 samples. Each sample of leaf litter was quickly collected in an area of less than half a squared meter over the ground, and then another liter of soil collected up to a depth of 10 cm (shovel length) within a m of the leaf litter sample. Each soil sample was examined in a white dish, within hours of collection. Because we observed failed predation attempts in the laboratory of pupae of 5 mm in minimum diameter, complete invertebrates smaller than that length in minimum diameter were counted. This criterion lead to count and identify 21 prey categories (mostly taxonomic orders and a larvae category). The trialed material was then preserved in alcohol. Using these datasets of prey diversity and abundance, we tested for effects of vegetation and microhabitat on both traits using a mixed model with a Poisson link function, in which sampling point entered as grouping variable.

Ecological distribution and density of syntopic lacertoid and BSLM species
We estimated each species' population density in ellipsoidal sampling units across all our five fieldwork localities. Following McCoy [27], ACG and an assistant raked the entire leaf litter surface of each sampling unit and the first 14 cm of sand existing under natural vegetation patches (i.e. bushes, small trees) and in similarly sized and shaped open patches between them. At each unit, we extracted and counted all the lizards found, and identified them to species level. Units were characterized as "shaded" or "exposed", depending on the presence of woody vegetation cover or just bare sand. We also labeled soil type as "no sand" if there was no sand between the leaf litter and the hard soil, or "sand", if there was a thick layer of sand of a least 5 cm under the leaf litter. To estimate each unit's area, we measured the two major perpendicular axes of each patch, and used those measures to solve the ellipsoid area equation. Using this procedure we sampled 277 units, raising data for all species of this study.
Due to the observed functional differences existing among morphs, we expected that the population density of BSLM species would depend more on sand availability, while lacertoids would be more generalistic with respect to sand, but still dependent on vegetation as refuge to the extreme thermal environment. We tested these expectations fitting a phylogenetic mixed model with gaussian link function, in which the density (abundance divided by plot size) of each lizard morph was the response variable, and habitat and sand availability were both binary fixed factors. Site entered as grouping variable.
Finally, to test for overall effects of morph evolution on the average density of each species, we fitted a phylogenetic mixed model (gaussian link function), using morph as fixed factor, and sampling site as random factor. For this test, we excerpted plots devoid of lizards.

Effects of morph and context on flight strategy
Across the three phylogenetic logistic mixed models, BSLMs escaped burrowing significantly more frequently than species of the lacertoid sister clade (PMM: trials=306, indi-viduals=198, species=10, range of z-scores=6.3-10.6, p values always <0.001, Table s1,a-c). Changing flight context, from being at bare sand to leaf litter had detectable effects for two of the three models. (PMM: trials=306, individuals=198, species=10, range of z-scores =(-4.56)-(-7.57), p values always <0.001, Table s1,a-c). Only the comparison of hiding versus running were not altered by the flight context (See table S1b). (PMM: trials=306, indi-viduals=198, species=10, z-value=0.4, p value=0.4, Table s1,a-c). In combination, these results suggest that BSLM maintain a rigid escape strategy while lacertoid species of gymnopthalmini lizards swiftly adjust their escape strategy to the environmental context ( Figure 2). Our non phylogenetic mixed model analysis detected similar effects and readers can re-run them using our supporting R scripts and data tables).

Effects of morph on sprint speed in different substrates
There was a significant interaction of morph and substrate (hard versus loose sand) (PMM gaussian: trials=972, species=10, slope=-17.54, SD=2.50, z-score=-6.85 p<0.001). BSLMs were significantly slower than lacertoids (PMM gaussian: trials=972, species=10, slope=-32.54, SD=2.27, z-score=-14.29 p<0.001), and opposite to them, sprinted faster on sand, Fig. 3)   Fig.3. Sprint speed of five lacertoid and five BSLM lizard species in different substrates. Each line connects the mean of maximum speeds of species measured across substrates. Violin width represents the frequency of species' values across the response variable.

Predatory performance and distribution of prey in nature
Our cricket hunting experiment showed clear effects of prey speed on hunting performance, indicating that healthy crickets were a difficult prey for species of both morphs (PMM poisson: trials=176, species=10, slope=1.19, SD=0.27, z-score=4.8 p<0.001). However, BSLM species consumed healthy crickets as well as lacertoid species (PMM poisson: trials=176, species=10, slope=0.53, SD=0.37, z-score=1.4 p=0. 15). For this experiment, there were not statistically detectable interaction between morph and prey type (PMM poisson: trials=176, species=10, slope=-0.22, SD=0.37, z-score=-0.59 p=0.55). See figure 4 (above) and table S3 for details. Fig .4. Above. Comparison of predatory performance between 5 BSLM and 5 lacertoid species eating healthy and leg injured crickets during a 24h hours experiment. Below. Results of an analogous experiment (134 observations) on the same species but comparing lizards' capability of hunting on surface prey versus prey located under sand (termites).
The termite hunting experiment also showed a significant interaction of morph and prey position (PMM poisson: trials=176, species=10, slope=-4.54, SD=0.33, z-score=-13.06 p<0.001), showing that buried termites were hard to access by lacertoids, while BSLM ate epigeal and buried prey equally well. See figure 4 (down) and table S4 for details.

Fig. 5.
Distribution of the number of prey types (above) and prey items (below) in different habitats (exposed, shaded) and microhabitats (soil surface, underground) in each of 33 sampling points, collected across 4 sites of northeastern Brazil.  Finally, lacertoid species did not systematically presented less dense populations than syntopic BSLM at the sampling units where they were found in syntopy (PMM gaussian: sampling points=98, Slope=-0.39, SD=0.83, z-score=-0.472, p=0.636).

Discussion
Many studies have focused on the variation of particular performance traits (The speed of ingesting prey, burrowing, or sprinting, [15][24] [23] or in the variation in microhabitat use among different species [8][28] [32]. Following Irschick and Losos [7], we integrated these different approaches to understand the evolution of limbless lizards more integratively. Our results suggest that at least among gymnophthalmini lizards, the evolution of BSLM implies in a loss of speed and flexibility in escape strategy, concomitant with their specialization in a burrowing strategy. These facts contrast with previous evidence, and the widespread view, that a higher sprint speed confers' individuals' fitness in lizards [5] [6]. This raises the question of how these a priori impairments in performance could have evolved. Risks of visually oriented predation and overheating are common at the open and sandy habitats where BSLM have most often evolved [24][25] [29]. As seen herein, sand and leaf litter typically slow down the sprint speed for lizard species of typical lacertoid species [36]. At these sandy habitats, the loss of grip and speed can lead to at least two opposite strategies, the evolution of morphological modifications that increase speed [37] and the acquisition of a slower, burrowing morphotype [28] [29]). These morphotypes have also been observed as more driven to bury, either when escaping predators, as shown by Australian lineages [24] or when avoiding thermal extremes, as found for Brazilian lineages [25] When it comes to foraging, gymnophthalmini BSLM species forage mostly during the night, using the sands' subsurface from where they surprise their prey [32]. Our results further show they can successfully attack challengingly fast prey in the surface, and that they can reach prey underground much better than syntopic lacertoid relatives. Thus, in an environment where sprint speed is difficulted, the evolution of slower morphotypes seems to be promoted when the environment offers opportunities to behaviorally compensate for a slower speed (ex. A soft substrate that allows sheltering and ambushing). In our study system, this might be the case, as sand not only allows BSLM's fixed flight strategy and access to a stock of underground prey, but also provides more options for thermoregulation [25].
The populational density of our studied species correlated with the availability of sand and vegetation cover in ways expected from existing theory and our behavioral experiments. First, the density of both morphotypes was affected by the presence of vegetation. A major factor for explaining lizards' density is access to food [12], and we found that, across sandy habitats of the Brazilian Caatinga prey is systematically concentrated in the leaf litter layers that exists under the vegetation. BSLM population density was arguably more dependent on sand availability than lacertoids, as we only observed for them a significant interaction of sand availability. Besides, vegetation impacted lacertoids' density more strongly (z-score of 13 for lacertoids vs 8 for BSLM), which is expected from their higher dependence of it to avoid thermal extremes [25].
BSLM's enhanced access to prey and to thermal refuges might lead to think that the evolution of these morphotypes may cause increases in population density. Our analyses provided mixed support for that hypothesis, though. While we did not find a systematic effect of morph on population density, the densest species were BSLM in all but one of our sampling sites. Specifically, the highlands of Gameleira, where sandy habitats are rarer. Yet, some of the BSLM species reached densities similar or lower than syntopic lacertoid sister taxa. These results suggest that an evolutionary process such as the acquisition of BSLM, despite involving a decrease in sprint performance and a more fixed escape strategy may help the new morphs to reach particularly high densities in sandy habitats, but that further factors may still restrict their population to levels similar to their sister taxa. This capacity to reach high population density within sandy habitats, supported by observations of other radiations [27], and the high functional dependence of BLSM on sand might explain why this morph has evolved so many times within Squamates [21], and its high proneness to vicariate following the fragmentation of sandy patches [38]. To know that for sure, integrative studies involving multiple radiations seem necessary.