Early burst of parallel evolution describes the diversification of gecko toe pads

Introduction

Together with describing the breadth of earth’s biodiversity, a core interest of biologists is understanding where the diversity came from and what processes generated it, with a particular interest in its repeatability (Gould, 1990). The likelihood of repeatability, or predictability, of evolution depends on the interplay between deterministic and stochastic processes. There is a rich body of literature describing how various evolutionary processes may or may not lead to the repeated evolution of traits (Lenormand et al., 2009; Losos, 2011; Bolnick et al., 2018; Buscher and Gorb, 2021). While the evolution of similar traits does not serve as sufficient evidence of deterministic selection, the repeated evolution of particular phenotypes invites further investigation to better understand how selective, stochastic, functional, or developmental factors may have generated an observed pattern.

Repeatedly evolved traits are often considered to be convergent, yet the history of “convergent” and similar terms is complicated. Losos (2011) defined convergent evolution as the “independent evolution of similar features in different evolutionary lineages” with later studies including the terms “parallel evolution” and “repeated evolution” (Bolnick et al., 2018; Stuart, 2019; Cerca, 2023). Cerca (2023) provided a general term, “repeated evolution”, to include the evolution of similar traits regardless of ancestral phenotypes. Bolnick et al. (2018) and Stuart (2019) advocated for a trait-based geometric set of definitions, one in which trajectories through trait space are compared. A taxa’s phenotypic variation represents a trajectory though trait space. By comparing the trajectories of different taxa, we can ask if taxa reside in parallel areas of trait space. Further still, if trajectories have a directionality, for example, change from an ancestral phenotype to a derived phenotype, we can compare evolutionary trajectories and hence quantify how parallel the changes were in different taxa. Under this definition, Bolnick et al. (2018) and Stuart (2019) suggest parallel evolution occurs when lineages evolve (change from ancestral to derived phenotypes) along similar directions, creating parallel trajectories of change. Cerca (2023) described such patterns as colinear trajectories. Bolnick et al. (2018) and Stuart (2019) defined convergent evolution, or confluent trajectories under Cerca (2023), as occurring when derived phenotypes are more similar in trait space than ancestral phenotypes, and divergent evolution, or anti-collinear trajectories following Cerca (2023), as occurring when derived phenotypes are more different than ancestral phenotypes. The above definitions of convergent/confluent, divergent/anti-collinear, and parallel/collinear evolution are the definitions we are adopting in this study.

Studies of repeatedly evolved traits often focus on ecologically relevant morphology, especially in the context of adaptive radiations, (e.g. limb length in Caribbean Anolis lizards (Losos, 2009), body and mouth shape in African cichlid fish (Muschick et al., 2012) and three-spine stickleback fish (Rundle et al., 2000)). In addition to these classic systems, we suggest the adhesive toe pads of lizards, especially geckos, are an underappreciated example of repeated evolution. Adhesive pads have evolved multiple times, including multiple clades of arthropods(Federle, 2006; Buscher and Gorb, 2021) and lizards. The pads of lizards, including geckos, anoles, and some skinks, are considered dry fibrillar adhesives; a novel example of naturally occurring biological nanotechnology (Autumn and Gravish, 2008). Using very small (<150 μm long) hair-like structures called setae (Ruibal and Ernst, 1965; Garner and Russell, 2021), pad bearing lizards adhere to surfaces using van der Waals forces and/or contact electrification (Autumn et al., 2002; Izadi et al., 2014). Origins of adhesive pads within squamates occurred in stem anoles (Russell, 2017), Prasinohaema and Lipinia skinks (Williams and Peterson, 1982), and multiple times in geckos (Gamble et al., 2012; Russell and Gamble, 2019), including strikingly similar adhesive tail pads (Bauer, 1998; Griffing et al., 2021), with multiple taxa secondarily losing toe pads as well over the clade’s 200-million-year history (Russell and Gamble, 2019). The similarities within and between gecko and anole toe pads have been discussed at the microscopic setal level (Ruibal and Ernst, 1965; Williams and Peterson, 1982; Russell, 2017; Russell and Garner, 2021) and at the macroscopic toe pad level (Russell, 1979; Han et al., 2004; Russell and Gamble, 2019). The breadth of gecko toe pad morphology includes many examples of presumably repeated morphologies, as evident in the naming of genera, for example Thecadactylus and Pseudothecadactylus, two distantly related genera with nearly identical external toe pad morphology (Brongersma, 1934; Russell, 1979; Han et al., 2004; Gamble et al., 2012; Russell and Gamble, 2019). Despite their external similarities, these genera have unique internal anatomies (i.e. paraphalanges) (Russell and Bauer, 1988). While toe pad morphology varies extensively across genera, morphologically similar yet distantly related genera are very common in Gekkota. Sister genera often have very different toe morphologies, while distantly related genera often share similar morphologies (Russell and Gamble, 2019). Conversely, within genera, toe pad morphology appears relatively consistent, with little variation across species.

Historically, gecko toe pad morphologies have been categorized into subjective toe pad classes, such as basal and distal toe pad classes (Riedel et al., 2021) or leaf, fan, and basal pad classes (Russell and Gamble, 2019), suggesting the repeated evolution of similar morphologies. To evaluate these patterns using quantitively defined morphologies and gain a better understanding of how gecko toe pad diversity evolved, we first determined if toe pad morphologies among geckos are morphologically convergent, parallel, or neither using quantitative approaches. We then developed an approach to evaluate models of trait evolution using high-dimensional morphological data to investigate the evolutionary patterns gecko toe pads likely evolved under.

To evaluate the similarities in toe pad morphologies across geckos, without solely relying upon assigned toe pad classes, we aimed to quantify toe pad shape using geometric morphometrics and phenotypic change vector analyses (Collyer and Adams, 2013; Bolnick et al., 2018; Stuart, 2019; Cerca, 2023). We hypothesize our geometric morphometric shape data will not only capture assigned toe pad classes, but that gecko toe pads will exhibit multiple examples of parallel/collinear morphology, assuming a similarly padless ancestor (Russell, 1979), suggesting each toe pad class evolved independently in each family, yet repeatedly across families (Figure 1). Specifically, we predict that trajectories within each family, from padless morphologies to distal morphologies will be parallel/collinear, with similar angles and magnitudes through morphospace. We also predict transitions from padless to our included incipient morphologies and transitions from padless to basal morphologies will all be parallel/collinear across families, yet the magnitudes from padless to incipient phenotypes will be shorter, while still collinear to transitions to basal pads (Figure 1) due to their apparent similarities in shape.

