Qualitative skeletal correlates of wing shape in extant birds (Aves: Neoaves)

Abstract

Background

Among living fliers (birds, bats, and insects), birds display relatively high aspect ratios, a dimensionless shape variable that distinguishes long and narrow vs. short and broad wings. Increasing aspect ratio results in a functional tradeoff between low induced drag (efficient cruise) and increased wing inertia (difficult takeoff). Given the wide scope of its functional effects, the pattern of aspect ratio evolution is an important factor that contributes to the substantial ecological and phylogenetic diversity of living birds. However, because the feathers that define the wingtip (and hence wingspan and aspect ratio) often do not fossilize, resolution in the pattern of avian wing shape evolution is obscured by missing information. Here I use a comparative approach to investigate the relationship between skeletal proxies of flight feather attachment and wing shape.

Results

An accessory lobe of the internal index process of digit II-1, a bony correlate of distal primary attachment, shows weak but statistically significant relationships to aspect ratio and mass independent of other skeletal morphology. The dorsal phalangeal fossae of digit II-1, which house distal primaries VIII and IX, also show a trend of increased prominence with higher aspect ratio. Quill knobs on the ulna are examined concurrently, but do not show consistent signal with respect to wing shape.

Conclusions

Although quill knobs are cited as skeletal correlates of flight performance in birds, their relationship to wing shape is inconsistent among extant taxa, and may reflect diverging selection pressures acting on a conserved architecture. In contrast, correlates of distal primary feather attachment on the major digit show convergent responses to increasing aspect ratio. In light of the diversity of musculoskeletal and integumentary mophology that underlies wing shape in different avian clades, it is unlikely that a single skeletal feature will show consistent predictive power across Neoaves. Confident inference of wing shape in basal ornithurine birds will require multiple lines of evidence, together with an understanding of clade-specific evolutionary trends within the crown.

Background

Bird wings represent an extreme modification of the tetrapod forelimb to meet the functional demands of flight. The most apparent change involves the elongation and stiffening of feathers on the forelimb to form an airfoil. The full complement of specially-modified flight feathers contributes more than 85% of total wing area in most extant birds [1], and the overall wing shapes that result from the aggregate of individual flight feather shapes determine critical aspects of aerodynamic function (Figure 1d): Long, narrow wings (high aspect ratio, AR) typically lower the drag associated with creating lift, leading to more efficient sustained flight, while increasing wing area for a given body mass (lowering wing loading) leads to more efficient soaring flight at lower speeds [2-4].

Figure 1

Bony morphology of feather attachment and composition of wing shape. (a) – (c) Major bony landmarks as visible on a rock pigeon (Columba livia). (a) Dorsal view of the distal forelimb. I: digit I; II-1: proximal phalanx digit II; II-2: distal phalanx digit II; III: digit III; IIP: internal index process; cmc: carpometacarpus; r: radius; u: ulna. (b) 1° attachment features on digits II-III, dorsal view. dpf: dorsal phalangeal fossae (primaries VIII-IX); fp: flexor process of digit III (primaries I-VI); iip1: internal index process (primaries IX-X); iip2: position of caudal lobe of internal index process (primary IX); tvi: tendon of ventral interosseous muscle; tfdm: tendon of flexor digiti minimi muscle. (c) 2° attachments on the ulna, caudal view. drp: dorsal remigial papillae; crp: caudal remigial papillae. (d) Schematic of flight feather attachment and contribution to wing shape in a pink-footed shearwater Puffinus creatopus (left) and a roadrunner Geococcyx californianus (right), taxa with roughly equivalent wing area but very different wing shapes and functions. Shading indicates groups of feathers with distinct attachment points. Note that distal primaries form the pointed wingtip in Puffinus, while the rounded wingtip in Geococcyx is formed by more proximal primaries. pd1°: predigital primary (II-2); md1°: middigital primaries (II-1); ad1°: addigital primary (III); cm1°: carpometacarpal primaries (III + carpometacarpus + tfdm); 2°: secondaries. Bird images: Puffinus courtesy M. Taylor (CC-by-SA 3.0) 2009; Geococcyx courtesy G. Kramer and USFWS, 2009.

