Gliding toward an understanding of the origin of flight in bats

Morphological data

We compiled linear measurements of 14 forelimb and 15 hindlimb traits (Fig. 1) from 231 extant mammals (1–2 specimens per species) across 21 orders (78% of mammalian orders) and 73 families (47% of mammalian families). Of these, 27 representative species of bats (48 specimens) were measured for the flyer group. These spanned both Yangochiroptera and Yinpterochiroptera suborders, five families, and fifteen genera ranging in mass from 5 to 700 g, thus broadly covering both the phylogenetic and physical diversity of Chiroptera. Linear measurements were obtained from specimens held at the Burke Museum of Natural History and Culture or from previously published work (Chen & Wilson, 2015; Grossnickle et al., 2020; Weaver & Grossnickle, 2020; Pevsner, Grossnickle & Luo, 2022). All specimens measured for this study were fully mature, determined by the closure of exoccipital–basioccipital and basisphenoid–basioccipital sutures on the cranium and ossification of all vertebrae and limb bones. Specimens were a mixture of male, female, and unknown sexes. We size-corrected all linear measurements by calculating log-shape ratios of each trait measurement (Mosimann, 1970) by log₁₀-transforming the mean trait value divided by the cube root of body mass prior to running all statistical analyses. Unfortunately, most specimens lacked accompanying body mass data; therefore, we compiled species-level body mass data from the PanTHERIA database (Jones et al., 2009) and the Handbook of the Mammals of the World book series (Wilson & Mittermeier, 2009).

We also collected the 29 measurements for four extinct bats (Onychonycteris finneyi, Icaronycteris index, Hassianycteris messelensis, and Tachypteron franzeni) using ImageJ (Schneider, Rasband & Eliceiri, 2012) on published full-skeleton images that include scale bars (Storch, Sigé & Habersetzer, 2002; Gunnell & Simmons, 2005; Simmons et al., 2008). We were only able to include these four extinct bats because we could only obtain a complete set of measurements from these species. O. finneyi, I. index, and H. messelensis have been recovered as stem bats in phylogenetic reconstructions and thus provide critical context for the evolution of powered flight in the absence of transitional bats in the fossil record (O’Leary et al., 2013; Rietbergen et al., 2023). Meanwhile, T. franzeni is generally thought to be a crown bat, so it is also interesting to include as the only non-stem fossil bat in our analyses (Storch, Sigé & Habersetzer, 2002; Ravel et al., 2016; Rietbergen et al., 2023). We estimated the body masses of extinct bats using a linear regression between body mass and the geometric mean of all linear measurements from the extant bats of our dataset, and used those body mass estimations to size-correct our fossil data again using log-shape ratios. We included fossil bats in our morphological dataset because, although extinct bats have a high degree of morphological similarity to extant bats (Norberg, 1985b; Bishop, 2008; Simmons et al., 2008), the fossil bat data may alter model estimations of the adaptive zone between flyers and gliders in an informative way. For example, the inclusion of fossils could make the flyer-glider trait optima more similar or more disparate by adding more phylogenetic history information and thus inform the understanding of the adaptive landscape underlying the evolution of mammalian flight.

To visualize the distribution of locomotor groups in the forelimb and hindlimb morphospaces, we performed principal component analyses (PCA) of the morphometric data using the princomp function in the R package stats (R Core Team, 2021).

Locomotor data

We categorized species into four locomotor groups: ground dwellers (n = 117), arborealists (n = 71), gliders (n = 16), and flyers (n = 27) following previous work (Chen & Wilson, 2015; Grossnickle et al., 2020; Weaver & Grossnickle, 2020; Pevsner, Grossnickle & Luo, 2022). The ground-dwelling mammals included diverse ambulatory, cursorial, saltatorial, and fossorial taxa. This ‘catch-all’ category was intentional as we are interested in investigating the arboreal to gliding to flying specialization trends only. Our arborealist category was also broad in that it includes both fully arboreal (tree-dwelling) and scansorial taxa. Despite these broad categories, ground-dwellers and arborealists show less variability in their trait optima for both the forelimb and hindlimb compared to gliders and flyers (see Results). Therefore, we do not view the ground-dwelling or arboreal optima as being pulled away from other groups in a forced way. Gliders included tree-dwelling mammals with derived morphologies (i.e., patagia and elongated limbs) for gliding locomotion (Grossnickle et al., 2020). Flyers included bats only.

Phylogenetic comparative analyses