We also considered what models of trait evolution would best capture our observed pattern of toe pad diversification using a modified approach calculating standardized phylogenetic independent contrasts and a node height test (Felsenstein, 1985; Freckleton and Harvey, 2006). If the diversification of gecko toe pad shape fit a pattern of Brownian Motion, i.e. a constraint rate of morphological diversification through time, this would suggest unconstrained and random evolution possibly via stochastic processes like drift or fluctuating selection (Harmon, 2018). Alternatively, the diversification of gecko toe pads may better fit an adaptive or deterministic model of trait evolution, such as early burst or late burst (Harmon et al., 2010) with the rate of morphological evolution slowing through time (early burst) or accelerating through time (late burst). Considering sister genera often fall into different toe pad classes, while distantly related genera often share similar morphologies, these patterns suggest high rates of morphological change early in the history of Gekkota. In addition, toe pad morphology appears relatively consistent within genera, suggesting slower rates within genera. As a result, we hypothesize that toe pad diversification will not fit Brownian motion but will instead be consistent with an early burst pattern of evolution, with rapid diversification early in the history of the clade (Simpson, 1944; Harmon et al., 2010), concurrent with the diversification of extant genera, fixing toe pad shape within genera. Together we anticipate a pattern of early fast parallel diversification from padless to repeated pad bearing morphologies. Overall, this pattern of toe pad evolution would suggest adaptation driving the initial diversification and repeated evolution of gecko toe pads.

Materials and methods

Data collection

We included species to capture toe pad diversity across all six limbed families of Gekkota: Diplodactlyidae (17 species), Carphodactylidae (2 species), Eublepharidae (3 species), Sphaerodactylidae (9 species), Phyllodactylidae (10 species), and Gekkonidae (52 species). We also included specimens of Sphenodon punctatus, a reptile species sister to all squamates, to provide further information as to the morphology of a likely padless ancestor. Our dataset included 250 specimens from 94 species (56 genera), with an average of 2.7 specimens per species (min = 1, max = 13, standard deviation = 2.8). The wide range of species allowed us to capture toe pad variation across and within genera. We strived to include representatives that captured as many toe pad origins or class transitions within each family as possible, also including multiple species within genera that contained transitions or origins. For example, we included species with and without toepads in Lucasium to capture the loss of pads in that genus and multiple species within Gonatodes and Cyrtodactylus to capture independent origins of pads in these genera (Russell and Gamble, 2019).

We photographed the ventral side of each specimen’s left or right hind foot using digital cameras or dissecting microscopes. Our dataset included live and preserved specimens. Preserved specimens were imaged while submerged in ethanol, pressed under a glass microscope slide to flatten their feet. Live animals were imaged with their feet pressed against a vertical pane of glass. Photographs of right hind feet were horizontally flipped prior to analyses so that our entire dataset appeared as left hind feet.

Although we quantitatively measured toe pad shape using geometric morphometrics, we also organized species into morphological classes to test how parallel, convergent, or divergent presumably similar morphologies were. We defined four toe pad classes: padless, distal, incipient, and basal (Figure 2), based on each species’ ventral scale morphology on their fourth rear toe. We developed these designations similar to previous work (Russell et al., 2015; Russell and Gamble, 2019; Riedel et al., 2021; Griffing et al., 2022). We defined our padless class as toes with no obvious adhesive setae and no exaggerated or dilated ventral scales. Padless species had granular or transverse ventral scales. Species with the appearance of toe pads were categorized into three classes: distal, incipient, or basal. Species with an isolated distal pair of pads (or single asymmetrical pad in the case of Sphaerodactylus) were placed in the distal category, even if they also had additional smaller adhesive scales more proximally on the toe. This category includes species often called “leaf toed geckos” (e.g. Phyllodactylus). We defined our incipient toe pad class as possessing moderately or rudimentarily exaggerated, widened, or dilated ventral scales situated near the base or middle of the toe, with mild dilation of the silhouette of the whole toe, similar to previously described frictional plates in some Gonatodes (Russell et al., 2015; Higham et al., 2016). Lastly, basal pads were considered all other well-developed pad morphologies. As a result, our basal class contains a wide range of morphologies, including species with many narrow lamellae covering much of toe, traditionally considered to have “basal pads” similar to Riedel et al. (2021). Our basal category also included species typically considered to be “fan toed” species, such as Ptyodactylus and Uroplatus (Russell and Gamble, 2019), despite the fact that their many narrow lamellae are concentrated at the distal end of their toe, and species with multiple pairs of adhesive scales across their toe (like Hemidactylus) that lack an exaggerated, clear, or isolated distal pair of pads.

Our four toe pad classes are based on different suspected evolutionary pathways to toe pads (Griffing et al., 2022). One pathway, likely underway in Heteronotia geckos (Riedel et al., 2021) is the evolution of adhesive scales at the very tip of the toe, on a pair of scales flanking the claw, leading to the evolution of morphologies in our distal category. Alternatively, adhesive pads may originate over a larger area of the toe (Griffing et al., 2022), possibly acting initially as “frictional plates” as observed in Gonatodes and Cyrtodactylus (Russell et al., 2015). This morphology is represented by our incipient toe pad class, which may lead to the evolution of more well-developed morphologies covering multiple scales, represented by our basal class.