Extant birds stand out from the other living flyers (including bats and insects) in their ability to maintain attached, efficient airflow in flapping and soaring flight with high aspect ratio wing shapes and high wing loadings, even at comparatively high Reynolds numbers [1-3,5]. As with many other organismal traits, the evolution of components of flight performance in extant birds, including wing shape, is best understood as gradual change along a continuum or a stepwise accretion of traits [6]. Thus a nuanced understanding of avian wing shape evolution requires information from fossils to constrain evolutionary rates and ancestral character states in deep time [7,8] as a potential explanation for the current ecological and phylogenetic diversity of birds [9].

Efforts to reconstruct evolutionary pattern in wing shape are hindered by the fact that the entire wing is often not preserved as a body fossil, leaving only the pectoral skeletal elements as clues to shape and function in extinct taxa. Inferences of wing shape and flight mechanics in extinct taxa have been primarily based on relationships between whole wing shape parameters and the length of skeletal elements, such as the ratio of humeral and ulnar length [10], or the length of major skeletal segments [11-13]. Cross-sectional properties [14,15] and bone tissue histology [16,17] of both proximal and distal wing elements have also been tested as proxies for wing shape and flight performance.

In contrast, distal limb surface morphology and its relationship to wing shape have a limited history of investigation. Bird wing morphology is unique among the known flying vertebrates (which also include bats and pterosaurs) in that the distal extent of the wing is not demarcated by skeletal features [18,19]. Rather than a membrane that attaches along the length of the pectoral skeletal elements, bird wings are mainly composed of primary (1°) and secondary (2°) flight feathers (remiges, sing. remex), that have restricted, focal attachment points to the underlying skeleton (Figure 1a-c). These attachment points function as mechanical links that transfer aerodynamic forces generated by flight feathers to the limb skeleton [20-24]. Despite their role in the mechanics of flight, remex attachments have not been as thoroughly studied as other potential correlates of flight behavior. Skeletal variability in proximal attachment points (ulnar quill knobs, Figure 1c) has been linked to flightedness [25]. The distal attachments of the 1° feathers, specifically the length of the internal index process (iip1, Figure 1b) have been suggested as correlates for aspect ratio [26,27].

The distal portion of a flapping wing produces the greatest aerodynamic force per unit area [28]. The length and shape of distal 1° feathers is thus expected to have a disproportionate effect on force production and induced drag when compared to 2° feathers and proximal 1° feathers. My initial hypotheses are: (a) bony correlates of distal 1° attachment are more prominent in birds with higher AR, and (b) bony correlates of proximal 1° and 2° attachment are more prominent in birds with lower AR.

The prominence of bony attachment features in general may vary with body size. Mass (M), as a proxy for size, is a relevant covariate for any study of skeletal morphology. My initial expectation is that all bony features in the study will be scored as more prominent with increasing mass.

I used phylogenetic distance-based Redundancy Analysis (RDA) [29,30] and variation partitioning (VARPART) [31,32] as means to explore relationships among the several variables simultaneously. RDA works much like a principal components analysis (PCA), which is commonly used to examine several morphological variables at once. RDA has the added benefit of explicitly including information from expected functional relationships. Where PCA reports patterns of greatest morphological variability, RDA reports patterns of greatest morphological variability that are correlated with a second set of variables, in this case AR and M. VARPART then provides a means of measuring the magnitude of the effect that the second set of variables has on morphology.

After initial exploration of the data, I employed phylogenetic ANOVA [33] to test for consistent relationships between qualitative bony indicators of proximal and distal remex attachment on one hand (Figure 1a-c, Table 1), and aspect ratio and body mass on the other (Figure 1d, Table 1). Phylogenetic ANOVA provides a more conventional and widely used approach to dichotomous hypothesis testing. The results of these tests directly address the question of whether individual categorically-scored bony characters provide a record of AR that might be interpreted for extinct taxa.

Table 1 Character scores, AR , and M for taxa used in this study

Methods

Remex attachment morphology