We used evolutionary models (Hansen, 1997; Butler & King, 2004) to examine the influence of locomotor modes on the evolution of forelimb and hindlimb morphologies in our extant-only and extant+extinct datasets. We first tested the hypothesis that locomotor modes influenced the evolution of the forelimb and hindlimb using multivariate evolutionary models on our extant datasets (Butler & King, 2004; Hansen, 1997) in the R package mvMORPH (Clavel, Escarguel & Merceron, 2015). To reduce the dimensions of our data, we conducted a principal component analysis (PCA) on our 14-trait forelimb dataset and 15-trait hindlimb dataset. We retained the first three principal components (PCs) of our forelimb dataset (91.8% of the forelimb variation) and hindlimb dataset (85.3% of hindlimb variation). We then fit six multivariate models to the scores of these three PCs. We fit a multivariate Brownian motion model (mvBM), which models trait variance in a random walk over time; a single-regime multivariate Ornstein-Uhlenbeck model (mvOU1), in which scores of each PC evolve towards a single trait optimum (θ parameter of OU models), giving a single mean optimum for all locomotor groups; and multivariate OU models (mvOUM) with a given number of optima corresponding to selective regimes (based on locomotor modes). To test our hypotheses on overall limb evolution, we fit a four-peak OU model (mvOUM_loc4) for distinct optima among ground-dweller, arborealist, glider, and flyer regimes. To test alternative hypotheses, we fit three three-peak models (mvOUM_loc3a-c) with distinct combinations of the four locomotor groups. The mvOUM_loc3a model tests for distinct optima among ground-dweller, arborealist, and a combined glider-flyer regimes. The mvOUM_loc3b model tests for distinct optima among ground-dweller, combined arborealist-glider, and flyer regimes. The mvOUM_loc3c model tests for distinct optima among a combined ground-dweller-arborealist, glider, and flyer regimes.

To account for uncertainty in phylogenetic topology, branch lengths, and ancestral character states (Fig. 1), we used the make.simmaps function in the R package phytools (Revell, 2012) to simulate 10 stochastic character maps across 1,000 trees randomly drawn from the posterior distribution of “completed” trees from Upham, Esselstyn & Jetz (2019), resulting in 10,000 character maps of locomotory regimes. We randomly sampled 250 trees and fit all six models to these trees. We assessed relative model support using small sample corrected Akaike weights (AICcW). Models with ΔAICc scores less than two were considered supported and interpreted as evidence of locomotor groups influencing species’ adaptive zones. We acknowledge that using a subset of PCs instead of the full dataset could bias our results (Adams & Collyer, 2018; Uyeda, Caetano & Pennell, 2015). However, our results remained consistent when we increased the number of PCs in our model-fitting analyses to the first four or five PCs. Further, we used simulations to assess whether we had adequate power to accurately distinguish between complex mvOU models from Brownian motion (Boettiger, Coop & Ralph, 2012). We also calculated the phylogenetic half-lives (ln(2)/ $\alpha$) of the best fitting model to assess the responsiveness of locomotor regimes to adaptive peaks (Hansen, 1997).

If the mvOUM_loc4 model is best-supported, we interpreted this as each locomotor regime having its own distinct adaptive optimum. This result alone does not provide evidence for the gliding origin of bats. However, if gliders are located at an optimum between arborealists and flyers (i.e., if arborealist < glider < flyer optima or arborealist > glider > flyer optima), there is support for an evolutionary “pathway” between the two locomotor modes relative to an ancestral arborealist or ground-dwelling state (i.e., one could theoretically bridge adaptive peaks via these pathways). Therefore, if the optima gradually change in value (as opposed to random changes in value) from arborealists to gliders to bats, our hypothesis that bats evolved from gliding ancestors is supported.

If the mvOUM_loc3a model is best-supported, we interpret this as gliders+flyers having a shared adaptive optimum and arborealists and ground-dwellers having their own distinct evolutionary optima. This model provides the clearest support for flyers evolving from gliding ancestors because lineages of both groups are evolving towards the same shared adaptive optima. Therefore, if this model best fit the data, our hypothesis is supported.

If the mvOUM_loc3b model is best-supported, we interpret this as arborealists+gliders having a shared adaptive optimum and flyers and ground-dwellers having their own distinct evolutionary optima. There is phylogenetic justification for this grouping because all gliders evolved from ancestrally arboreal groups (e.g., gliding squirrels evolved from tree squirrels; see tree in Fig. 1). These groups are similar to the 4-group model in that if arborealists+gliders have an optimum between those of ground-dwellers and flyers, there is support for an evolutionary “pathway” between the two locomotor modes relative to an ancestral ground-dwelling state.

Finally, if mvOUM_loc3c model is best-supported, we interpret this as ground-dwellers+arborealists having a shared adaptive optimum and flyers and gliders having their own distinct evolutionary optima. As with the mvOUM_loc4 and mvOUM_loc3b models, if gliders have an optimum between those of ground-dwellers+arborealists and flyers, there is support for an adaptive “pathway” between the two locomotor modes relative to an ancestral arborealist or ground-dwelling state, providing evidence for our hypothesis that flyers evolved from gliding ancestors.

