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  <cit> , and the overall wing shapes that result from the aggregate of individual flight feather shapes determine critical aspects of aerodynamic function : long, narrow wings  typically lower the drag associated with creating lift, leading to more efficient sustained flight, while increasing wing area for a given body mass  leads to more efficient soaring flight at lower speeds  <cit> .figure  <dig>
bony morphology of feather attachment and composition of wing shape.
–  major bony landmarks as visible on a rock pigeon .
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.
1° attachment features on digits ii-iii, dorsal view.
dpf: dorsal phalangeal fossae ; fp: flexor process of digit iii ; iip1: internal index process ; iip2: position of caudal lobe of internal index process ; tvi: tendon of ventral interosseous muscle; tfdm: tendon of flexor digiti minimi muscle.
2° attachments on the ulna, caudal view.
drp: dorsal remigial papillae; crp: caudal remigial papillae.
schematic of flight feather attachment and contribution to wing shape in a pink-footed shearwater puffinus creatopus  and a roadrunner geococcyx californianus , 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 ; md1°: middigital primaries ; ad1°: addigital primary ; cm1°: carpometacarpal primaries ; 2°: secondaries.
bird images: puffinus courtesy m. taylor  2009; geococcyx courtesy g. kramer and usfws,  <dig>
extant birds stand out from the other living flyers  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  <cit> .
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  <cit> .
thus a nuanced understanding of avian wing shape evolution requires information from fossils to constrain evolutionary rates and ancestral character states in deep time  <cit>  as a potential explanation for the current ecological and phylogenetic diversity of birds  <cit> .
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  <cit> , or the length of major skeletal segments  <cit> .
cross-sectional properties  <cit>  and bone tissue histology  <cit>  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  in that the distal extent of the wing is not demarcated by skeletal features  <cit> .
rather than a membrane that attaches along the length of the pectoral skeletal elements, bird wings are mainly composed of primary  and secondary  flight feathers , that have restricted, focal attachment points to the underlying skeleton .
these attachment points function as mechanical links that transfer aerodynamic forces generated by flight feathers to the limb skeleton  <cit> .
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  has been linked to flightedness  <cit> .
the distal attachments of the 1° feathers, specifically the length of the internal index process  have been suggested as correlates for aspect ratio  <cit> .
the distal portion of a flapping wing produces the greatest aerodynamic force per unit area  <cit> .
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:  bony correlates of distal 1° attachment are more prominent in birds with higher ar, and  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 , 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   <cit>  and variation partitioning   <cit>  as means to explore relationships among the several variables simultaneously.
rda works much like a principal components analysis , 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  <cit>  to test for consistent relationships between qualitative bony indicators of proximal and distal remex attachment on one hand , and aspect ratio and body mass on the other .
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  <dig>
character scores,
ar
, and
m
for taxa used in this study
taxon
drp
crp
dpf
iip1
iip2
fp
ar
m
accipiter cooperii
accipiter gentilis
accipiter nisus
accipiter striatus
anas crecca
anas georgica
anas penelope
anas platyrhynchos
anhinga anhinga
anser albifrons
anser anser
anser indicus
ardea cinerea
ardea herodias
caracara plancus
cathartes aura
chionis albus
ciconia abdimii
coccyzus americanus
columba palumbus
diomedea exulans
falco columbarius
falco mexicanus
falco peregrinus
falco rusticolus
falco sparverius
fregata magnificens
fulmarus glacialis
gavia immer
gavia stellata
gyps africanus
gyps rueppellii
haematopus ostralegus
limnodromus griseus
megascops asio
mergus serrator
oceanites oceanicus
pachyptila desolata
pachyptila turtur
passer domesticus
pelecanoides georgicus
pelecanoides urinatrix
pelecanus erythrorhynchos
pelecanus occidentalis
pelecanus onocrotalus
pelecanus rufescens
phaethon aethereus
phaethon lepturus
phaethon rubricauda
phalacrocorax aristotelis
phalacrocorax atriceps
phalacrocorax auritus
phalacrocorax carbo
phasianus colchicus
platalea ajaja
plegadis chihi
puffinus huttoni
puffinus lherminieri
puffinus nativitatis
puffinus pacificus
puffinus tenuirostris
rallus longirostris
rissa tridactyla
stercorarius parasiticus
sterna maxima
sula dactylatra
sula leucogaster
sula sula
tyto alba
uria aalge
zenaida macroura
ichthyornis dispar
eocypselus rowei
parargornis messelensis
methods
remex attachment morphology
several skeletal features have been implicated in the attachment of flight feathers in past studies  <cit> .
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  <cit>  for six remigial attachment characters  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 =  <dig> extant taxa .table  <dig>
categorically scored bony characters included in this study
character
lz <dig> #
soft tissue attachment
related remiges
dorsal remigial papillae a-e
caudal remigial papillae a-e
dorsal phalangeal fossae b,d,f
internal index process presence/absence f
internal index process shape f
flexor process g
characters related to feather attachment have been identified by several sources .
corresponding characters described in  <cit>   included as an aggregate reference.
a <cit> ; b <cit> ; c <cit> ; d <cit> ; e <cit> ; f <cit> ; g <cit> ; h <cit> .
published scores for similar characters were also used to score three extinct taxa: parargornis messelensis  <cit> , eocypselus rowei  <cit> , and ichthyornis dispar  <cit> .
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  that are homologous with extant examples.
taxon scores for attachment morphology were submitted to principal coordinates analysis  using a hamming distance metric with cailliez correction for negative eigenvalues  <cit> .
this analysis and all subsequent analyses were performed using r  <dig> .
<dig>  <cit> .
pco provides a semi-quantitative view of a set of categorical characters that is amenable to linear modeling  <cit> .
correlations between the original character scores and pco scores provides a projection of the original categorical variables on the pco axes  <cit>  analogous to loadings in principal component analysis .
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 , wing span , and wing area  for the study taxa were taken from published sources  <cit> .
span and area were recombined into aspect ratio .
ar and m were loge transformed to achieve normality .
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  approach  <cit> .
a subset of  <dig>  trees containing the study taxa was obtained from a published avian supertree set  <cit> .
a separate phylogenetic covariance matrix was derived from the topology and branch lengths of each tree  <cit> .
phylogenetic signal  for each tree was jointly estimated for all variables  <cit> .
estimated values of λi were tested against λi =  <dig> and λi =  <dig> 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  <cit>  for each tree i. because ar and m have heterogeneous units, they were ranged to  after centering.
centered and ranged variables were then pgls transformed using the inverse square root of each phylogenetic covariance matrix .
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− <dig> 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  <cit>  and rda  <cit> .
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  and the residual components .
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  <dig>  trees provides a distribution of adjusted r <dig> 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 <dig> 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  <cit>  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 , 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 .
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  <cit> .
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 .
eigenvalues for pco axes are presented in table  <dig>  relationships between the unconstrained pco space and the original binary or categorically scored bone morphology characters are presented in table  <dig> figure  <dig>
biplots of bone morphology characters on successive pco axes.
insets show bone feature loadings on principal coordinate axes  1- <dig>   3- <dig>   5- <dig>  and  7- <dig>  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.
eigenvalues of pco axes used to represent categorically scored bony characters
relative eigenvalue
cumulative eigenvalue
pco1
pco2
pco3
pco4
pco5
pco6
pco7
pco8
drp
crp
dpf
iip1
iip2
fp
this matrix was used to transform the loadings of pco axes in the rda ordination into a biplot of the original morphological variables .
phylogenetic signal
joint estimates of λi for the study data are significantly different from both zero and one  for all phylogenetic trees considered in the study, and are normally distributed around a mean λi =  <dig>  .
variable loadings in the redundancy analysis  ordination space for all trees are clustered around similar values .figure  <dig>
histogram showing distribution of joint
λ
i
values estimated for the tree sample.
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 .
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  and the loadings of ar and m along the first axis as significantly different from random .
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 .table  <dig>
varpart results for single representative tree, with the range observed in the tree sample
model
adjusted r
2
upper adj.
r
2
lower adj.
r
2
ar
m
ar + m
ar | m
m | ar
ar ∩ m
**significant by permutation test at p ≤  <dig> .
*significant by permutation test at p ≤  <dig> .
summary of varpart and rda results.
plot of taxon rda scores.
hue and shading of points denotes aspect ratio.
silhouettes show representative wing shapes for anser, cathartes, sula, and fregata.
inset biplot representing pgls correlations among response matrix y  and explanatory matrix x  in the canonical rda axes.
note the strong positive correlation between aspect ratio  and distal primary attachment features .
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.
–  proximal phalanges of digit ii  from rallus and fregata in dorsal view, show morphological disparity at the extremes of the first rda axis.
representative tree topology with an ancestral character state reconstruction of ar—hues and shading correspond to .
rda axis eigenvalues for single representative tree
axis
eigenvalue
relative eigenvalue
**significant by permutation test at p ≤  <dig> .
although m has a non-random effect on the ordination of bony morphology , mass primarily loads along the second canonical axis , which is only marginally distinct from a random rearrangement of the original pco axes .
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  is associated with significantly different ar and m .
as scored for this analysis, presence or absence of the flexor process of digit iii  and the prominence of dorsal phalangeal fossae  do not show significant relationships to ar  or m .
modeled ar for extinct taxa are inconclusive
fitted ar for eocypselus, parargornis , and ichythyornis  were all close to the value inferred for the root of the neoavian tree in the representative ordination .
these values represent an underestimation for eocypselus   <cit> , and an overestimation for parargornis   <cit> .
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 .
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  but not dorsal phalangeal fossae , 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  and proximal  1° feather attachments.
the fourth pco axis accounts for variability in dorsal remigial papillae  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 iip <dig>  used in this study obscures a broad taxonomic range of shape and relative size in this feature.
the ‘caudal lobe’   <cit>  is only found in taxa with very prominent internal index processes , 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 , and presumably exert greater stress on their attachment sites, resulting in more prominent bony features .
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   <cit> .
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 .
biplots of the original morphological variables on the non-canonical  axes from rda  show that crp, dpf, and drp show independent patterns of variation on the first, second, and third non-canonical axes, respectively.
distal limb pneumatization  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  <cit> .
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 .
in a similar manner, dpf of hyperpneumatized taxa  show prominent excavations associated with air sacs surrounding follicles of the distal 1° remiges.figure  <dig>
non-canonical axes from rda, representing residual variation in bony morphology that is not related to
ar
or
m
.
plots show non-canonical axes  1- <dig>  and  3- <dig>
microct surfaces and slices through dorsal remigial papillae .
fregata  and
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 , prominent drp are not associated with pneumatization of the distal forelimb—the same is also true of the dpf of rissa, stercorarius, and sterna  <cit> .
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  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  as a predictor of ar when phylogeny was taken into account  <cit> .
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  <cit> , 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 .
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- <dig> 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  <cit> , cross-sectional properties  <cit> , and bone tissue organization  <cit> .
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  <cit> .
availability of supporting data
the data sets and r language scripts supporting the results of this article are available in dryad .
additional file
additional file 1:
hieronymus flight feather table s <dig> xls
: excel spreadsheet of taxa included in the study.
coding abbreviations: lze—livezey & zusi  exemplar taxon; lzc—livezey & zusi  exemplar congener; ost—coded from osteological specimen.
institutional abbreviations: cm—carnegie museum of natural history.
wing data reference abbreviations: w’77—warheit ; hb’99—hertel & ballance ; p’01—pennycuick ; p’08—pennycuick .
abbreviations
1primary feather
2secondary feather
ad1addigital primary feather
cm1carpometacarpal primary feathers
cmccarpometacarpus
crpcaudal remigial papillae
dpfdorsal phalangeal fossae
drpdorsal remigial papillae
fpflexor process digit iii
imanus digit i
ii-1proximal phalanx of manus digit ii
ii-2distal phalanx of manus digit ii
iiimanus digit iii
iipinternal index process
iip1internal index process
iip2caudal lobe of internal index process
mbody mass
md1middigital primary feathers
pcoprincipal coordinates analysis
pd1predigital primary feather
pglsphylogenetic generalized least-squares
rradius
araspect ratio
rdaredundancy analysis
tfdmtendon of flexor digiti minimi muscle
tvitendon of ventral interosseous muscle
uulna
varpartvariation partitioning
competing interests
the author declares that he has no competing interests.

 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, and the overall wing shapes that result from the aggregate of individual flight feather shapes determine critical aspects of aerodynamic function: long, narrow wings typically lower the drag associated with creating lift, leading to more efficient sustained flight, while increasing wing area for a given body mass leads to more efficient soaring flight at lower speeds.figure bony morphology of feather attachment and composition of wing shape. – major bony landmarks as visible on a rock pigeon. 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. 1° attachment features on digits ii-iii, dorsal view. dpf: dorsal phalangeal fossae; fp: flexor process of digit iii; iip1: internal index process; iip2: position of caudal lobe of internal index process; tvi: tendon of ventral interosseous muscle; tfdm: tendon of flexor digiti minimi muscle. 2° attachments on the ulna, caudal view. drp: dorsal remigial papillae; crp: caudal remigial papillae. schematic of flight feather attachment and contribution to wing shape in a pink-footed shearwater puffinus creatopus and a roadrunner geococcyx californianus, 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; md1°: middigital primaries; ad1°: addigital primary; cm1°: carpometacarpal primaries; 2°: secondaries. bird images: puffinus courtesy m. taylor 2009; geococcyx courtesy g. kramer and usfws, extant birds stand out from the other living flyers 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. 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. thus a nuanced understanding of avian wing shape evolution requires information from fossils to constrain evolutionary rates and ancestral character states in deep time as a potential explanation for the current ecological and phylogenetic diversity of birds. 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, or the length of major skeletal segments. cross-sectional properties and bone tissue histology 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 in that the distal extent of the wing is not demarcated by skeletal features. rather than a membrane that attaches along the length of the pectoral skeletal elements, bird wings are mainly composed of primary and secondary flight feathers, that have restricted, focal attachment points to the underlying skeleton. these attachment points function as mechanical links that transfer aerodynamic forces generated by flight feathers to the limb skeleton. 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 has been linked to flightedness. the distal attachments of the 1° feathers, specifically the length of the internal index process have been suggested as correlates for aspect ratio. the distal portion of a flapping wing produces the greatest aerodynamic force per unit area. 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: bony correlates of distal 1° attachment are more prominent in birds with higher ar, and 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, 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 and variation partitioning as means to explore relationships among the several variables simultaneously. rda works much like a principal components analysis, 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 to test for consistent relationships between qualitative bony indicators of proximal and distal remex attachment on one hand, and aspect ratio and body mass on the other. 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 character scores, ar, and m for taxa used in this study taxon drp crp dpf iip1 iip2 fp ar m accipiter cooperii accipiter gentilis accipiter nisus accipiter striatus anas crecca anas georgica anas penelope anas platyrhynchos anhinga anhinga anser albifrons anser anser anser indicus ardea cinerea ardea herodias caracara plancus cathartes aura chionis albus ciconia abdimii coccyzus americanus columba palumbus diomedea exulans falco columbarius falco mexicanus falco peregrinus falco rusticolus falco sparverius fregata magnificens fulmarus glacialis gavia immer gavia stellata gyps africanus gyps rueppellii haematopus ostralegus limnodromus griseus megascops asio mergus serrator oceanites oceanicus pachyptila desolata pachyptila turtur passer domesticus pelecanoides georgicus pelecanoides urinatrix pelecanus erythrorhynchos pelecanus occidentalis pelecanus onocrotalus pelecanus rufescens phaethon aethereus phaethon lepturus phaethon rubricauda phalacrocorax aristotelis phalacrocorax atriceps phalacrocorax auritus phalacrocorax carbo phasianus colchicus platalea ajaja plegadis chihi puffinus huttoni puffinus lherminieri puffinus nativitatis puffinus pacificus puffinus tenuirostris rallus longirostris rissa tridactyla stercorarius parasiticus sterna maxima sula dactylatra sula leucogaster sula sula tyto alba uria aalge zenaida macroura ichthyornis dispar eocypselus rowei parargornis messelensis methods remex attachment morphology several skeletal features have been implicated in the attachment of flight feathers in past studies. 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 for six remigial attachment characters 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 =  extant taxa.table categorically scored bony characters included in this study character lz # soft tissue attachment related remiges dorsal remigial papillae a-e caudal remigial papillae a-e dorsal phalangeal fossae b,d,f internal index process presence/absence f internal index process shape f flexor process g characters related to feather attachment have been identified by several sources. corresponding characters described in included as an aggregate reference. a; b; c; d; e; f; g; h. published scores for similar characters were also used to score three extinct taxa: parargornis messelensis, eocypselus rowei, and ichthyornis dispar. 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 that are homologous with extant examples. taxon scores for attachment morphology were submitted to principal coordinates analysis using a hamming distance metric with cailliez correction for negative eigenvalues. this analysis and all subsequent analyses were performed using r.. pco provides a semi-quantitative view of a set of categorical characters that is amenable to linear modeling. correlations between the original character scores and pco scores provides a projection of the original categorical variables on the pco axes analogous to loadings in principal component analysis. 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, wing span, and wing area for the study taxa were taken from published sources. span and area were recombined into aspect ratio. ar and m were loge transformed to achieve normality. 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 approach. a subset of trees containing the study taxa was obtained from a published avian supertree set. a separate phylogenetic covariance matrix was derived from the topology and branch lengths of each tree. phylogenetic signal for each tree was jointly estimated for all variables. estimated values of λi were tested against λi =  and λi =  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 for each tree i. because ar and m have heterogeneous units, they were ranged to after centering. centered and ranged variables were then pgls transformed using the inverse square root of each phylogenetic covariance matrix. 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 and rda. 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 and the residual components. 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 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 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, 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. 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.  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. eigenvalues for pco axes are presented in table  relationships between the unconstrained pco space and the original binary or categorically scored bone morphology characters are presented in table  figure biplots of bone morphology characters on successive pco axes. insets show bone feature loadings on principal coordinate axes 1- 3- 5- and 7- 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. eigenvalues of pco axes used to represent categorically scored bony characters relative eigenvalue cumulative eigenvalue pco1 pco2 pco3 pco4 pco5 pco6 pco7 pco8 drp crp dpf iip1 iip2 fp this matrix was used to transform the loadings of pco axes in the rda ordination into a biplot of the original morphological variables. phylogenetic signal joint estimates of λi for the study data are significantly different from both zero and one for all phylogenetic trees considered in the study, and are normally distributed around a mean λi = . variable loadings in the redundancy analysis ordination space for all trees are clustered around similar values.figure histogram showing distribution of joint λ i values estimated for the tree sample. 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. 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 and the loadings of ar and m along the first axis as significantly different from random. 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.table varpart results for single representative tree, with the range observed in the tree sample model adjusted r 2 upper adj. r 2 lower adj. r 2 ar m ar + m ar | m m | ar ar ∩ m **significant by permutation test at p ≤ . *significant by permutation test at p ≤ . summary of varpart and rda results. plot of taxon rda scores. hue and shading of points denotes aspect ratio. silhouettes show representative wing shapes for anser, cathartes, sula, and fregata. inset biplot representing pgls correlations among response matrix y and explanatory matrix x in the canonical rda axes. note the strong positive correlation between aspect ratio and distal primary attachment features. 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. – proximal phalanges of digit ii from rallus and fregata in dorsal view, show morphological disparity at the extremes of the first rda axis. representative tree topology with an ancestral character state reconstruction of ar—hues and shading correspond to. rda axis eigenvalues for single representative tree axis eigenvalue relative eigenvalue **significant by permutation test at p ≤ . although m has a non-random effect on the ordination of bony morphology, mass primarily loads along the second canonical axis, which is only marginally distinct from a random rearrangement of the original pco axes. 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 is associated with significantly different ar and m. as scored for this analysis, presence or absence of the flexor process of digit iii and the prominence of dorsal phalangeal fossae do not show significant relationships to ar or m. modeled ar for extinct taxa are inconclusive fitted ar for eocypselus, parargornis, and ichythyornis were all close to the value inferred for the root of the neoavian tree in the representative ordination. these values represent an underestimation for eocypselus, and an overestimation for parargornis. 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. 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 but not dorsal phalangeal fossae, 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 and proximal 1° feather attachments. the fourth pco axis accounts for variability in dorsal remigial papillae 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 iip used in this study obscures a broad taxonomic range of shape and relative size in this feature. the ‘caudal lobe’ is only found in taxa with very prominent internal index processes, 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, and presumably exert greater stress on their attachment sites, resulting in more prominent bony features. 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. 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. biplots of the original morphological variables on the non-canonical axes from rda show that crp, dpf, and drp show independent patterns of variation on the first, second, and third non-canonical axes, respectively. distal limb pneumatization 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. 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. in a similar manner, dpf of hyperpneumatized taxa show prominent excavations associated with air sacs surrounding follicles of the distal 1° remiges.figure non-canonical axes from rda, representing residual variation in bony morphology that is not related to ar or m. plots show non-canonical axes 1- and 3- microct surfaces and slices through dorsal remigial papillae. fregata and 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, prominent drp are not associated with pneumatization of the distal forelimb—the same is also true of the dpf of rissa, stercorarius, and sterna. 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 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 as a predictor of ar when phylogeny was taken into account. 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, 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. 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.  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- 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, cross-sectional properties, and bone tissue organization. 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. availability of supporting data the data sets and r language scripts supporting the results of this article are available in dryad. additional file additional file 1: hieronymus flight feather table s xls: excel spreadsheet of taxa included in the study. coding abbreviations: lze—livezey & zusi exemplar taxon; lzc—livezey & zusi exemplar congener; ost—coded from osteological specimen. institutional abbreviations: cm—carnegie museum of natural history. wing data reference abbreviations: w’77—warheit; hb’99—hertel & ballance; p’01—pennycuick; p’08—pennycuick. abbreviations 1primary feather 2secondary feather ad1addigital primary feather cm1carpometacarpal primary feathers cmccarpometacarpus crpcaudal remigial papillae dpfdorsal phalangeal fossae drpdorsal remigial papillae fpflexor process digit iii imanus digit i ii-1proximal phalanx of manus digit ii ii-2distal phalanx of manus digit ii iiimanus digit iii iipinternal index process iip1internal index process iip2caudal lobe of internal index process mbody mass md1middigital primary feathers pcoprincipal coordinates analysis pd1predigital primary feather pglsphylogenetic generalized least-squares rradius araspect ratio rdaredundancy analysis tfdmtendon of flexor digiti minimi muscle tvitendon of ventral interosseous muscle uulna varpartvariation partitioning competing interests the author declares that he has no competing interests. 
an accessory lobe of the internal index process of digit ii- a bony correlate of distal primary attachment, shows weak but statistically significant relationships to aspect ratio and mass independent of other skeletal morphology. however, because the feathers that define the wingtip 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. 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. quill knobs on the ulna are examined concurrently, but do not show consistent signal with respect to wing shape. 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 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. increasing aspect ratio results in a functional tradeoff between low induced drag and increased wing inertia. avesfeatherosteologyforelimb evolutionphylogenetic comparative methodsquill knobsissue-copyright-statement© the author 2015 among living fliers, birds display relatively high aspect ratios, a dimensionless shape variable that distinguishes long and narrow vs. short and broad wings. the dorsal phalangeal fossae of digit ii- which house distal primaries viii and ix, also show a trend of increased prominence with higher aspect ratio. in contrast, correlates of distal primary feather attachment on the major digit show convergent responses to increasing aspect ratio.
 among living fliers, 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 and increased wing inertia. 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 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.  an accessory lobe of the internal index process of digit ii- 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- 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.  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. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords avesfeatherosteologyforelimb evolutionphylogenetic comparative methodsquill knobsissue-copyright-statement© the author 2015 