Several skeletal features have been implicated in the attachment of flight feathers in past studies [20-24,26]. These features have often been employed in systematics, and as such are known from a broad taxonomic range of living birds. Published binary and multistate categorical scores coded at genus level [34] for six remigial attachment characters (Table 2) were chosen to represent bony features known to be associated with remex attachment. Published scores were then used as a basis to score several additional taxa for a pooled sample of n = 71 extant taxa (Additional file 1).

Table 2 Categorically scored bony characters included in this study

Published scores for similar characters were also used to score three extinct taxa: Parargornis messelensis [37], Eocypselus rowei [38], and Ichthyornis dispar [39]. Primary feathers are well-preserved for Parargornis and Eocypselus, which allows estimation of AR from the fossils. Ichthyornis represents an early occurrence of two feather attachment features (internal index process and dorsal phalangeal fossae, Table 2) that are homologous with extant examples.

Taxon scores for attachment morphology were submitted to principal coordinates analysis (PCO) using a Hamming distance metric with Cailliez correction for negative eigenvalues [29]. This analysis and all subsequent analyses were performed using R 2.15.2 [40]. PCO provides a semi-quantitative view of a set of categorical characters that is amenable to linear modeling [29]. Correlations between the original character scores and PCO scores provides a projection of the original categorical variables on the PCO axes [41] analogous to loadings in principal component analysis (PCA). The PCO scores for each taxon form the response matrix Y. Extinct taxa were included in this portion of the analysis.

Wing shape and body mass

Measurements of mass (M), wing span (b), and wing area (S) for the study taxa were taken from published sources [1,42,43]. Span and area were recombined into aspect ratio (AR, = b ²/S). AR and M were log_e transformed to achieve normality (tested with Shapiro-Wilk W, p _normal = 0.64 & 0.61, respectively). These two variables form the explanatory matrix X. Estimates of AR from Parargornis and Eocypselus were excluded from this matrix and from subsequent steps.

Phylogenetic context

This study accounted for phylogenetic error covariance using a phylogenetic generalized least squares (PGLS) approach [44]. A subset of 1,000 trees containing the study taxa was obtained from a published avian supertree set [9]. A separate phylogenetic covariance matrix was derived from the topology and branch lengths of each tree [45]. Phylogenetic signal (λ _i ) for each tree was jointly estimated for all variables [46]. Estimated values of λ _i were tested against λ _i  = 0 and λ _i  = 1 by a likelihood-ratio chi-squared test with one degree of freedom. Estimated values of λ _i were used to scale off-diagonal branch lengths in phylogenetic covariance matrices.

All variables were phylogenetically centered around an estimate of ancestral character state [47] for each tree i. Because AR and M have heterogeneous units, they were ranged to (−1, 1) after centering. Centered and ranged variables were then PGLS transformed using the inverse square root of each phylogenetic covariance matrix (C _i ^-1/2). Because the explanatory and response matrices X and Y are multiplied in the course of VARPART and RDA analyses, this transformation is equivalent to multiplication by C ⁻¹ in standard PGLS approaches. The transformed variables can be analyzed by standard parametric linear approaches while accounting for expected phylogenetic error covariance among related species.

Statistical analysis

Relationships between bone morphology and wing shape were initially evaluated using VARPART [31,32,48] and RDA [29,48]. In this study, VARPART partitions the variation in bony morphology into components that are ‘explained’ by correlation with AR and M, and a residual component that is not correlated to either. RDA then provides a PCA-like view of the components of bony morphological variation that are correlated to AR and M (the canonical axes) and the residual components (the non-canonical axes).

VARPART and RDA were run with an iteration to incorporate each of the i trees in the sample, to accommodate uncertainty in phylogeny. Using the set of 1,000 trees provides a distribution of adjusted R² scores and axis loadings. Comparison within this distribution serves as a sensitivity analysis of VARPART and RDA to assumptions of tree topology and branch length. As the ordination results are not nested, results for individual trees are not directly comparable. VARPART and RDA results for a single tree with adjusted R² values closest to the means for all terms was taken as a representative ordination.