We also examined the influence of locomotor modes on the evolution of each individual forelimb and hindlimb trait in our extant-only and extant+extinct datasets using univariate evolutionary models. For each of the 29 traits, we fit a Brownian motion model (BM), a single-optimum Ornstein-Uhlenbeck model (OU1), and four multiple-peak OU (OUM) models. The OUM models included OUM_loc4, OUM_loc3a, OUM_loc3b, and OUM_loc3c, featuring the same locomotor regimes as the multivariate mvMORPH models. All models were fit using the R package OUwie (Beaulieu et al., 2012). To maintain consistency with the mvMORPH analyses, we applied a single selection parameter ( $\alpha$) and rate parameter ( $\sigma$²) across all selective regimes while allowing optima (θ) to vary between regimes (Hansen, 1997; Butler & King, 2004). For traits in which an OUM model was the best-supported model, we computed 95% confidence intervals using parametric bootstrapping of the model parameters with the OUwie.boot function in the R package OUwie (Beaulieu et al., 2012). We randomly sampled 100 of our character maps and performed 10 bootstrap replicates per stochastic character map sampled.

We repeated the above multivariate and univariate model fitting approaches on a merged dataset that includes extant species and the four fossil species. We first generated a MCC tree from the 1,000 trees randomly drawn from the posterior distribution of “completed” trees from Upham, Esselstyn & Jetz (2019) and manually placed the four extinct bats on the phylogeny as a polytomy with extant bats at 69.8 MYA, which was the average branch age of bats based on the posterior distribution of “completed” trees. Although Tachypteron franzeni is often hypothesized to fall within crown bats (Storch, Sigé & Habersetzer, 2002; Ravel et al., 2016; Rietbergen et al., 2023), the exact relationship with extant bats remains unresolved (e.g., Jones, Beard & Simmons, 2024). Therefore, we conservatively include T. franzeni in our polytomy instead of within extant bats. Regardless, the position of T. franzeni is unlikely to change the model results because all extinct taxa show fairly derived traits (compared to gliders, arborealists, and generalist species included) (e.g., T. franzeni lies consistently with both extant and extinct bats in morphospace and Fig. S1). We then simulated 50 stochastic character maps on the phylogeny and fit all four models across the 50 stochastically mapped trees to account for uncertainty in the ancestral character states.

Lastly, additional processes outside the scope of our investigation among arboreal, gliding, and flying transitions may also influence limb evolution, resulting in evolutionary shifts in limb morphology that are not associated with the evolution of gliding or flying. These hidden processes are not captured by our hypothesis-testing framework described above. Therefore, we also used PhylogeneticEM to determine if additional shifts of limb evolution occur along branches in a pattern not predicted by our locomotor categorical scheme. The PhylogeneticEM approach does not use a priori regime assignments (i.e., locomotor mode assignments are not incorporated into analyses). We used a scalar OU model that infers the full evolutionary rate matrix and accounts for correlations within multivariate datasets (i.e., PC1–PC3 for both forelimb and hindlimb datasets) using the R package PhylogeneticEM version 1.4.0 (Bastide et al. 2018).

Forelimb morphospace and evolutionary model-fitting

PC1 explains 81.5% of the forelimb variation and separates flyers (both extant and extinct) from the remaining three locomotor modes (Fig. 2A; Table S1A). PC1 also slightly separates gliders, particularly dermopterans, from arborealists and ground-dwellers. This axis describes the relative lengths of the armwing (humerus, radius, and ulna) and the handwing (metacarpals and phalanges). PC1 varies across relatively shorter (−PC1) to elongate (+PC1) bones, revealing a trend within the dataset in the elongation of multiple forelimb bone elements from arborealists to gliders to flyers. PC2 explains 7.8% of the morphological variation and primarily describes the relative widths of the long bones and hands. This variation of wide to narrow bones mostly occurs within locomotor groups, especially within the diverse ground-dweller group. PC3 explains 2.9% of the morphological variation with +PC3 dominated by metacarpal length, ulna olecranon length, and scapula height and −PC3 dominated by intermediate phalanx length and humerus distal width. Fossil bats overlap with extant bats in morphospace, which is consistent with the lack of transitional forms in the bat fossil record.

The mvOUM_loc3b model, which assigns arborealists and gliders to the same regime, was the best-fitting model for the extant multivariate forelimb dataset (PCs 1–3) (AICcW = 0.99, Table 1A) with an average phylogenetic half-life of 15.2 million years (myr) among PCs 1–3. Distinct optima were supported for these three locomotor groups, with PC1 and PC2 optima scores increasing in value from ground-dweller to arborealist+glider to flyer hindlimbs (Fig. 2C). Repeating the multivariate model-fitting with the addition of the four extinct flyer species yielded a mvOUM_loc4 model as the best fit (mAICcW = 0.99; Table S6A) with an average phylogenetic half-life of 127.2 myr among PCs 1–3. However, the optima exhibited unrealistic magnitudes (Table S6B).