After assigning each species to a toe pad class, we used geometric morphometrics to capture toe silhouette and ventral scale shape of each individual specimen. We examined the fourth hind toe (often the longest toe) of each imaged foot, following the work of previous studies (Schulte et al., 2004; Howell et al., 2022; Michaud et al., 2023), placing landmarks using the TPSdig software (Rohlf, 2010). Our landmark scheme included 11 fixed landmarks and 10 curves connected to fixed landmarks at each end, with 13 semi landmarks along each curve (Figure 3). The landmarks were chosen to highlight functionally relevant features of the toe that are repeatable across padless and pad bearing species. Landmarks 1 and 2 were placed on the left and right side of the base of the toe where the toe meets the palm. Landmarks 3 and 4 were placed on the left and right side of the toe adjacent to the widest scale or pair of scales. Landmark 5 was placed at the distal end of the toe at the ventral base of the claw. Landmarks 6-11 focused on the widest three consecutive ventral scales. In pad bearing species, dilated adhesive scales are termed lamellae or scansors. Landmarks 6 and 7 were placed on the left and right side of the most distal of the three widest scales. Landmarks 8 and 9 were placed on the left and right side of the middle of the three widest scales. Landmarks 10 and 11 were placed on the left and right side of the most proximal of the three widest scales. Curves between fixed landmarks were drawn to outline the silhouette of the toe (landmarks 1 to 3, 3 to 5, 5 to 4, 4 to 2) and to outline the three widest scales (landmarks 6 to 7, 7 to 6, 8 to 9, 9 to 8, 10 to 11, and 11 to 10). Each curve was represented initially by 15 evenly spaced semilandmarks, but the first and last semilandmark from each curve was removed due to their overlap with existing fixed landmarks, resulting in each curve having 13 semilandmark. In total, each image initially had a total of 161 landmarks, but this was reduced to 141 landmarks per specimen. Our decision to capture the shape of the three largest consecutive lamellae/scales represents a compromise. It is not feasible to quantify all of the ventral scales/lamellae on each specimen, which vary extensively in quantity across specimens. While our three focal scales are likely not developmentally homologous across all our focal species, they are repeatable. In addition, there is some evidence that the distal scales of geckos with pads may in fact be developmentally homologous (Griffing et al., 2022).

For all phylogenetic related analyses (generating PCA phylomorphospace plots, phylogenetic Procrustes ANOVA analyses, and calculating phylogenetic independent contrasts) we used the ultrametric fully resolved DNA-based squamate phylogeny from Tonini et al. (2016), trimmed to only include our focal species (Supplementary Figure S1). We had morphological data for six species that were not found in the Tonini et al. (2016) phylogeny (Aristelliger expectatus, Sphaerodactylus inaguae, Phyllodactylus gerrhopygus, Goggia microlepidota, Goniurosaurus orientalis, and Rhynchoedura eyrensis). In these situations, closely related tips in the tree were renamed to match our focal species. When our sampling contained a single species in its genus, which was often the case, the name swap had no phylogenetic impact on our analyses. When considering genera in which we sampled multiple species, and no other information as to which available tips were most closely related to our focal taxa, we chose tips to rename that maximized the phylogenetic distance between focal species.

We imported our landmarking into Rstudio (v. 2022.12.0.353, R v. 4.1.2). We used the R packages geomorph (v. 4.0.1, Adams and Otarola-Castillo, 2013) and Morpho (v. 2.9, Schlager, 2017) to manipulate and analyze geometric morphometrics data. We conducted a Procrustes alignment, optimizing semilandmark locations using minimizing bending energy, and conducted a variety of diagnostic analyses to assure our specimens were accurately landmarked (Supplementary Figure S2).

Analyses

Our initial analyses included a species-level principal components analysis (PCA) to visualize the general variation across our shape dataset, using the “groupPCA” function in the Morpho package (Figure 4). A between-group PCA was chosen to emphasize differences between species (i.e., the designated “groups”) while minimizing the impact of within-species variation and our uneven sampling. We also included a projection of our species-level phylogeny in our PCA plots, estimated the proportion of total variation captured by the first four PC axes, and projected the shape each axis captured, represented by the shape one standard deviation above or below the mean shape of each axis.

To investigate how accurately our a priori toe pad class assignments captured variation in our shape dataset, we conducted a series of analyses including a phylogenetic Procrustes ANOVA, correlating species mean toe pad shape with our assigned toe pad classes, while accounting for the non-independence of species. We also used a cross-validation approach, as part of the “groupPCA” function in Morpho, conducting a between-group PCA using our toe pad classes as groups, to evaluate how accurately specimens could be assigned to their proper toe pad class and statistically comparing the shapes of specimens assigned to each class. Lastly, we used a morphological clustering approach. We estimated Procrustes distances between specimens (a Procrustes distance matrix) via the “gpagen” geomorph function and performed a hierarchical cluster analysis, using Ward’s 1963 clustering criterion to build a dendrogram based on morphological similarities (Supplementary Figure S4; Murtagh and Legendre, 2014).

Parallel morphologies

Given the apparent similarities in toe pad morphologies across families, our goal was to quantify these similarities. There is strong evidence geckos repeatedly evolved and lost toe pads and as a result, some taxa likely never had pads (ancestrally padless) while other extant padless taxa had ancestors with pads that lost them (secondarily padless) (Gamble et al., 2012; Russell and Gamble, 2019). Ideally analyses would compare extant padded taxa (incipient, distal, or basal morphologies), to padless ancestral morphologies but it can be difficult to precisely estimate which clades within Gekkota represent independent origins (Gamble et al., 2012; Russell and Gamble, 2019). Ancestral state reconstructions in systems with directional evolution that lack fossil data can be difficult (Slater et al., 2012). In addition, we have very few direct observations of ancestral gecko morphologies. although there are examples of fossil geckos, including specimens in amber (Arnold and Poinar, 2008; Daza et al., 2014; Sherratt et al., 2015; Daza et al., 2016), there are very few prior to the diversification of extant families with information about soft tissue like toe scalation.

We started our investigation of parallel toe pad morphologies with the assumption that each of the four pad bearing families (Diplodactylidae, Sphaerodactylidae, Phyllodactylidae, and Gekkonidae) initially exhibited a padless morphology. Instead of trying to simulate ancestral morphologies, we treated extant padless morphologies as a stand-in for each family’s ancestral padless morphology. Dollo’s law of irreversibility suggests secondarily padless clades may not re-evolve padless morphologies identical to those of ancestrally padless geckos and other lizards (Russell and Gamble, 2019). To confirm if secondarily padless geckos external morphology is similar to ancestrally padless geckos, we examined the diversity across padless geckos (plus Sphenodon). We conducted a phylogenetic Procrustes ANOVA, correlating species mean toe pad shape for all padless species (20 species, including Sphenodon) with their family taxonomic rank, while accounting for the non-independence of species. We also conducted a species-level between-group PCA to identify major axes of morphological variation across padless species (Figure 5), again including a projection of a trimmed phylogeny (Tonini et al., 2016), an estimation of the variation captured by each PC axis, and projected toe shapes along each axis one standard deviation above and below the mean shape of each axis.