Eigenvectors from RDA were used as a basis to plot taxon scores in constrained reduced space [15,30] on the canonical axes. In a sense, this step takes taxon scores from an unconstrained ordination that reflects patterns of shared variability in bony morphology (PCO), and projects those scores onto a new set of axes that reflect patterns of shared variability in bony morphology that are correlated with AR and M (RDA). The constrained reduced space provides a concise summary of variation in bone morphology that is directly tied to variation in wing shape. Extinct taxa were placed on this ordination space using their PCO scores.

One-to-one relationships between skeletal features and AR or M identified from the ordination were further tested in phylogenetic context using phylogenetic ANOVA [33,46].

Results

Representation of categorically scored bony morphology by Principal Coordinates

The first twelve PCO axes are used to represent bony morphology for subsequent analyses. Of these axes, the first eight each represent greater than 1% of total variability (Figure 2). Eigenvalues for PCO axes are presented in Table 3. Relationships between the unconstrained PCO space and the original binary or categorically scored bone morphology characters are presented in Table 4.

Figure 2

Biplots of bone morphology characters on successive PCO axes. Insets show bone feature loadings on Principal Coordinate Axes (a) 1-2, (b) 3-4, (c) 5-6, and (d) 7-8. drp: dorsal remigial papillae; crp: caudal remigial papillae; dpf: dorsal phalangeal fossae; iip1: presence/absence of internal index process; iip2: shape of internal index process; fp: flexor process digit III.

Table 3 Eigenvalues of PCO axes used to represent categorically scored bony characters

Table 4 Projection of the original morphological characters onto the first eight PCO axes

Phylogenetic signal

Joint estimates of λ _i for the study data are significantly different from both zero and one (p < 0.001) for all phylogenetic trees considered in the study, and are normally distributed around a mean λ _i  = 0.67 (Figure 3). Variable loadings in the redundancy analysis (RDA) ordination space for all trees are clustered around similar values (Figure 4).

Figure 3

Histogram showing distribution of joint λ _i values estimated for the tree sample.

Figure 4

Biplots showing the distribution of bone and wing shape character loadings for the tree sample. drp: dorsal remigial papillae; crp: caudal remigial papillae; dpf: dorsal phalangeal fossae; iip1: presence/absence of internal index process; iip2: shape of internal index process; fp: flexor process digit III; AR: aspect ratio; M: mass.

Aspect ratio is linked to variability in skeletal features of the manus

VARPART and RDA identified significant components of variation in bony feather attachment morphology that are best explained by AR (R² = 0.065, Table 5). Although the coefficient of determination is low, this reflects the relationship of AR to all of the included skeletal features simultaneously. Permutation tests identified both the first canonical axis (Figure 5a,b) and the loadings of AR and M along the first axis as significantly different from random (p ≤ 0.05, Table 6). These results, along with the loadings shown in Figure 5b, show that although the statistical model has poor predictive power, there is an association between AR and more prominent bony feather attachment features in the hand. This association is convergent across multiple neoavian lineages (Figure 5a,d).

Table 5 VARPART results for single representative tree, with the range observed in the tree sample

Figure 5

Summary of VARPART and RDA results. (a) Plot of taxon RDA scores. Hue and shading of points denotes aspect ratio. Silhouettes show representative wing shapes for Anser, Cathartes, Sula, and Fregata. (b) Inset biplot representing PGLS correlations among response matrix Y (solid rays) and explanatory matrix X (dashed rays) in the canonical RDA axes. Note the strong positive correlation between aspect ratio (AR) and distal primary attachment features (iip2). drp: dorsal remigial papillae; crp: caudal remigial papillae; dpf: dorsal phalangeal fossae; iip1: presence/absence of internal index process; iip2: shape of internal index process; fp: flexor process digit III. (c) – (d) Proximal phalanges of digit II (II-1) from Rallus and Fregata in dorsal view, show morphological disparity at the extremes of the first RDA axis. (e) Representative tree topology with an ancestral character state reconstruction of AR—hues and shading correspond to (a).

Table 6 RDA axis eigenvalues for single representative tree