When we examined each forelimb trait individually using univariate models, we found that the OUM_loc4 model was the best-fitting model for six of the 14 forelimb traits (AICcW = 0.51–1.00), OUM_loc3b was the best-fitting model for five traits (AICcW = 0.36–0.65), OUM_loc3c was the best fitting model for metacarpal length (AICcW = 0.63), and OUM_loc3a was the best-fitting model for ulna olecranon length (AICcW = 0.58; Table S2A). The OU1 model was the best fit for the proximal phalanx width of the third manual digit (AICcW = 0.51). Armwing (i.e., humerus and radius) length traits have optima (from the best-fitting OUM_loc4 models) that gradually increase from ground-dwellers to arborealists to gliders to flyers; parametric bootstrapping revealed little to no overlap in the 95% confidence intervals of optima between most locomotor groups (Figs. 3B and 3C). Optima of armwing widths (from the best-fitting OUM_loc3b models) indicated that flyers exhibit relatively wider humeri at the distal end, proximal end, and mid-shaft compared to all other locomotor groups (Figs. 3G–3I; Table S3A). Flyers exhibit relatively higher optima for handwing length traits (i.e., longer third manual metacarpal, proximal phalanx, and intermediate phalanx) than all other locomotor groups (Figs. 3D–3F; Table S3A). For each, glider optima were intermediate between arborealist and flyer optima. For handwing width traits, flyers exhibit relatively higher OUM_loc3b optima for metacarpal widths compared to all other locomotor groups (Fig. 3J), but no differences in proximal phalanx width as evidenced by OU1 support for this trait. An adaptive optimum for a flyer-glider combined group in the OUM_loc3a models was supported for ulna olecranon length only (Table S2A). In total, however, glider optima (in combination with arborealists) were between the 95% confidence intervals of ground-dweller and flyer optima for all traits except scapula length/width, humerus proximal width, and metacarpal width (Figs. 3A, 3G, 3J; Table S3A). For these last trait optima, glider+arborealist optima were statistically the same as that of ground-dwellers while flyers exhibited higher optima than all other locomotor groups (Figs. 3A, 3G, 3J; Table S3A).

Adding fossil flyers into forelimb model-fitting analyses yielded overall slightly more support for OUM_loc3a models compared to extant-only modeling; scapula height and humerus proximal width demonstrate support for a flyer-glider combined group (Table S4A). OUM_loc4 and OUM_loc3b models were still widely supported as the best-fitting model for six traits and five traits, respectively; OU1 was still the best-fitting model for proximal phalanx width (Table S4A). Overall, the length and width optima values for the forelimb were very similar between extant+extinct and extant only-modeling (Table S5A).

Hindlimb morphospace and evolutionary model-fitting

PC1 explains 47.5% of hindlimb variation and more positive PC1 scores describe the relative elongation (+PC1) of the fibula, metatarsal, and pedal proximal and intermediate phalanges (Fig. 2B; Table S1B). PC1 slightly separates gliders and flyers, with gliders falling towards +PC1 and flyers falling towards −PC1, indicating glider hindlimb autopodia are overall more elongated than flyers’. Arborealists and ground-dwellers are spread across PC1. PC2 explains 17.8% of the variation, where +PC2 is associated with relatively lengthened pedal proximal phalanges, intermediate phalanges, and ischium and −PC2 with relatively narrower pedal proximal and intermediate phalanges. PC2 loosely separated ground-dwellers, arborealists, and the glider-flyer group. PC3 explains 11.3% of the morphological variation with +PC3 dominated by ilium length, femur length, and fibula length and −PC3 dominated by metatarsal width, intermediate phalanx width, and proximal phalanx width. These loadings show that variation between the four locomotor groups mainly occurs throughout the hindlimb autopodia region. It should be noted that many bats have incomplete fibulae, so perhaps the reason that fibula length had a higher contribution than tibia length to the differences across PCs 1–3 was due to this anatomical difference in bats; regardless, since fibula and tibia length loadings were closely correlated, it appears that in this PCA the zeugopodal length of mammals also contributes substantially to overall morphological differences out of the traits measured (Table S1).

The mvOUM_loc4 model was the best-fitting model for the extant multivariate hindlimb dataset (PCs 1–3) (AICcW = 0.92, Table 1B) with an average phylogenetic half-life of 0.05 myr among PCs 1–3. Distinct optima were supported for the four locomotor regimes, with PC1 and PC2 optima scores decreasing in value from ground-dweller to arborealist to glider to flyer hindlimbs, albeit with overlap (Fig. 2D). Adding in the four extinct flyer species yielded again the mvOUM_loc4 model as the best fit, and very little change in the optima values (Table S6) with an average phylogenetic half-life of 15.2 myr among PCs 1–3.