Confident that we can justify comparing extant padless taxa to extant pad bearing taxa, we conducted a series of trajectory analyses to compare shape differences across assigned groups (Collyer and Adams, 2013). Using nested categories (i.e., toe pad classes within families), our first set of comparisons considered shape between padless and pad bearing classes across families (Figure 6), testing if the differences between toe pad classes are consistent across families, pointing to parallel evolution. For example, our analyses quantified the distance (magnitude) and direction (angle) in morphospace between padless diplodactylid geckos and diplodactylid geckos with basal pads. This was repeated for each family, quantifying the distance and direction between padless and basal pad bearing taxa within each family. These family-specific morphological differences (magnitudes and angles) were then statistically compared across families to determine if the differences between padless and basal pads were significantly different in magnitude and/or angle. This was then repeated for each toe pad class (padless versus incipient and padless versus distal, Figure 6). Our second set of trajectory analyses considered comparisons within families, evaluating how each toe pad class differed morphologically from padless morphologies (Figure 7). These analyses again quantified the magnitude and angle. We compared padless versus basal, padless versus distal, and padless versus incipient in Diplodactylidae, Sphaerodactylidae, and Gekkonidae; whereas we only compared padless versus basal and padless versus distal in Phyllodactylidae, due to Phyllodactylidae not containing any examples of incipient pads.

All of our trajectory analyses relied on fitting linear models that used a randomized residual permutation procedure with 10,000 permutations; the “lm.rrpp” function from the RRPP R package (Collyer and Adams, 2018). For our first set of trajectories analyses, our linear models predicted toe pad shape using family, pad class, and their interaction as independent variables. Our second set of trajectory analyses predicted toe pad shape using pad class, an assigned pair designation, and their interaction as independent variables. Using the output from our linear models, we used the “trajectory.analysis” function to assign repeated nested factors (e.g. pad classes within families) and statistically compare the distance (magnitude) and trajectory (angle) in morphospace between toe pad classes across families. Over the course of the above tests, we made a total of 50 comparisons. In order to control for an expected number of false positives (Zelditch et al., 2012), we used the “p.adjust” function in R, implementing the false discovery rate method (Benjamini and Hochberg, 1995; Whitlock and Schluter, 2015) to calculate adjusted p-values, termed q-values.

Trait evolution

To investigate how toe pad shape may have evolved across geckos, we developed a new application of the node height test (Freckleton and Harvey, 2006), calculating standardized phylogenetic independent contrasts using high dimensional shape data (Felsenstein, 1985). A node height test uses standardized phylogenetic independent contrasts as estimates of the rate of morphological change at each node in a phylogeny. Our modified approach allows for the estimation of standardized phylogenetic independent contrasts using shape data instead of univariate morphological data. Traditional standardized phylogenetic independent contrasts are calculated by dividing the difference in trait values between sister taxa by the square root of their summed branch lengths. To estimate standardized phylogenetic independent contrasts via shape data, instead of calculating the difference in trait values, we estimated Procrustes distances between the shapes of sister taxa. Procrustes distance is the square root of the sum of squared differences in the positions of the landmarks in two shapes (Rohlf and Slice, 1990; Dryden and Mardia, 1998). We then treated our estimated Procrustes distance the same as a traditional independent contrasts approach would treat regular contrasts, dividing by the square root of the summed branch lengths. Following traditional phylogenetic independent contrasts methods, after calculating a node’s contrast, the relevant sister taxa were collapsed into a single branch extending from the existing node by the product of the two condensed branch lengths divided by their sum. Similar to traditional phylogenetic independent contrasts calculations, each condensed branch was assigned a new trait value as a weighted mean of the two previous sister taxa’s trait values, weighted by their respective branch lengths. In our case, we calculated weighted average shapes for each condensed branch by dividing each previous tip shape by their respective branch length, summing them, and then dividing by the sum of the inverse branch lengths, again mirroring traditional phylogenetic independent contrasts calculations.

Using the above-described modified approach, we calculated standardized phylogenetic independent shape contrasts using a trimmed version of the previously mentioned ultrametric binary phylogeny from Tonini et al. (2016), focusing solely on Gekkotan taxa (excluding Sphenodon). We then conducted a node height test, using a linear regression with natural log transformed standardized phylogenetic independent contrasts as our dependent variable and node heights above the root as our independent variable, asking if the slope of the resulting line was significantly different than zero and calculating a 95% confidence interval (Figure 8). By comparing standardized phylogenetic independent contrasts against the depth of each node in the phylogeny, a significant change in the rate of evolution can be detected. No change in contrasts through time imply a constant rate of evolution, consistent with Brownian motion. A decrease in contrasts closer to the present implies a slowdown in evolution, consistent with an early burst model of evolution, and lastly, an increase in the contrasts would suggest an. increase in the rate of evolution and a late burst model of evolution (Freckleton and Harvey, 2006). It is worth noting we chose to include taxa that maximized taxonomic and morphological diversity. Our dataset includes most of the examples we are aware of within-genus gains or losses of adhesive pads or changes in toe pad class. However, this decision limited our ability to thoroughly sample genera with similar morphologies. Our sampling scheme omits many within-genera nodes with little morphological change (i.e. slow rates of evolution), resulting in a conservative test for evaluating a model of early burst evolution.

Results

PCA and pad classification

To better understand the general variation across our toe pad shape dataset, we performed a species-level principal components analysis (Figure 4). Our first PC axis captured 75.2% of our dataset’s variation, with PC2 capturing an additional 12.7%, PC3 capturing 5.2%, and PC4 capturing 1.5%, for a total of 94.6% across our first four PC axes (see Supplementary Figure S3 for PC 3 and 4). The projected shapes along each axis suggest that PC1 and PC2 clearly captured our assigned toe pad classes. PC1 captured the degree of toe dilation and the location of the three biggest scales, with narrow toes and low small scales captured by low PC1 values versus exaggerated pads and scales having higher values, distinguishing between padless/incipient pads and basal/distal pads. PC2 captured the narrowness/roundness of the focal scales and how spread out and uniform they were, clearly distinguishing between the small or narrow scales/lamellae of many padless and basal species versus the enlarged and often rounded scales on many incipient and distally padded specimens. PC3 and PC4 appeared to focus on how notched or divided the three focal scales/lamellae were, with lower values of both PC3 and PC4 capturing more divided/notched scales (see Supplementary Figure S3 for PC 3 and 4). Our padless and incipient species along with many of our basal species do not have notched or divided scales, and hence they tended to appear towards the upper right of the plot whereas most distal species have divided pads and so these specimens filled the lower left of the plot. The included phylogenetic projection in the plot includes many parallel or converging phylogenetic branches, highlighting gecko’s abundance of repeated evolution.