Although M has a non-random effect on the ordination of bony morphology (Table 5), mass primarily loads along the second canonical axis (Figure 5b), which is only marginally distinct from a random rearrangement of the original PCO axes (p = 0.07).

Shape of the caudal margin of digit II is a correlate of aspect ratio and mass

The presence of a caudal lobe of the internal index process (iip2, Figure 1b) is associated with significantly different AR and M (p = 0.001 & p = 0.013, respectively; phylogenetic ANOVA using the representative tree with 1,000 simulations and Holm multiple comparisons adjustment). As scored for this analysis, presence or absence of the flexor process of digit III (fp) and the prominence of dorsal phalangeal fossae (dpf) do not show significant relationships to AR (p = 0.8 & p = 0.4, respectively) or M (p = 0.07 & p = 0.8).

Modeled AR for extinct taxa are inconclusive

Fitted AR for Eocypselus, Parargornis (7.3), and Ichythyornis (7.8) were all close to the value inferred for the root of the neoavian tree in the representative ordination (7.6, Figure 5d). These values represent an underestimation for Eocypselus (inferred AR of 8.2 from feather impressions) [38], and an overestimation for Parargornis (inferred AR of 5.3) [37]. The weak fit of AR to the full RDA model, which contains all of the morphological features, does not allow any conclusions to be drawn regarding wing shape in Ichthyornis.

Discussion

Variation in feather attachment morphology

Loadings of the PCO axes identify prominent patterns of variation in bony flight feather attachment features in living birds (Figure 2a-d). The first PCO axis captures a general trend of linked presence or absence among most of the bony characters. The second PCO axis shows a pattern of divergence in which several taxa possess caudal remigial papillae (crp) but not dorsal phalangeal fossae (dpf), and vice versa, which points to a trade-off between the prominence of 2° and distal 1° feather attachments, respectively. The third PCO axis captures a similar tradeoff between distal (iip2, dpf) and proximal (fp) 1° feather attachments. The fourth PCO axis accounts for variability in dorsal remigial papillae (drp) that does not co-occur with any of the above patterns. The remaining axes account for residual variation that occurs as exceptions to the abovementioned trends.

Phalangeal morphology and distal primary attachment

The simple presence/absence character score for iip1 (Figure 1b) used in this study obscures a broad taxonomic range of shape and relative size in this feature. The ‘caudal lobe’ (iip2) [34] is only found in taxa with very prominent internal index processes (IIP), and thus may also serve as a crude proxy for the length of the IIP, as defined in the strict sense. Because the IIP and its caudal lobe are attachment sites for ligaments that hold primaries VIII, IX and X to the major digit, the significant relationship between these features and AR suggests a relatively straightforward interpretation: longer distal 1° feathers contribute to increased span, drive up AR (Figure 1d), and presumably exert greater stress on their attachment sites, resulting in more prominent bony features (Figure 5d).

The relationship between IIP size and AR has been suggested to follow from a hypothesized role of the IIP as a guide for the tendon of the ventral interosseus muscle (tvi, Figure 1b) [26]. By this interpretation, increased IIP length would result in increased moment arm for adduction of the distal phalanx and its attached predigital 1° feather. However, the potential performance benefits of increased moment arm for predigital 1° adduction in flapping flight are unclear, and the role of the IIP as a remigial ligament attachment site, although supported by anatomical data, has not previously been considered. Further anatomical and histological studies will undoubtedly shed more light on the functional significance of the IIP.

Ulnar remigial papillae

The prominence of ulnar remigial papillae does not show a clear relationship to AR or M (Figure 5b). Biplots of the original morphological variables on the non-canonical (residual) axes from RDA (Figure 6) show that crp, dpf, and drp show independent patterns of variation on the first, second, and third non-canonical axes, respectively. Distal limb pneumatization (‘hyperpneumatization’ [49]) is a possible confounding variable for these features, as many of the hyperpneumatized taxa in the study also show prominent drp and dpf. In many birds with distal limb pneumaticity, subcutaneous airsacs lie adjacent to the follicles of 1° and 2° remiges at their attachments [49]. Cross-sections of drp in these taxa reveal substantial erosion rooms within the papillae themselves that are independent of the medullary cavity of the ulna (Figure 7). In a similar manner, dpf of hyperpneumatized taxa (Figure 5d) show prominent excavations associated with air sacs surrounding follicles of the distal 1° remiges.