When we examined each hindlimb trait individually using univariate models, we found that the OUM_loc4 model was the best-fitting model for three of 15 traits, all of which were hindlimb long bone lengths (AICcW = 0.33–0.60; Table S2B). Specifically, gliders have relatively longer tibia and femoral optima compared to flyers and the other locomotor groups (Figs. 4B and 4C, Table S3B ). Fibula length was the third hindlimb trait supported by OUM_loc4, but fibulae are rudimentary in most bats. The OUM_loc3a model best fit three traits (AICcW = 0.36–0.52), all of which were widths: femur shaft width, metatarsal width, and pedal proximal phalanx width (Table S3B). Parametric bootstrapping of these OUM_loc3a-supported models revealed that the glider-flyer group had relatively narrower femur shafts, metatarsals, and pedal proximal phalanges optima than the other locomotor groups (Figs. 4F, 4I, 4J; Table S3B). The OUM_loc3b model was the best-fitting model for four of the 15 traits (AICcW = 0.44–0.71): ilium, ischium, metatarsal, and intermediate proximal phalanx length (Table S2B). Flyers had significantly smaller relative ilium and ischium length optima than other locomotor groups. Gliders+arborealists had metatarsal and intermediate proximal phalanx optima intermediate between flyers and ground-dwellers (Figs. 4D and 4E; Table S3B). The OUM_loc3c model best fit three traits (AICcW = 0.42–0.48): pelvis length, femur distal width, and tibia proximal width (Table S2B). In the latter two traits, gliders exhibited intermediate width optima between flyers (which were smaller) and ground-dwellers+arborealists (which were larger). For pelvis length, like ilium and ischium length, there was no clear trend in the direction of the optima values across locomotor groups (Table S3B).

Adding fossil flyers into hindlimb model-fitting analyses yielded only one difference compared to extant-only modeling: pelvis length was now best fit with a OUM_loc3b instead of OUM_loc3c model (Table S4B). The length and width optima for the hindlimb were very similar between extant+extinct and extant-only modeling (Table S5B).

Evolutionary shifts in limb morphology

The PhylogeneticEM model indicated seven evolutionary shifts in the forelimb and 12 shifts in the hindlimb (Fig. 5). These shifts occur along the branches of named clades. In the forelimb, evolutionary shifts occurred along monotremes, three clades of “moles” (i.e., Notorycitidae, Chrysochloridae, and Talpidae), Chiroptera, Antilopinae, and Primatomorpha. In the hindlimb, evolutionary shifts occurred along Cingulata, Tenrecidae, Chiroptera, Primatomorpha, Hystricognathi rodents, Dipodidae, and marsupials. Further shifts in hindlimb evolution occur within marsupials including Notoryctidae, Peramelidae, Myrmecobiidae, Potoroidae, and Macropodidae.

Forelimb trait evolution

We found evidence that mammalian glider forelimb morphologies are evolving toward adaptive peaks that are intermediate between those of arborealists and flyers. Our multivariate forelimb results and univariate length trait results provide support for the hypothesis that selective pressures for flight influenced the evolution of flyer forelimbs from glider-like forelimbs. Using multivariate evolutionary models, we found a three optima model with distinct selective regimes for ground-dwellers, arborealists+gliders, and flyers as the best fit to the first three PCs of our data (Table 1A). The arborealist+glider optima are intermediate between those of flyers and ground-dwellers (Fig. 2C). PC1 is mostly explained by the variation in the forelimb bone lengths and PC2 by their widths, creating a coarse picture that forelimb bones likely become relatively longer and possibly wider from ground-dwellers to arborealists+gliders to flyers. Results from our univariate evolutionary models are consistent with these findings, revealing that glider forelimb optima are intermediate between flyer and arborealist optima for length traits (Figs. 3A–3F), but only partially for width traits (Figs. 3G–3J) across both the armwing and handwing (See Table 2 for a breakdown of whether each univariate model supports our hypothesis). The phylogenetic half-lives of the best multivariate model (15.2 myr) and mean of best univariate models of each trait (18.4 myr) are short relative to Chiroptera (ca. 70 MYA), providing support that the forelimb traits are strongly pulled toward distinct locomotor peaks across the adaptive landscape. Furthermore, variation in phylogenetic half-lives from univariate models (Table 2) suggest that these traits are pulled towards distinct peaks either over different adaptive slopes or at different rates.