While our assigned toe pad classes were clearly captured by PC1 and PC2, there were some areas of overlap including two Sphenodon punctatus specimens within the cloud of specimens with distal pads. The three largest scales on these specimens were more distal on the toe than most of the other padless species, though similarly sized scales cover Sphenodon’s entire toe. In addition, padless Coleonyx elegans also overlapped with specimens exhibiting distal pads, due to a pair of large non-adhesive scales flanking each side of their claw, resembling a distal pair of adhesive scales. There were also a incipient pad bearing specimens from Cyrtodactylus loriae, Quedenfeldtia moerens, and Naultinus elegans that appeared lower on PC2, grouping with other padless specimens, due to their narrow, yet still relatively exaggerated scales.

We next investigated how well our assigned toe pad classes captured patterns in our shape data. Our phylogenetic Procrustes ANOVA found our toe pad categories significantly predicted species mean shape (p = 0.001). Our groupPCA cross-validation analyses successfully classified specimens into the correct class 92% of the time. This analysis also compared the shape of specimens assigned to each class. Specimens assigned to the incipient class were significantly different in shape compared to specimens assigned to the basal class (p = 0.01). The incipient class vs. the distal class, padless vs. basal, and padless vs. distal classes were also all significantly different (p = 0.0001), whereas specimens in the padless and incipient classes were not significantly different (p = 0.13).

Using Ward’s 1963 clustering criterion (Murtagh and Legendre, 2014), we generated a dendrogram, clustering specimens based on their toe shape (Supplementary Figure S4). These results recapitulate our PCA findings (Figure 4). Specimens generally clustered by toe pad class, with the deepest clades capturing the same groupings seen in PC1 and PC2. Each of our four assigned toe pad classes clustered into nearly monomorphic clades. Padless and incipient clades grouped together into a larger clade and distal and basal pad clades formed a larger clade. While specimens of the same species often clustered together, specimens did not generally cluster by family. Each of the four largest clades, comprised primarily of a single toe pad class, included examples from each family. As also observed in our PCA (Figure 4), select Cyrtodactylus loriae, Quedenfeldtia moerens, and Naultinus elegans specimens grouped with padless taxa instead of incipient specimens and padless Sphenodon punctatus and Coleonyx elegans specimens fell in the distal pad clade. Interestingly, our Hemidactylus specimens grouped in the distal clade, with taxa that have additional proximal pads (e.g., Strophurus, Oedura, Amalosia, Afroedura). In our PCA analyses (Figure 4), when compared to all the other specimens assigned to the basal class, Hemidactylus specimens fell lowest on PC1 and highest on PC2, nearest to specimens in the distal class. This clustering makes sense given the paired toe pads of Hemidactylus (but lacking a prominent distal pair), and our landmarking scheme’s inclusion of only the three largest scales or pairs of scales.

Padless taxa comparisons

Our phylogenetic Procrustes ANOVA, correlating mean toe shape of padless taxa with family rank found that family assignment could not significantly predict toe shape (p = 0.8). This suggests ancestrally padless taxa (Sphenodon, Carphodactylidae, and Eublepharidae (Russell and Gamble, 2019)) and padless taxa from the four other pad bearing families, of which most are presumably secondarily padless, are not significantly morphologically different. Our between-species PCA considering padless taxa (Figure 5) captured 91% of the variation across padless species with only two PC axes (PC1: 84.2% PC2: 7.0%). PC1 captured how distal the three largest ventral scales were. Many of the specimens at the low end of PC1, having more distally chosen scales, were the same specimens we found clustering with distally padded taxa in our clustering analyses (Supplementary Figure S4, e.g. Sphenodon, Nactus, Coleonyx, and Goniurosaurus). PC2 distinguished between taxa with transverse ventral scales (high values) and granular ventral scales (low values). Specimens high on PC2, having transverse scales, were often the same specimens found to cluster with taxa exhibiting incipient pads (Supplementary Figure S4, e.g. Gonatodes, Phyllurus, Pristurus, Cyrtopodion).

Parallel morphologies

We conducted two sets of trajectory analyses, comparing the magnitude and angle between padless and pad bearing morphologies across families (Figure 6) and within families (Figure 7). To visualize our results, each comparison is plotted in the same between-group principal component morphospace, using family + toe pad class as our PC groups (Figures 6, 7). This representation of morphospace is very similar to our between-species PCA (Figure 4), capturing similar dimensions of variation and patterns across toe pad classes.

Our first set of analyses compared the trajectories of toe pad shape between families to determine if the trajectories (magnitude or angle) were statistically different. When considering the differences between padless and basal classes across families (Figure 6A), we found all of our magnitudes to be similar (q-values over 0.05) except for Phyllodactylidae versus Gekkonidae (q = 0.005) and Phyllodactylidae versus Sphaerodactylidae (q = 0.005) using adjusted p-values (i.e. q-values). The morphological differences between padless and basal morphologies were similar across families except when considering Phyllodactylidae due to padless Phyllodactylidae being at the far low end of PC1, resulting in large morphological differences between padless Phyllodactylidae and Phyllodactylidae geckos with basal pads. This difference was significantly larger than the differences between padless and basal Gekkonidae or the differences between padless and basal Sphaerodactylidae. In addition, none of the angles between padless and basal classes were significantly different across families (q-values over 0.05), illustrating parallel differences between padless and basal morphologies across families.

When comparing the differences between padless and distal morphologies across families (Figure 6B), we again found the magnitudes of our comparisons to be similar (q-values over 0.05) other than the comparison of Phyllodactylidae versus Gekkonidae (q = 0.025) and Phyllodactylidae versus Sphaerodactylidae (q = 0.025) which were both statistically different. These results again illustrate general similarity in the difference between padless and distal morphologies, although the distance between Phyllodactylidae basal and distal morphologies is significantly larger than distance between basal and distal morphologies in other families. In addition, our padless and distal angle comparisons were all similar (q-values over 0.05) again displaying parallel morphological differences across families. When comparing padless to incipient morphologies across families (Figure 6C), we found all angles and all magnitudes to be similar (q-values over 0.05), again illustrating parallel differences between padless and incipient pads across families.