Figure 6

Non-canonical axes from RDA, representing residual variation in bony morphology that is not related to AR or M . Plots show non-canonical axes (a) 1-2, and (b) 3-4.

Figure 7

MicroCT surfaces and slices through dorsal remigial papillae ( drp ). (a-b) Fregata (CM S13606) and (c-d) Cathartes, showing open spaces isolated from the main medullary cavity of the ulna. Cathartes specimen is embedded in plastic for histological sectioning—feather pulp cavities appear as void spaces in the slice. crp, caudal remigial papillae; 2°, pulp cavity of secondary remex.

A single interpretation of quill knobs cannot be extended to address all taxa included in the present study. In some taxa (e.g. Ardea), prominent drp are not associated with pneumatization of the distal forelimb—the same is also true of the dpf of Rissa, Stercorarius, and Sterna [49]. Although the morphological pattern of quill knobs scored in this study suggests that the overlying system of feather attachment ligaments is homologous among extant birds, the patterns of bone growth and resorption that produce different bony character states are clearly homoplastic. Quill knobs may appear similar in these independent occurrences, but their presence may very well serve different purposes in different clades.

Wing shape and many-to-one mapping

The coefficient of determination for AR found by VARPART (Table 5) is noteworthy despite its low value for two reasons. First, the preservation of any signal that can be detected through the noise and low resolution of several anatomically independent, categorically scored characters across the diversity of neoavian birds is encouraging. The relationships between phalangeal morphology and AR detected by RDA warrant further investigation with more precise measures of morphology.

Second, although the relationship between discrete bony attachment characters and AR is weak, it is not a great deal weaker than previous models that employed brachial index (BI, humerus length/ulna length) as a predictor of AR when phylogeny was taken into account [50]. Difficulty in predicting wing shape from skeletal features in general may be a reflection of distinct functional relationships within different neoavian clades. If this is the case, predictive power may only be found by considering more phylogenetically restricted samples.

Extant birds may also display a pattern of many-to-one mapping [51], in which several distinct combinations of element and feather lengths may result in functionally equivalent wings with similar AR. Wing span, as a component of AR, is largely a summation of the lengths of humerus, ulna, manus, and distal 1° feathers (Figure 1d). The low predictive power of both the current RDA model and BI may reflect the fact that each model excludes important skeletal information from either the proximal or distal segments of the wing. The consilience of results from this study and previous work points in two directions for future studies: look for relationships in more restricted regions of the avian tree, and begin with skeletal data that address all four of the variable-length segments in the avian wing.

Conclusions

Although the ulnar papillae have been the most frequently cited correlates of feather attachment, their link to overall wing shape is problematic, and likely driven by divergent trends in different clades of neoavian birds. Homoplasy of prominent ulnar papillae may be better thought of as a reflection of conserved ligament architecture, rather than a convergent flight-related functional signal.

The internal index process of digit II-1 and associated attachment structures have previously been overlooked as correlates of flight feather attachment, but show a clear relationship to aspect ratio and the shape of the distal primaries. This relationship is apparent even at the low level of precision available in binary and categorical character scores. A more fine-grained approach to quantifying the size and shape of distal primary attachment features would provide a further test of the sensitivity of phalangeal morphology as a correlate of wing shape.

Feather attachment features provide a novel line of evidence for wing shape and function, alongside previously established relationships such as element length indices [10-13], cross-sectional properties [14,15], and bone tissue organization [16]. Combining evidence from feather attachment morphology with other forms of skeletal evidence, and restricting the phylogenetic scope of modeling efforts, may yet yield accurate skeletal predictors of wing shape to shed light on the evolution of flight mechanics within crown neoavian birds as a dynamic, divergent process [52].

Availability of supporting data

The data sets and R language scripts supporting the results of this article are available in Dryad (http://dx.doi.org/10.5061/dryad.2p31f).