Flyers exhibit trait optima indicating evolution of relatively longer forelimb bones compared to all other locomotor groups (Figs. 2A, 3A–3F) particularly in the handwing, which exhibited rapid evolution as a necessary adaptation for bats (Vaughan, 1970; Sears et al., 2006; Amador, Simmons & Giannini, 2019). Increased bone lengths may enhance the surface area of the patagium for greater lift (Norberg, 1985b). Flyers also displayed evidence of evolution toward relatively wider bones in the armwing compared to all other locomotor groups (Figs. 3G–3J) consistent with findings that bat humeri and radii exhibit significantly larger diameters compared to nonvolant mammals (Swartz & Middleton, 2008). These relatively wider proximal forelimb bones combined with thinner yet highly mineralized walls provide resistance to bending and torsion during wing flapping (Swartz, Bennett & Carrier, 1992; Papadimitriou, Swartz & Kunz, 1996). Furthermore, the transition from relatively more robust proximal ends of the humerus to relatively less robust distal ends may also reflect greater muscle mass closer to the center of rotation of the wing joints, which would greatly decrease flight inertial power and facilitate economical flight (Norberg, 1985a; Swartz, 1997). In the handwing, we found flyers exhibit relatively thicker third metacarpals (Fig. 3J) that corresponded with its relatively increased length (Fig. 3D) compared to nonvolant mammals. In contrast, the proximal phalanx length of the third manual digit was relatively longer compared to nonvolant mammals whereas its width was best supported with an OU1 model, consistent with previous work showing that the proximal phalanx diameter is not significantly smaller in bats compared to nonvolant mammals (Swartz & Middleton, 2008). As a result, bats exhibit relatively thick bone cortices to combat greater bending moments associated with this relatively increased elongation (Swartz, Bennett & Carrier, 1992; Papadimitriou, Swartz & Kunz, 1996; Swartz, 1997). Increased cortical thickness along with the reduction of mineralization of the handwing work together to resist bending failure while simultaneously maintaining flexibility during flying (Papadimitriou, Swartz & Kunz, 1996).

In gliders, we found evidence that forelimb bone length optima are largely intermediate between ground-dwellers and flyers. Gliders tend to exhibit relatively greater length optima in most forelimb bones compared to ground-dwelling and arboreal mammals (Thorington & Heaney, 1981; Thorington & Santana, 2007; Grossnickle et al., 2020; Linden et al., 2023; Rickman et al., 2023) but do not reach the lengths of flyer forelimb bones. Nevertheless, like flyers, the longer bones expand the surface area of the patagium (Thorington & Heaney, 1981; Thorington & Santana, 2007) and enable further adjustment of their patagium during flight, leading to more efficient aerodynamics, increased agility, and maneuverability (Thorington, Darrow & Anderson, 1998; Zhao et al., 2019). Forelimb bone widths in gliders, while clearly intermediate between those of ground-dwellers and flyers for some traits (Figs. 3H, 3I) share an optimum with arborealists (shown by the OUM_loc3b support for these traits; Figs. 3G–3J) and overlap with ground-dwellers for other trait optima (Figs. 3G, 3J). Together, these results suggest gliders may maintain high relative widths while continuing to elongate forelimb bones, supporting previous findings that the evolutionary elongation of glider forelimbs is not also accompanied by slimming of the bones (Grossnickle et al., 2020).

Hindlimb trait evolution

We found mixed evidence that mammal glider hindlimb morphologies are intermediate to those of arborealists and flyers. Using our multivariate modeling approach, we found distinct optima for ground-dwellers, arborealists, gliders, and flyers that best fit the first three PCs of our hindlimb data (Table 1B). The trait loadings indicate that the autopodial bones dictate the distribution of PC scores in PCs 1 and 2 (Table S1B). However, there is much overlap among the optima of ground-dwellers, arborealists, and gliders in morphospace (Fig. 2D).

The univariate models uncovered further patterns: 1) width optima of hindlimb bones were consistently narrower in flyers and gliders compared to arboreal mammals (Figs. 4F–4J), and 2) there was no pattern from arboreal to glider to flyer in the hindlimb length optima. Univariate modeling revealed that the widths of hindlimb bones tend not to differ between flyers and gliders (Figs. 4F–4J) because the widths of hindlimb bones in flyers and gliders exhibited either a single optimum (i.e., OUM_loc3a) separate from ground-dwelling and arboreal mammals (Figs. 4F, 4I, 4J) or an overlap in confidence intervals of optima between flyers and gliders (Figs. 4G, 4H). This shared flyer-glider adaptive zone suggests that flyers and gliders exhibit narrower hindlimb optima, especially in the autopod region (Figs. 4F–4J). In contrast, univariate modeling revealed distinct optima in long bone lengths of the hindlimb among the four locomotor regimes, in which gliders exhibited relatively longer femoral and tibial optimal lengths compared to flyers, arborealists, and ground-dwellers (Figs. 4B and 4C). Unlike bat forelimbs, bat long bones in the hindlimb (i.e., the femur and tibia) are not the longest compared to the other locomotor regimes; in fact, the confidence intervals of optimal lengths largely overlap between flyers and arborealists (Table S3B). Furthermore, flyers exhibit relatively shorter metatarsals but relatively longer pedal intermediate phalanges compared to all other locomotor regimes, whereas these autopodial lengths do not differ between gliders and arborealists (Figs. 4D and 4E). The different evolutionary patterns of hindlimb bones demonstrate the unique locomotor diversity among the four regimes. Although gliders and flyers are both aerial, these behaviors require different functional morphologies. In gliders, relatively longer hindlimbs are needed to increase the surface area of the patagium in tandem with relatively longer forelimbs to generate lift during gliding (Thorington & Heaney, 1981). In contrast, in bats, the relative lengthening of the forelimb supports most of the wing membranes whereas the hindlimb, specifically the ankle, serves as an anchor to the wing membrane. Bat hindlimbs are still used for flight by helping modulate lift, drag, or pitch (Cheney et al., 2014).