Our second set of trajectory analysis compared differences between toe pad categories within families (Figure 7). We found the same relationships in each family. Differences from padless to basal and padless to distal were all statistically similar in magnitude (q = 1) but different in angle (q < 0.02 for all comparisons). Distal and basal pads are similarly diverged from padless morphologies but diverged in different directions. Conversely, trajectories from padless to incipient were always statistically different in both magnitude and angle when compared to other trajectories (q < 0.02 for all comparisons). Incipient pads are less diverged as compared to distal or basal pads, but also diverged in their own direction, notably not following similar trajectories as basal pads.

Trait evolution

Our node height test analyses asked if our shape data suggested an increase, decrease, or constant rate of evolution of toe shape across time by comparing shape-based standardized phylogenetic independent contrasts across time, treating these contrasts as an estimate of the rate of shape evolution at each node. Our analyses found a significant decrease in contrasts through time (p = 0.035) suggesting an early burst model of trait evolution (Figure 8). It is worth noting the prevalence of lightly colored branches 20 to 70 million years above the root, after the diversification of most families of geckos and during the diversification of extant genera, consistent with the many suspected within-family, yet between-genera transitions between pad bearing and padlessness (Russell and Gamble, 2019). Notable specific results include the dark node/branches connecting our two focal Sphaerodactylus species. Given the relatively deep age of the Sphaerodactylus genus and little morphological diversity within the genus, we estimated a very slow rate of evolution for these branches. We also observed three nodes with relatively fast rates of evolution. The node shared by our two Lucasium species, with one species having distal pads and the other being secondarily padless, suggested a high rate of evolution for the genus. We also observed a high rate of evolution for node leading to Homonota and Phyllodactylus, which again represents either a gain or loss of pads, depending on the reconstruction (Russell and Gamble, 2019). Lastly, we observed a notably high rate of morphological evolution for the node and branches connecting the clade containing Microgecko (padless) and Cnemaspis from the Indian subcontinent (incipient pads) to the rest of the Gekkonidae family. The weighted mean shape of the Cnemaspis/Microgecko clade was notably different in shape than the weighted average shape of the rest of Gekkota, although with such short branch lengths connecting these two clades, this different may have been subtle. In addition, different phylogenetic topologies are sometimes observed with these taxa (Russell and Gamble, 2019).

Discussion

Our results suggest that gecko toe pads have a history of repeated parallel evolution, with multiple examples of the same distinct morphologies (basal, incipient, and distal pads) evolving from presumably similar padless ancestors (Russell, 1979). In addition, our results suggest this repeated evolution occurred following an early burst pattern of evolution, slowing down over time. Below we discuss methodological considerations of our analyses, possible toe pad transitions and evolutionary processes responsible for our observed morphological patterns.

Methodological considerations

Our study included images from both live and preserved specimens. Specimen fixation and preservation is known to shrink lizard specimens up to 4% (Vervust et al., 2009; Maayan et al., 2022). If the preserved specimens included in our study shrunk isometrically, these changes would affect our results, due to the size independence of geometric morphometric analyses. However, this is likely not the case. Soft tissue is often more dramatically affected during preservation. Despite non-isometric shrinkage having likely occurred in our preserved specimens, these changes likely had minimal effect on our results. Our data quantified relative scale/lamellae shape and location on the toe, across very different morphologies. The variation introduced by combining live and preserved specimens is much less than the variation we see across genera. In addition, we do not expect fixation/preservation to change the relative shape and location of scales/lamellae on the toe, even if it does shrink the toe as a whole.

Toe straightness or curvature, an artifact of specimen preparation prior to collecting our images was not clearly captured by major PC axes, as it has been in previous work collecting similar data (Howell et al., 2022; Michaud et al., 2023). This is likely due to our use of a between-species groupPCA (and not a specimen-level PCA) and the variation in our dataset, which includes a wide variety of toe morphologies across families. By using a species-level PCA, our analysis emphasizes variation across species, whereas specimen-specific variation such as specific specimens being curved or straight, was treated as within-group noise. The wide breadth of morphological variation covered by our dataset may have also swamped any subtle sample-specific variation like toe curvature.

Our toe pad categorizations warrant further discussion. Geckos are often organized by their toe pad shapes into basal or distal categories (Gamble et al., 2012; Russell and Gamble, 2019; Riedel et al., 2021), with previous studies typically focusing on the overall location of the adhesive pad on the toe. Our toe pad classifications differed slightly, instead focusing the presence of an isolated distal pair of pads to designate a distal pad. This resulted in particular taxa being categorized differently than previous studies. Under our designations, we considered fan toed species (Uroplatus and Ptyodactylus) to have basal pads since they lack a pronounced distal pair of pads, even though the pad as a whole is concentrated at the distal end of their toe. Similarly, species like Rhoptropus, Rhoptropella, and Lygodactylus also have pads concentrated at the tip of their toes yet are considered to have basal pads by previous studies. This inconsistency is why we chose to focus on the presence/absence of an exaggerated distal pair of pads to assign toe pad class. In addition to our toe pad class assignments, our landmark assignments also include considerations. The landmarks capturing ventral scale morphology were chosen for their functional similarities and repeatability, not their homology, although some argument could be made for the distal most scales of pad bearing species being developmentally homologous (Griffing et al., 2022).

Prior to our trajectory analyses, we considered the morphological variation across our padless taxa to justify using them as representations of padless ancestral morphologies. The inclusion of Sphenodon puncatus, as an outgroup to all lizards, provided valuable information regarding ancestral padless morphologies. We found our focal padless species to be statistically similar in shape, with shape not correlated with family. The similarity of padless morphologies suggests we could view extant padless morphologies as a proxy for padless ancestral morphologies. It is worth noting that this approach did result in some caveats to our conclusions. Some of our extant padless taxa were likely secondarily padless (Russell and Gamble, 2019; Zhuang et al., 2019) and hence would represent transitions from padded to padless and not vice versa. In addition, such strong similarities across all padless taxa limited the possible outcomes of our trajectory analyses. Our toe pad trajectory comparisons across and within families were limited to divergent/anti-collinear or parallel/collinear geometric models (Bolnick et al., 2018), since all of our trajectories started at roughly the same place in morphospace; a padless morphology. Lastly, while our results allude to parallel evolution, i.e. parallel transitions from padless ancestors to derived pad bearing morphologies, our analyses were limited to evaluating parallel differences in morphospace (Bolnick et al., 2018; Stuart, 2019). We did not directly consider morphologies over time and cannot comment on other definitions of parallel evolution such as genetic parallelism, although there is evidence of repeatedly evolved developmental pattens across lizard toe pads (Griffing et al., 2022). However, as noted earlier, the generality of padless morphologies supports the idea that extant padless morphologies are externally similar to ancestrally padless morphologies. These results also suggests that external anatomy may violate Dollo’s law, although previous studies of internal anatomy highlight anatomical differences between primary and secondarily padless taxa (Russell and Gamble, 2019; Zhuang et al., 2019) supporting Dollo’s law is still valid at some levels during the revolution of padless morphologies. Studies considering of the revolution of whole digits in lizards similarly highlight the context dependence of Dollo’s law (Kohlsdorf and Wagner, 2006; Galis et al., 2010; Siler and Brown, 2011; Wagner et al., 2018).

Another consideration of our trajectory analyses is our prediction that the trajectories from padless to pad bearing morphologies would be parallel across families. In other words, morphological differences would have statistically similar angles through morphospace. The statistical tests we used as part of our trajectory analysis, which evaluate if our data has enough evidence to reject such a null, consider being similar as the null hypothesis. As a result, it is difficult to determine if our results present strong evidence for parallel differences or lack the power to reject our null hypothesis. This shortcoming was discussed in Bolnick et al. (2018). Despite these shortcomings, we believe that this is the best statistical test currently available, but improved statistical methodologies are needed.

Lastly, considering our node height test results, we found that the evolution of toe pad shape appears to have slowed down, with the rate of morphological change being slowest in our most recent nodes, most of which are genus-level nodes. In choosing which species to include in our study, we primarily focused on capturing morphological and taxonomic diversity across Gekkota. As a result, we included only a few morphologically similar species within a genus but always included morphologically divergent species within genera. For example, we included padless and pad bearing species within Lucasium, Gonatodes, and Cyrtodactylus, while alternatively, there are many other genera with little morphological variation, like Sphaerodactylus, of which we only included two species out of a possible ~100 described species. These decisions likely skewed our contrast data by including more recent young nodes with high rates of morphological evolution, resulting in a conservative test of early burst evolution. If we had included additional species, more thoroughly sampling morphologically consistent genera, this would have resulted in the addition of more young nodes with small contrasts (i.e. slow rates of evolution), and therefore an even stronger early burst pattern with a steeper negative slope. It is also worth noting that our phylogenetic independent contrasts were natural log transformed, compressing differences between large values and exaggerating differences between small values.

Toe pad transitions

While we did not directly model evolutionary transitions from padless to pad bearing or between toe pad classes, we can explore the topic. Our trajectory analyses only compared padless morphologies against other toe pad classes. We did not consider trajectories (i.e. transitions) between pad bearing classes, although previous authors have commented on such transitions (Russell and Gamble, 2019). We included fan toed geckos (Uroplatus and Ptyodactylus) within our basal category. Their toe pad morphology is similar to other basal pad bearing geckos, for example Gehyra and Hemiphyllodactylus, both with distally bifurcated or notched pads comprised of many narrow lamellae. This similarly could suggest a shared developmental pathways, or fan toed genera having evolved from taxa with similar basal pads. Inferring such evolutionary transitions is difficult because it relies on phylogenetic hypotheses. In Gamble et al. (2012) Uroplatus is sister to Ancylodactylus (i.e. African Cnemaspis that have incipient pads) and Ptenopus (padless), within a clade sister to Narudasia (padless), suggesting a weak shared ancestry between fan toes and incipient pad morphologies. This topology is similar to the phylogeny we used (Tonini et al., 2016). Yet in Russell and Gamble (2019), Uroplatus is sister to Afrogecko and Christinus, suggesting fan toes have a closer affinity to distal pad morphologies. Similarly, Ptyodactylus is sister to Thecadactylus (basal pads) in Tonini et al. (2016) and Gamble et al. (2012) but is sister to Haemodracon (distal pads) in Russell and Gamble (2019), again providing conflicting relationships between fan toed taxa and other pad classes. With only two extant fan toed genera, broad conclusions can be difficult to draw.

In addition to considering the origins of fan toed geckos, our study introduced the incipient toe pad class, which we suggest could represent an early stage in the evolutionary trajectory of basal pads. Riedel et al. (2021) describe a contemporary example of the origin of a distal toe pad morphology in Heteronotia, which complements the developmental patterns observed in Sphaerodactylus (Griffing et al., 2022) and may serve as a model for the evolution and development of geckos with distal pairs of pads. Alternatively, the pads of lizards with basal pads tend develop across a large area of the toe, with additional modifications later in development (Griffing et al., 2022). This pattern complements a hypothesis in which basal pads may evolve from an intermediate incipient pad like morphology, as observed in Gonatodes (Russell et al., 2015; Higham et al., 2016). While our trajectory analyses didn’t find incipient pads to be along the same trajectory as basal pads, we might not expect the process of evolution to be so linear. Ancestral traits need to be functional and provide a fitness advantage in and of themselves and so the evolution of basal pads, if hypothesized to move from a padless ancestor, through an incipient like morphology, to extant observed basal pad morphologies may not have followed a linear path through morphospace in order to maintain/maximize performance while following genetic lines of least resistance (McGlothlin et al., 2018; Rothier et al., 2022). We believe continued efforts to produce well supported species-level phylogenetic hypotheses while also incorporating embryological (Griffing et al., 2022) and eventually developmental genetic data will greatly help resolve these questions.

Adaptive origins