Overall, gliders appear to be evolving relatively lengthened hindlimb bones while bats are evolving more gracile hindlimb bones. Considering these results, we reject our prediction that flyer and glider hindlimbs overlap in their adaptive optima. This finding may point to the fact that modern gliders are completely distinct evolutionarily and potentially even morphologically from the gliding ancestors of bats. These results may highlight the pitfalls of modeling modern gliders as proxies for ancestral ones instead of invalidating our hypothesis that bats evolved from ancestral gliders. Further support of these hypotheses is provided by the relatively short phylogenetic half-lives of the best multivariate model (0.05 myr) and average of the best univariate models of each trait (15.3 myr) compared to the age of Chiroptera (ca. 70 MYA).

Incorporating fossil bats

The morphospace positions of fossil bats in the forelimb PCA confirm extinct bats’ morphological similarity to extant bats in the degree of forelimb specialization (Fig. 2A). Amador, Simmons & Giannini (2019) modeled aerodynamic traits instead of skeletal traits and found that the Eocene bat Onchonycteris finneyi occupies an area of aerofoil morphospace between extant gliders and bats. In our morphospace, O. finneyi is the closest fossil bat to gliders across PC1 (Fig. S1A), congruent with Amador, Simmons & Giannini’s (2019) findings. Tachypteron franzeni is the only extinct species included that is hypothesized to be part of an extant bat family (Emballonuridae or Miniopteridae), but its forelimb morphology lies outside the range of living bats along PC2 (Fig. S1A) (Storch, Sigé & Habersetzer, 2002). A caveat is that no species of either crown Emballonuridae or Miniopteridae were included in this study.

Extinct flyers occupy the same region of hindlimb morphospace as extant flyers, but three of four species overlap clearly with the region of morphospace of non-flyers, a higher proportion compared to that of extant bats (Fig. 2B). In this hindlimb morphospace, T. franzeni is the extinct species that most clearly overlaps with extant flyers (Fig. S1B). Based on its intermediate limb traits, Simmons et al. (2008) suggested that O. finneyi was likely more proficient at climbing than extant bats and may have employed quadrupedal locomotion and under-branch hanging. Interestingly, O. finneyi and Desmodus rotundus (the common vampire bat) are almost indistinguishable in the hindlimb morphospace (Fig. S1B). Vampire bats are known for their adept quadrupedal locomotion and climbing (e.g., Riskin, Bertram & Hermanson, 2005), so O. finneyi’s similar morphology may lend support to more walking-capable early bats. However, the linear measurements we used here may not reflect the complex adaptations required for these behaviors (e.g., ninety-degree joint rotations for under-branch hanging, Riskin, Bertram & Hermanson, 2005), so the placement of these species in morphospace leaves many possibilities open for ancestral locomotor behaviors.

Including extinct flyers in our model-fitting analyses as the closest morphologies we have to flyer-glider intermediaries did not significantly alter estimations of the adaptive zone between flyers and gliders; the inclusion of extinct species appeared to make the flyer-glider trait optima marginally more similar for the forelimbs and had no change for the hindlimbs (Table 2). The addition of extinct species to the multivariate model-fitting analysis shifted support from mvOUM_loc3b to mvOUM_loc4, suggesting that the glider optimum is potentially distancing itself from arborealists’ and towards flyers’ (Table S6A). Unfortunately, the extant+extinct mvOUM_loc4 model predicted extremely negative optima values, which represent unrealistic morphologies that are potentially artifactual, making it difficult to gain more insight into this result (Table S6B). The addition of extinct species to univariate forelimb model-fitting analyses yielded overall slightly more support for OUM_loc3a models compared to extant-only modeling (Table S4A). The addition of extinct species to hindlimb model-fitting analyses yielded overall no change in support for the mvOUM_loc4 model or OUM_loc3a models compared to extant-only modeling (Tables S6, S4B). Regardless, the length and width trends are the same as above between locomotor groups across the fore- and hindlimb (Table S5).

Limb trait optima and the adaptive landscape