Our observed slowing parallel evolution could have been generated by selection/adaptation, stochastic processes, or constraints (e.g., developmental or genetic), although the adaptive value of adhesive toe pads appears intuitive, allowing geckos and other pad bearing lizards to navigate their habitat in unique ways. Toe pad morphology in lizards has been linked to performance, habitat use, and selective events in a variety of contexts (Losos, 1990; Irschick et al., 2006; Winchell et al., 2018; Donihue et al., 2020; Michaud et al., 2023). However, considering adhesive toe pads as key innovations has received mixed conclusions with toe pads facilitating arboreal niches (Miller and Stroud, 2021) and the anole adaptive radiation (Burress and Muñoz, 2021) but not higher evolutionary diversification rates in Australasian diplodactyloid geckos (Garcia-Porta and Ord, 2013).

Our observed early burst pattern of evolution could represent the filling of trait space as proposed by Simpson (1944), with each toe pad class representing an ecomorphological adaptive zone. This would also explain the repeated parallel evolution of similar morphologies, with different Gekkotan lineages moving towards the same adaptive peaks or zones. If this were the case, it would suggest that our morphological data may also fit a multi-regime Ornstein–Uhlenbeck (OU) model of trait evolution (Uyeda and Harmon, 2014), in which each toe pad class would be an optima and our observed slowdown in the rate of morphological evolution would be due to lineages approaching optima. To our knowledge, the only approach of fitting OU models to geometric morphometrics data is to limit analyses to specific principal component axes, fitting a PC axis to its own model of trait evolution. We imagine this approach could have limited value, as the analysis would only consider the aspects of morphology captured by a particular PC axis, which may or may not be biologically informative (Howell et al., 2022). While the idea of gecko toe pad classes representing discrete adaptive peaks is an intriguing concept, future studies considering the ecomorphology of toe pad classes will need to navigate the nuanced definitions of convergence/parallel evolution, considering repeatedly evolved similarities in biomechanics, performance, and habitat correlations of different toe pad morphologies. Such studies could strengthen future investigations of toe pads as key innovations (Miller et al., 2023). Even if repeatedly evolved toe pad morphologies appear similar externally and macroscopically, they would not necessarily be similar at the micro scale (i.e. setae) or internally and hence could operate differently (Russell, 1979; Russell and Gamble, 2019; Zhuang et al., 2019). Setal morphology and performance vary extensively across pad bearing lizards (Peattie, 2007; Hagey et al., 2014; Garner and Russell, 2021), between padless and pad bearing groups (Russell et al., 2015; Higham et al., 2016; Riedel et al., 2021), and with habitat use (Garner et al., 2022), but few studies have considered how setal morphology or micromechanics vary across toe pad classes.

To provide further challenges, the adaptive scenario that drove the repeated evolution of our observed extant morphologies may have been lost to time. Our results suggest that, early in the history of geckos, around 50-100 million years ago, toe pad shape evolved quickly (Figure 8), diversifying (and converging) into the morphological classes we see today, distinguishing closely related genera, and generating similarities between distantly related taxa. Then, more contemporarily, evolution appears to have slowed down. This complements our anecdotal observation of morphological variation within genera typically being slight variations of the same theme. If today’s observed toe pad classes diversified as Gekkotan families were diversifying, possibly acting as key innovations and opening a series of new adaptive zones, should we expect modern, extant morphologies to still correlate with those same dimensions of microhabitat that originally spurred their evolution? Those ecological correlations may have eroded over time, with deep radiations of toe pad morphologies now acting as constrains, with lineages losing their evolutionary lability, solidifying developmental pathways, and locking genera into specific morphologies. Meanwhile, environments continued to change, continents drifted, and lineages continued to diversify, slowing down the rate of toe pad shape evolution. If that is the case, we might expect contemporary within-genus ecological changes and diversification to have eroded any ecomorphological patterns connecting ancient toe pad class diversification to microhabitat. Instead, we may expect contemporary within-genus ecomorphological patterns. In this situation, species within each genus would have diversified by developing variations on the genus’ particular type of toe pad to thrive in different ecological contexts.

Additionally, our observed patterns may have also be affected by nonadaptive processes. The embryological development of gecko and anole toe pads appear to be strikingly similar (Griffing et al., 2021; Griffing et al., 2022), which is somewhat surprising for independently evolved morphologies. These developmental similarities allude to a nonadaptive explanation for the repeated evolution of Gekkotan toe pad classes, with similar morphologies possibly being the result of squamates’ tendency to modify their digit developmental pathways along the same “paths of least resistance” (McGlothlin et al., 2018; Rothier et al., 2022). This again suggests constraints may be an important factor in the diversification and parallelism of gecko toe pad evolution.

We strongly encourage the continued study of gecko toe pad macroevolution and believe geckos and other pad bearing lizards are an excellent system to study the interplay between ecomorphology, biomechanics, key innovations, adaptive radiations, and developmental constraints (Bauer, 2019). The rising fields of gecko toe pad development (Griffing et al., 2022) and genomics (Pinto et al., 2023) are currently laying the foundation for exciting future studies investigating the underlying factors contributing to our observed patterns of squamate toe pad diversity.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Ethics statement

Ethical approvals, various IACUC and wildlife collection permits, were received for our work including live animals in accordance with local legislation and institutional requirements.

Author contributions

JM: Data curation, Investigation, Methodology, Validation, Writing – original draft, Writing – review & editing. TH: Conceptualization, Formal Analysis, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This project was conducted using equipment funded by the Mississippi University for Women Science and Mathematics department and the Mississippi State Space Grant Consortium.

Acknowledgments

We wish to acknowledge the natural history museums that provided access to specimens for this project (American Museum of Natural History, Museum of Comparative Zoology, California Academy of Science, Queensland Museum, University of Kansas Biodiversity Institute, Museum of Vertebrate Zoology, Carnegie Museum of Natural History, Chicago Field Museum, Bernice Pauahi Bishop Museum, and the University of Michigan Museum of Zoology) and the people and organizations that facilitated our access to live specimens (Jon Boone, Eben Gering, Alyssa Stark, Kiwi Park Queenstown, Orana Wildlife Park, National Geographic Society and the Waitt Institute (#W216-12), BEACON Center for the Study of Evolution in Action, and the Mississippi University for Women). We would also like to thank Bailey Howell, Maggie Haines, Lauren Cirino, and four reviewers for comments on our manuscript and the Society of Integrative and Comparative Biology conference for the opportunity to share and discuss preliminary results with fellow scientists.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2023.1334870/full#supplementary-material

References