Our results provide insight on the potential adaptive landscape connecting arboreal, gliding, and flying locomotion (Fig. 6). We found that glider forelimb optima fall intermediate between those of arborealists and flyers in their lengths but only partially in their widths. Except for long bone lengths, glider hindlimb optima are intermediate between ground-dweller/arborealist and flyer optima or even overlap with flyer optima. Together, these results emphasize that selection on key length/width traits may be pulling glider traits towards a flying adaptive zone yet that certain glider traits (e.g., femur and tibia lengths) are derived specifically for gliding morphology. Bishop (2008) proposed that the high degree of convergence among gliders suggests functional constraints on the evolution of powered flight from gliding; however, more recent work found that gliding lineages display incomplete convergence (species retain unique morphologies even while evolving some similarities; Leal, Knox & Losos, 2002; Stayton, 2006). For gliders, there is little evidence that gliding is under strong selective pressures toward a single-optimum gliding morphology (Grossnickle et al., 2020). The clustering of bat species in forelimb morphospace (Fig. 2A) may represent a narrow adaptive peak within which forelimb morphologies are highly conserved. In contrast, the broad glider morphospace may represent a gradually inclined selection plane, with gliders still far from a peak (Grossnickle et al., 2020). This broad adaptive slope could represent the foothill of a steep flyer peak such that gliders are very slowly evolving towards a flyer peak (Fig. 6A). The evolutionarily oldest extant glider lineages, dermopterans (colugos) and anomalurids (scaly-tailed squirrels), have evolved farthest from ancestral arborealists (Grossnickle et al., 2020) although some extinct gliders may have experienced relatively greater change from ancestral arborealists (e.g., Meng et al., 2017; Luo et al., 2017). Here, the dermopteran species Cynocephalus volans and Galeopterus variegatus are the gliders closest to the flyer forelimb region of morphospace (Figs. 2A; S1A), suggesting a trajectory in morphospace from gliding to flying. Further, colugos are known as the “mitten gliders” because of their interdigital membranes, which, along with their common carpal morphology with bats (Simmons & Geisler, 1998), reflects hypotheses that interdigital webbing was key to evolving powered flight (Anderson & Ruxton, 2020; but see Gardner & Dececchi, 2022 for an alternate view).

Our results also suggest that both flyers and gliders exhibit broad adaptive optima (i.e., having large confidence intervals) in the hindlimb due to gliders sharing an optimum or overlapping optima with flyers, or being intermediate between non-flyers and flyers for many traits (Figs. 2D, 4; Amador, Simmons & Giannini, 2019). These broad adaptive optima could represent overlapping peaks for gliders and flyers separated only by a shallow fitness valley driven by the gliders’ elongated femur and tibia lengths (Fig. 6B). The hypothetical adaptive landscape is therefore valuable because the transition from quadrupedalism to hindlimb suspension has been hypothesized to enable forelimbs to evolve more freely since they would no longer be weight-bearing structures (Bishop, 2008). Indeed, there is both quadrupedal suspension and hindlimb suspension in dermopterans, resulting in a tradeoff in the ability to walk quadrupedally upright (Fujiwara, Endo & Hutchinson, 2011; Granatosky, 2018). In the PC1 and PC2 hindlimb morphospace, G. variegatus is the only glider that overlaps with the flyer morphologies (Figs. 2B; S1B), reflecting a hindlimb adaptive optimum potentially more easily attained than that of flyers’ specialized forelimbs (Fig. 6).

Our modeling approach offers a first step towards support for the hypothesis that selection may more strongly pull certain limb traits towards a flyer adaptive zone. In particular, gliders’ intermediate forelimb length and overlapping hindlimb width optima between that of arborealists and flyers could be the traits driving lineages across fitness valleys due to their greater effect in, for example, the functionality of the gliding/flying apparatus. Short phylogenetic half-lives relative to the crown age of Mammalia provide additional evidence of distinct locomotor peaks. However, to complete the tentative adaptive landscape presented here (Fig. 6), these results could be strengthened by more comprehensive coverage of the lengths/widths across more mammal species and biomechanical studies tying these traits to specific functions and selective pressures.

Importantly, we acknowledge that the evolutionary transitions from arboreal to gliding to flying are not the only processes that dictate the adaptive landscape of mammalian limb evolution. Additional processes not explicitly tested here may also influence the evolution of the forelimb and hindlimb. Using PhylogeneticEM, we found that evolutionary shifts in limb morphology occur along the branches of named clades (Fig. 5). In the forelimb, we find evidence that these shifts are loosely structured based on locomotor mode. Antelopes (Antilopinae), colugos and primates (Primatomorpha), and bats (Chiroptera) exhibited evolutionary shifts in the forelimb, and all exhibit adaptations to cursorial, arboreal, and flying locomotor modes, respectively (Polly, 2007). Furthermore, three of the shifts occur on branches along three clades (marsupial moles, Notoryctidae; golden moles, Chrysochloridae; and moles, Talpidae) that independently evolved subterranean fossoriality. Their adaptations for stout, robust forelimbs are well documented (Polly, 2007). Evolutionary shifts in the hindlimb can also be partially explained by distinct locomotor modes including arboreal, flying, and saltatorial adaptations in colugos and primates (Primatomorpha), bats (Chiroptera), and jerboas (Dipodidae). Five shifts in hindlimb morphology also occur within marsupials, a result that may be unsurprising considering that marsupials display greater hindlimb disparity compared to placentals (Pevsner, Grossnickle & Luo, 2022). Overall, these results highlight that confounding effects from ecological adaptations and phylogenetic structure can result in clade-specific shifts across the adaptive landscape of mammals and other vertebrates (Uyeda & Harmon, 2014; Felice et al., 2019; Law, 2021; Law et al., 2022; Law, Hlusko & Tseng, 2024).