Target Taxa and Measurement Data

Measurement data for modern anatids were taken from skeletal specimens of 787 individuals representing 103 modern species or subspecies (Appendix Table 10 and  Supplemental Material Table S13). Taxonomy follows del Hoyo and Collar (2014). Subspecies were distinguished for the Branta canadensis–hutchinsii complex to accommodate great intraspecific size variation but were pooled for other species. All these species and subspecies are referred to as “species” for convenience. They were assigned to 2 groups defined a priori: the flightless group, which consists of 4 species, T. leucocephalus, T. pteneres, T. brachypterus, and A. aucklandica; and the volant group, which includes all other species. Measurements of limb bones of another flightless species, A. nesiotis, were taken from Williams (2015b) and were used for a posteriori evaluations of classification rules. Some specimens of Tachyeres spp. were identified from osteological keys given by Livezey and Humphrey (1992). Fulton et al. (2012) considered that the volant Falkland Islands population formerly identified as T. patachonicus was probably conspecific with the parapatric flightless population of T. brachypterus, which together would form a partially flightless species. For the purposes of the present study, individuals of T. patachonicus from the Falkland Islands were excluded from analysis and T. brachypterus was treated as a totally flightless species. Only skeletally mature specimens, as assessed by surface textures of long bones (Tumarkin-Deratzian et al. 2006, Watanabe and Matsuoka 2013), were included in the analysis. Specimens of captive individuals were measured, but all bones showing signs of pathology were excluded from the analysis. When available, the means of dimensions of left and right sides were used as individual values. For each species and for each variable, a species mean value was calculated from all available individual values, including those from individuals lacking some dimensions (i.e. incomplete skeletons), pooling sexes and localities. The modern sample covers all extant subfamilies and tribes recognized by Kear (2005) and spans almost the entire range of body size of modern anatids, from the smallest Nettapus auritus (mean body mass 266 g) to the largest Cygnus buccinator (mean body mass 11,110 g; body mass data from Dunning 2008).

Skeletal measurements were also taken on 365 specimens of 11 species of fossil anatids (Appendix Table 11 and  Supplemental Material Table S13). Measurement data for 5 additional species were taken from the literature and included in the analyses (Worthy and Holdaway 2002, Matsuoka et al. 2004, Mayr and Pavia 2014, Williams et al. 2014, Williams 2015b), although it should be noted that these measurements may have been taken in different ways. Several putatively flightless fossil anatids, including Garganornis ballmanni, Megalodytes morejohni, Cygnus falconeri, Chelychelynechen quassus, Thambetochen xanion, and Talpanas lippa, were not included in this study because of the unavailability of measurements of most major limb bones (at least from published materials).

Measurements were taken with a Mitutoyo digital caliper to 0.01 mm and rounded to ±0.1 mm for bones up to ∼200 mm, and with a metal ruler with a finite end and rounded to ±1 mm for those >200 mm. The following 7 dimensions were measured (Figure 1; anatomical terminology follows Baumel and Witmer 1993): carinal height, the dorsoventral depth of carina sterni measured along the cranial margin (pila carinae); humerus length, the greatest length from caput humeri to condylus ventralis; ulna length, the greatest length from olecranon to condylus dorsalis ulnae; carpometacarpus length, the greatest length from trochlea carpalis to facies articularis digitalis major; femur length, the greatest length from trochanter femoris to condylus lateralis; tibiotarsus length, length from facies articularis lateralis (not from crista cnemialis cranialis) to condylus lateralis; and tarsometatarsus length, greatest length from eminentia intercotylaris to trochlea metatarsi III. The area of the sternal carina and the length of the entire manus might better represent functional units than the carinal height and carpometacarpus length, respectively, but the former measurements are rarely available for fossil species, so ones were employed as proxies. In some instances (e.g., in Ptaiochen pau), the sternal carina was virtually absent (height ∼0 mm). In such individuals, the carinal height was arbitrarily set to 0.1 mm to allow log transformation.

FIGURE 1.

Osteological measurements taken in this study. (A) Humerus in caudal view. (B) Ulna in dorsal view. (C) Carpometacarpus in dorsal view. (D) Tibiotarsus in cranial (left) and lateral (right) views. (E) Femur in cranial view. (F) Tarsometatarsus in dorsal view. (G) Sternum in left lateral view. (H) Sternum in cranial view. See text for detailed definition of measurements. Drawings based on Aythya valisineria (USNM 288639).

Multivariate Allometry

The interspecific variation of limb dimensions in Anatidae was investigated with principal component analysis (PCA) from a variance-covariance matrix of log-transformed variables. It is of interest to this study to examine whether there is a general trend of deviation of flightless anatids from volant ones in limb proportion (shape), which is represented by the principal components (PCs) other than the first one (which represents the size component; Jolicoeur 1963a, 1963b, Klingenberg 1996). All statistical analyses were conducted in R 3.2.0 (R Core Team 2015).

In order to estimate PCs while incorporating phylogenetic nonindependence, the phylogenetic principal component analysis (pPCA; Revell 2009) was conducted. Pagel's (1999) lambda statistic was used to adjust the effect of phylogenetic correction and was optimized with maximum likelihood estimation. Hypothetical ancestral states of log-transformed variables at nodes of the phylogenetic tree were reconstructed with the function fastAnc in the package phytools (Revell 2012), which reconstructs ancestral states utilizing Felsenstein's (1985) contrast algorithm (see also Rohlf 2001). The reconstruction was performed for all variables on a tree scaled with the optimal Pagel's lambda, and corresponding pPC scores were calculated. Although only species with phylogenetic information can be included in pPCA, pPC scores for those without phylogenetic information were plotted a posteriori for visual presentation. pPCA was performed with the package phytools (Revell 2012).

The topology of the working phylogeny was taken mainly from Gonzalez et al. (2009), and some taxa not sampled by those authors were grafted according to several other phylogenies (Sorenson et al. 1999, Paxinos et al. 2002, Fulton et al. 2012, Liu et al. 2014, Mitchell et al. 2014). Nodes with incongruent branching orders were collapsed into a polytomy. All branch lengths were first calculated with Grafen's (1989) method and were subsequently scaled using the function chronos in the package ape (Paradis et al. 2004), which estimates a chronogram with the penalized likelihood, using each branch length as the number of substitutions per character, with the smoothing parameter λ = 1 and age constraints taken from previous studies (Gonzalez et al. 2009, Fulton et al. 2012, Mitchell et al. 2014). Because all extinct taxa represented in the tree were quite young in geological age (Late Pleistocene–Holocene), branch lengths were not scaled for their tip ages. The working phylogeny includes 76 modern and 4 fossil species ( Supplemental Material Figure S1). Grafting of missing taxa onto a phylogeny based on taxonomic information is often used in the literature, but it is not employed here because the method is suboptimal in estimating phylogenetic signals of traits and, possibly, evolutionary correlations between traits (Rabosky 2015).

Linear Discriminant Analysis

Classification rules of the linear discriminant analysis (LDA) of skeletal dimensions were constructed from measurement data of modern anatids and subsequently applied to fossil anatids. LDA extracts a linear combination of variables that maximizes the discrimination between each pair of groups, assuming multivariate normal distributions with a common variance-covariance matrix across groups. Group assignments of observations are based on posterior probabilities, relative likelihoods of assignments given prior probabilities for groups. Values of the linear combination can be used to represent classification results, when adjusted so that zero corresponds to the discriminant boundary (so-called discriminant score). In this study, 2 groups (volant and flightless) based on 7 skeletal dimensions were considered, and 93 species means (89 volant and 4 flightless) with complete sets of variables were used to construct rules.

Raw variables were log transformed before the analysis to conform to allometric analyses. Although there is slight evidence of asymmetric distribution of data points, LDA is known to be reasonably robust to the deviation from multivariate normality (McLachlan 1992) and gives a quite stable classification rule (Hastie et al. 2009).

The sample prior proportions (89/93 for volant vs. 4/93 for flightless species) were used as prior probabilities in the analysis, which resulted in adding 3.10 [= log(89 / 4)] for all discriminant scores (not included in the constant term shown in tables). This practice renders inference of flightlessness more conservative than when equal prior probabilities (0.5 vs. 0.5) are used. Such a practice may be justified when observations are assumed to come from groups of different sizes (see also Angielczyk and Schmitz 2014). The extant proportions of Anatidae (∼160 volant vs. 5 flightless) appear to be less realistic than the sample proportions, because the former values are likely to have been biased by selective extinctions of flightless species in the Holocene (e.g., Duncan and Blackburn 2004).

The precision of classification rules was assessed by the bootstrap cross-validation (BCV) error rate estimation (Fu et al. 2005) with 10,000 bootstrap replications. The bootstrap resampling was conducted in a stratified manner (i.e. the subsample number for each group in a bootstrap replicate was fixed as in the original sample). The confidence intervals of discriminant coefficients were constructed from bootstrap percentiles that were constructed in the same manner. In order to examine the effect of intraspecific variation in classification results, all available individual values were evaluated with the constructed rules. When possible, species means with incomplete sets of variables, which were not used to construct rules, were also evaluated.

A statistical classification method that takes phylogenetic correlation between organisms into account has been proposed by Motani and Schmitz (2011). However, its statistical properties and underlying models are not well understood, and there seem to be several methodological issues in the application. Thus, I decided not to employ this method in the present study.

Multivariate Allometry

The pPCA was first performed with the dataset of 7 log-transformed variables for both modern and fossil anatids. The inspection of pPC plots clearly indicated that several fossil species have extreme skeletal proportions and perform as outliers in the estimation of pPCs, masking general patterns of interspecific variation (results not shown). In order to more closely inspect variations between modern volant and flightless species, the dataset including only modern species was analyzed (Table 1 and Figure 2). pPC1 was strongly positively correlated with all 7 variables, explaining 92.0% of the total variance, and hence can be interpreted as the general skeletal size axis, whereas other pPCs can be interpreted as shape axes. pPC2 is negatively correlated with the carinal height and lengths of 3 wing bones and positively with 3 leg bones and, thus, can be interpreted as the contrast between pectoral and pelvic elements. pPC3 can be interpreted as an additional contrast between ulna length vs. carinal height and femur length. Modern flightless species were characterized by high pPC2 scores and low pPC3 scores, indicating that, not unexpectedly, they have relatively small wings for their size. The results of ancestral state reconstructions at nodes differed substantially, depending on whether fossil species were included in the analysis, probably because the fossil species included had extreme skeletal proportions (Figure 2). It is notable that modern flightless species occupy a distinct region of the morphospace from volant ones (Figure 2), especially in the “shape” subspace (i.e. subspace perpendicular to pPC1). Furthermore, the ancestral state reconstruction demonstrated that such characteristic skeletal proportions in flightless anatids had been obtained independently in 2 lineages, Anas and Tachyeres (Figure 2). These results suggest that discrimination of fossil crown-group flightless anatids from skeletal proportions is feasible, as long as the assumption that the form–function relationship is common within Anatidae holds.

TABLE 1.

Results of a phylogenetic principal component analysis (pPCA) including 75 modern anatids. Coefficients of first 3 pPCs and phylogenetic means (p-mean) of log-transformed variables are shown, as well as eigenvalues and percent of phylogenetic variance explained (PVE) by each pPC. Abbreviations of variables: CAR = carinal height; HUM = humerus length; ULN = ulna length; CMC = carpometacarpus length; FEM = femur length; TIB = tibiotarsus length; TMT = tarsometatarsus length.

FIGURE 2.

Phylogenetic principal component (pPC) plots of 7 log-transformed variables, based on modern species alone. (A) pPC1 vs. pPC2. (B) pPC3 vs. pPC2. Solid branching lines connect hypothetical ancestral states at nodes reconstructed with states in modern species, and broken ones are those reconstructed with both modern and fossil species. pPC scores for species without phylogenetic information were also plotted a posteriori. Note that some fossil species have extreme pPC scores, which are shown in insets. Legend: circles = modern volant species; triangles = modern flightless species; crosses = fossil species.

Linear Discriminant Analysis

Full model. An LDA discriminant rule with all the 7 log-transformed variables (full model) clearly separated volant and flightless anatids (Figure 3A). All training cases, or modern species means, were correctly assigned to their original group with posterior probabilities >0.97, except for volant T. patachonicus, which was assigned to the flightless group with a posterior probability of 0.81, resulting in an apparent error rate of 1.08% (= 1/93). Classification results of the leave-one-out cross-validation were identical to those of the original analysis, only T. patachonicus being misidentified. The estimated BCV error rate was 1.13%. With the classification rule, most of 477 individuals of modern species were correctly assigned to their group, except for 2 out of 4 individuals of T. patachonicus. Despite these misclassifications, discriminant scores did not overlap between volant and flightless groups, for both species means and individual values (Supplemental Material Table S1). Altering prior probabilities from the sample values of 0.95:0.05 (= 89:4) to the equal 0.50:0.50 and to 0.97:0.03 (= 160:5) prior probabilities would move the discriminant boundary by 3.10 and −0.36, respectively, and does not qualitatively affect assignments of any modern species means.

FIGURE 3.

Summary of discriminant scores in linear discriminant analysis. (A) The full model. (B) The 4-variable model. Distributions of discriminant scores corresponding to species means of modern anatids are shown as histograms. Distributions of scores calculated for individuals are shown by rugs at the bottom. Flightless species are denoted by darker gray. Discriminant score for the species mean of Tachyeres patachonicus is indicated by a black triangle.

Discriminant coefficients were large positive values for the ulna and carpometacarpus, large negative values for the humerus and femur, and relatively small in magnitude for the remaining 3 variables (carinal height, tibiotarsus, and tarsometatarsus length; Table 2). The 95% bootstrap confidence intervals (CIs) of coefficients for the last 3 variables overlapped zero, whereas those of the other 4 variables did not (Table 2).

TABLE 2.

Coefficients of variables used in linear discriminant analysis. For each of the full and 4-variable models, coefficients are shown with bootstrap 95% confidence intervals in parentheses. Note that the constant does not include the adjustment term for prior probabilities. Abbreviations are defined in Table 1.

Model selection. HAIC was evaluated for all possible combinations of 7 variables that include more than 1 variable (120 models). The results of model selection are summarized in Table 3. The full model ranked 16th with a ΔHAIC value of 6.64 (w_i = 0.009), showing high redundancy among variables. All models with substantial model weights (w_i > 0.01) included ulna length and femur length, along with several or none of the other variables, which suggests the importance of these 2 variables in assessing flight ability in anatids, as confirmed by the Akaike weights and model-averaged scaled discriminant coefficients of variables (Table 4). The carpometacarpus length and humerus length follow the 2 variables, whereas the other 3 variables were much less important for these measures.

TABLE 3.

Results of model selection by HAIC (see text). Doubled negative log likelihood (−2 log L), doubled estimated bias (2 b), ΔHAIC, model weight (w_i), and model formula are shown for those with HAIC smaller than that of the full model. Abbreviations are defined in Table 1.

TABLE 4.

Variable weights and averaged coefficients based on ΔHAIC (see text). Averaged coefficients and those scaled with within-group standard deviation of variables are shown for the ordinary linear discriminant analysis. Abbreviations are defined in Table 1.

Across all models with more support than the full model (Table 3), signs and magnitudes of discriminant coefficients were quite similar, except that the signs were variable for tibiotarsus and tarsometatarsus length. Classification results for all species mean and individual values of both modern and fossil species were also consistent across those models wherever applicable, except for some individuals of T. patachonicus and the species mean of Anas marecula, which generally lie near the discriminant boundary. Models with the smallest HAIC within models including each of 4 variables, 3 variables, and 2 variables are referred to, respectively, as the “4-variable model” (humerus, ulna, carpometacarpus, and femur length), “3-variable model” (ulna, carpometacarpus, and femur length), and “2-variable model” (ulna and femur length) in the following analyses (Table 3).

Size correction. A discriminant rule of size-corrected LDA based on all 7 variables (size-corrected full model) was constructed from scores of PC2 to PC7 extracted from the within-group variance-covariance matrix of the variables as described above. The discriminant rule of the size-corrected full model was quite similar to that in the ordinary LDA, as indicated by the small angle between 2 discriminant coefficient vectors (Table 5). This is because the contrast vector between the group mean vectors in the original space happened to be nearly perpendicular to the within-group PC1 axis in this case. Group assignments in the size-corrected full model were identical to those in the ordinary full model for all species means (only T. patachonicus is misclassified) and for all individuals except one individual of T. patachonicus, which lies near the discriminant boundaries in both rules (Supplemental Material Figure S2A). The estimated BCV error rate was 1.14%.

TABLE 5.

Coefficients of size-corrected models in linear discriminant analysis. For each of the full and 4-variable models, coefficients are shown with bootstrap 95% confidence intervals in parentheses. The angle between the discriminant coefficient vectors in the ordinary and size-corrected LDAs (degrees) is shown in parentheses at the heading. Note that the constant does not include the adjustment term for prior probabilities. Abbreviations are defined in Table 1.

As for the full model, the size-corrected LDAs corresponding to reduced models were quite similar to corresponding ordinary LDAs, as indicated by small angles between discriminant vectors (Table 5). For the 4-variable model, all group assignments of modern species mean and individual values in size-corrected LDA were identical with those in the ordinary LDA ( Supplemental Material Figure S2B). The estimated BCV error rate was 0.97% for the size-corrected 4-variable model. Additional information on size-corrected LDA is given in  Supplemental Material Tables S3– S5.

Inference of fossil species. LDA discrimination rules constructed with species means of modern anatids were applied to species means and individuals of fossil anatids to infer the presence–absence of flight ability. For species means, discriminant scores and posterior probability for the assignment of flightless group in the full, 4-variable, 3-variable, and 2-variable models are shown in Table 6. When available for prediction, species means of Cnemiornis calcitrans, C. gracilis, Branta hylobadistes, B. rhuax, Chenonetta finschi, P. pau, Thambetochen chauliodous, Chendyes lawi, Shiriyanetta hasegawai, Anas chathamica, and “Annaka Short-winged Swan” were invariably assigned to the flightless group with high posterior probabilities (∼1.00), whereas Mergus milleneri was always assigned to the volant group with high posterior probabilities (>0.99). The species mean for A. marecula was assigned to the volant group in the 3-variable model, but to the flightless group in the others, with relatively ambiguous posterior probability values.

TABLE 6.

Classification results of species mean values of fossil anatids. For each species and each model, a discriminant score calculated from species means and posterior probability for the assignment to the flightless group is shown. A positive discriminant score and a posterior probability <0.50 indicate assignment to the volant group, whereas a negative score and a posterior probability >0.50 indicate assignment to the flightless group. All values were adjusted for the sample prior probabilities of 89/93 vs. 4/93. Posterior probabilities >0.995 are shown as “∼1.00.”

Only limited sets of variables are available for Bambolinetta lignitifila, Cayaoa bruneti, and Chendytes milleri, precluding application of most highly supported models. The flight ability of these species was assessed by ad hoc models with the highest HAIC support within sets of applicable models (Table 7). As a result, B. lignitifila was assigned to the volant group with a discriminant score well within the range of the volant group's score (4.12) and a high posterior probability (0.98). Cayaoa bruneti and Chendytes milleri were assigned to the flightless group with high posterior probabilities (∼1.00).

TABLE 7.

Classification results for Bambolinetta lignitifila, Cayaoa bruneti, and Chendytes milleri, including variables available, the model formula of the ad hoc model with highest HAIC support within those applicable, corresponding model weight (w_i′), discriminant score calculated from species means, and posterior probability (pp) for the assignment to the flightless group. For further information, see Tables 1 and 6.

For species for which more than one associated skeleton was available, results of the group assignment for individual values are summarized in Table 8. Individuals of Branta rhuax, P. pau, and T. chauliodous were consistently assigned to the flightless group. One individual of B. hylobadistes was, unlike other individuals of the same species, assigned to the volant group in the 4-variable model (Table 8).

TABLE 8.

Classification results for individuals of fossil anatids. For each species and each model, the range of discriminant scores is shown, with individual numbers available for the prediction with the model in parentheses. All values were adjusted for the sample prior probabilities of 89/93 vs. 4/93.

Resampling experiment. Resampling experiments were conducted for the full, 4-variable, 3-variable, and 2-variable models to examine possible effects of varying numbers of variables. First, resampling was performed by randomly choosing a single individual value independently for every element to form a hypothetical distribution of species “means” for each species, from which resampled distributions of discriminant scores were obtained (Figure 4; Supplemental Material Table S7). The resampled distributions were widely dispersed and in some cases spanned the discriminant boundary, although the dispersion tends to be smaller for models with fewer variables (Table 9). Dispersion of the resampled distributions of discriminant scores also varies among species and is relatively small for Aythya marila, Netta rufina, Melanitta deglandi, and Tachyeres spp., whereas it is quite large for Biziura lobata.

FIGURE 4.

Distribution of discriminant scores from resampled experiments. (A) The full model. (B) The 4-variable model. (C) The 3-variable model. (D) The 2-variable model. Discriminant scores were calculated from a resampled distribution; a single individual value for every element was randomly chosen 10,000 times for each species. Distributions of resampled species means are shown as box plots, whereas original individual values for each species are indicated by gray circles. Vertical positions stand only for graphical purpose. Note that Chendytes lawi has extreme discriminant scores, which are shown in insets. Abbreviations for species labels: An. c. = Anas chathamica; An. m. = Anas marecula; As. r. = Anser rossii; As. a. = Anser albifrons; Ay. m. = Aythya marila; Bi. l. = Biziura lobata; Br. s. = Branta sandvicensis; Cd. l. = Chendytes lawi; Ch. f. = Chenonetta finschi; Cn. c. = Cnemiornis calcitrans; Cn. g. = Cnemiornis gracilis; Cy. c. = Cygnus columbianus; Ml. d. = Melanitta deglandi; Mr. m. = Mergus merganser; Ne. r. = Netta rufina; Sh. h. = Shiriyanetta hasegawai; Ta. b. = Tachyeres brachypterus; Ta. pa. = Tachyeres patachonicus; Ta. pt. = Tachyeres pteneres.

TABLE 9.

Scaled standard deviation of discriminant scores from resampling experiments. For each species, a single individual value was randomly chosen for every element 10,000 times, and corresponding discriminant scores were calculated for each model. For comparison among models, standard deviations of resampled discriminant scores were scaled for within-group standard deviation of discriminant scores corresponding to original species means in each model. Only results for ordinary LDA are shown, because those for size-corrected LDA were almost identical to them. For other statistics, see  Supplemental Material Table S7.

Discriminant scores for fossil species known only from unassociated bones were compared with the resampled distributions. More specifically, resampled distributions were treated as null distributions for the test of flightless hypothesis of such fossil species. In all models considered, discriminant scores for species means of Cnemiornis calcitrans, C. gracilis, Chendytes lawi, and Shiriyanetta hasegawai were much smaller than the resampled distributions of modern volant species, indicating that small discriminant scores for those species are unlikely to have been extracted from populations similar to volant species by chance alone (Figure 4). By contrast, discriminant scores for Anas chathamica and Chenonetta finschi, which unambiguously lie in the region assigned to the flightless group in most models applicable (the full model was not applicable to these species), were within the resampled range of volant Biziura lobata in the 4-variable and 3-variable models (Figure 4B, 4C). In order to more closely examine these species, another set of resampling was conducted for each species, by resampling as many individual values for every element as were available for those species and taking their means, which is expected to reflect a similar degree of precision of estimation of species means in the fossil species. Results with the 4-variable model for A. chathamica and C. finschi, whose data were taken from Williams (2015b) and Worthy (1988), respectively, are shown in Figure 5, along with those with ad hoc models for Chendytes milleri and Cayaoa bruneti. For C. finschi, the sample size in this study was too small to allow further comparisons; hence, data reported by Worthy (1988) for 2 localities (Graveyard layer 2 and Martinborough) were examined separately instead. For each species, perhaps except for C. finschi from Graveyard layer 2, the species mean was clearly lying below resampled distributions of any modern volant species compared, indicating that the species mean is unlikely to have come from distributions similar to modern volant species by sampling error due to the lack of associated skeletons.

FIGURE 5.

Distribution of discriminant scores from resampled experiments for Anas chathamica (A) and Chenonetta finschi (B) based on the 4-variable model, and for Chendytes milleri (C) and Cayaoa bruneti (D) based on ad hoc models. Discriminant scores were calculated from the simulated distribution for each fossil species, which is a resampled distribution of 10,000 mean vectors of randomly chosen individual values of the same number as available for every element of that species. For C. finschi, values from 2 localities (Graveyard layer 2 and Martinborough) studied by Worthy (1988) are shown, and the smaller set of sample sizes was used. Distribution of resampled species means are shown as box plots. Vertical positions are for graphical purposes only. Note that Chendytes milleri has extreme discriminant scores, which are shown in insets. Abbreviations for species labels: Ca. b. = Cayaoa bruneti; Cd. m. = Chendytes milleri (see Figure 4 for other abbreviations).

Results for inference of fossil species and resampling with size-corrected LDAs were quite similar to those with ordinary LDAs (see  Supplemental Material Figures S3 and  S4 and  Tables S8– S12).

Applicability and Limitation of the Method

Modern anatids are diverse both taxonomically and morphologically. Some species can vertically take off from the water or from the land, whereas others need to taxi for a long distance on the water before taking off (Raikow 1973). In addition, although most anatids are suited to foraging on the water surface, some specialize in feeding on the ground, and some others regularly dive into the water column to obtain food items (e.g., Hughes and Green 2005). Given the morphological and ecological diversity of modern anatids, one might doubt whether there could be a universal criterion to assess flight ability in fossil anatids. However, it is notable that the 2 lineages of modern flightless anatids, the diminutive A. aucklandica and the large-bodied, diving Tachyeres spp., appear to have undergone similar modifications of limb skeletal proportions, as demonstrated by the ordination of pPCA and ancestral state reconstruction (Figure 2); in short, the flightless species have relatively short wing elements, especially distal ones, compared to their ancestral states. In addition, a few sedentary insular species with disparate ecologies, including Branta sandvicensis and Mergus australis (see Miller 1937, Livezey 1989), appear to have shifted in a similar direction. These facts suggest the presence of a certain component in limb proportion that responds to the evolution of reduced flight ability, which justifies the assessment of flight ability in fossil anatids from limb proportion.

LDA models constructed from skeletal measurements in modern anatids gave reasonable separation and classification results for the volant and flightless groups. In all models that attained substantial support from HAIC (w_i > 0.01) as well as in the full model, all species means were correctly assigned to their original group except for T. patachonicus, whose discriminant scores tend to be distributed around discriminant boundaries. Misclassifications of T. patachonicus were not completely unexpected, however, since that species is known to have the highest wing loading among volant anatids and to be “partially flightless,” in that some of the heaviest males are (at least temporarily) not capable of becoming airborne (Humphrey and Livezey 1982). Tachyeres patachonicus would have been under selection for large body size and small wings, which resulted in its specialized morphological traits (Livezey and Humphrey 1986). Apart from the interpretability of misclassified cases, it should be noted that the inference based on LDA is only in a statistical sense; the actual functional threshold for flightlessness, if one exists, may lie somewhere else than the discriminant boundary. Although the separation of ranges of discriminant scores between groups may be satisfactory for a descriptive purpose, it can raise a practical problem in prediction when a discriminant score for a fossil anatid falls between the 2 ranges. Although the discriminant rule would strongly suggest that the species would be flightless in such cases, the inference should be made with caution since the position of the actual functional boundary is not clear. The inference of the actual boundary is inevitably prohibited by the scarcity of the training sample or by the low diversity of modern flightless anatids that form the basis for construction of discriminant rules. Nevertheless, the discriminant rules presented here incorporated most of the information available at present and would represent an objective basis for inferring the presence–absence of flight ability in fossil anatids.

In the models considered above, discriminant scores for individuals were distributed around corresponding species mean values and mostly assigned to the same group with them (see also Figure 4), which suggests that classification results for individuals of fossil anatids would be largely consistent within species, except when they lie near discriminant boundaries. For fossil species known only from isolated bones, the resampling procedure described above can be used to assess the reliability of assignments based on species means by evaluating the extent of possible errors due to intraspecific variation.

Classification results for both species means and individuals in models with smaller HAIC values than the full model were almost consistent with those in the full model, indicating that little information is lost by dropping some variables (carinal height, tibiotarsus length, and tarsometatarsus length) and that these models are more preferable than the full model on the grounds of parsimony. Indeed, the estimated BCV error rates for models with smaller HAIC values (range: 0.86–1.09%, mean = 0.96%) were slightly lower than that for the full model (1.13%), which suggests that the reduced models may perform better in prediction. One of the largest advantages of reduced models, along with avoidance of overfitting to training sample, is that the inferences are less likely to be affected by sampling errors (intraspecific variations) in estimating species means when only a small number of isolated specimens are available for measurement, as is true for many fossil anatids. This is because errors in estimating the species mean for each element are amplified by multiplications of discriminant coefficients and are summed in calculating discriminant scores. When more variables are included, more errors are included in the calculation. Such inflations of discriminant scores were exemplified by the resampling experiments, in which dispersions of discriminant scores were larger in models with more variables (Figure 4 and Table 9;  Supplemental Material Figure S3 and  Tables S7 and  S10).

One peculiarity observed in the resampling experiment was that dispersions of resampled distributions of discriminant scores were quite large for Biziura lobata compared to other species considered (Figure 4 and Table 9), whereas those for original individual values were not so large. It is likely that this is largely due to the pronounced sexual dimorphism of body size in B. lobata, which is among the most extreme within Anatidae (Livezey and Humphrey 1984, McCracken et al. 2000); when the intraspecific variation in body size and, hence, in morphological variables is large, the dispersion of a linear combination or proportion of resampled individual measurements is expected to be large.

Some previous studies have explored correlations between morphological and ecological traits by explicitly controlling for the effect of size (e.g., Simons 2010, Hughes 2013). In the present study, size-corrected LDAs were conducted by removing allometric size represented by PC1 of the within-group variance-covariance matrix for each variable set. The resultant classification rules of size-corrected LDAs were quite similar to those of the corresponding ordinary LDAs, largely because the contrast vector between the 2 group means happened to be nearly perpendicular to the size axis. The motivation for the allometric size correction was that skeletons may change in proportions while varying in size in order to retain similar performance, which is expected to be reflected in the allometric scaling of skeletal proportions. Unfortunately, however, current observations are likely to be largely constrained by the availability of the modern sample, and it is not clear whether the allometric scaling trends observed in the sample and the apparent threshold perpendicular to it hold for birds with larger or smaller body size than observed in the current sample, which spans almost the entire range of body-size variation in modern Anatidae. For example, Northcote (1982) estimated the body mass of an extinct swan (Cygnus falconeri) as ∼16 kg from the dimensions of tarsometatarsus, assuming elastic similarity to recent C. cygnus and considered that such a large anatid might have a limited flight ability or even be flightless. In the absence of appropriate modern species of comparable size, it is not clear whether such an assertion based solely on body mass can be justified or whether statistical classification methods like LDA constructed from modern samples can be applied to such a large species. Another problem is exemplified by a small extinct species, Anas marecula (see below). Inferences of flight ability based on estimation of kinematic parameters of flight (e.g., Campbell and Tonni 1983, Pennycuick 1996, Ksepka 2014) have the same fundamental difficulties of extrapolations, among others (Alexander 2003). Issues concerning size-related variation in flight ability and extrapolations clearly require further theoretical and empirical studies.

Functional Interpretation

Functional interpretations of discriminant coefficients and relative importance of variables would be useful in assessing functional bases underpinning discriminant rules. Discriminant coefficients in this case correspond to hypothetical marginal effects of changes in skeletal dimensions while other dimensions are held constant. In almost all models, the discriminant axes represent a contrast between relative sizes of wing and leg elements, and flightless species were characterized by relatively small wing elements (especially the ulna and carpometacarpus) and large leg elements, as expected from the results of PCA. It might be less intuitive that the coefficients for humerus length are negative (in contrast to the other wing elements) whenever it is included in models with ulna length (see also Table 2 and  Supplemental Material Table S2). Statistically speaking, such a contrast in signs of coefficients for 2 variables with apparently similar roles might arise as an artifact due to the high correlation between the humerus length and ulna length (i.e. when they are redundant). However, the fact that many models including both of these variables attain relatively large supports from HAIC (Table 3) suggests that they are not completely redundant in discriminating the 2 groups. Rather, it is more likely that the contrast between humerus and ulna lengths is associated with a functional constraint in skeletal proportions of wings in volant birds. The ratio between the humerus (upper arm) and the ulna (forearm) is known to be correlated with various morphological traits relevant to flight mechanics, and to show a strong phylogenetic signal (Rayner and Dyke 2003, Nudds et al. 2004, 2007). Rayner and Dyke (2003) postulated that a relatively short humerus is beneficial in reducing wing inertia during flapping flight, and that there are constraints in the position of the wrist joint in the upstroke and various morphological traits in the shoulder joint complex that are affected by a modification of the ratio. It is likely that the ratio is under such constraints imposed by flapping flight in volant anatids, whereas it is free to vary in flightless anatids that have been released from those constraints. Relative shortenings of distal wing elements, especially the ulna and radius, have been repeatedly observed in flightless anatids (Livezey and Humphrey 1986, Livezey 1990, 1993c), as well as in other avian groups (e.g., Gadow 1902, Livezey 1992, 1993a, Middleton and Gatesy 2000, Smith 2011). Given that skeletal proportions of modern flightless anatids appear to have been shifted in similar directions, it is conceivable that functional demands other than flight (e.g., combating or steaming) or developmental constraints may affect skeletal proportions of flightless species. At present, however, little information is available to further discuss effects of such factors on avian skeletal proportion.

Among the 3 variables taken from leg bones (femur, tibiotarsus, and tarsometatarsus lengths), femur length is by far the most important by means of both model weights and magnitudes of scaled coefficients (Table 4). It can be concluded that the tibiotarsus and tarsometatarsus lengths are largely redundant in discriminating the 2 groups when femur length is taken into account. It is notable that femur length has a different scaling pattern than other limb bones and shows strong negative allometry in relation to the length of the other 2 leg bones (Table 1). Zeffer et al. (2003) found that femur length is the most strongly correlated with body mass among the lengths of the 3 leg elements (both log transformed). The avian femur in general is tightly connected to the trunk and plays an important role in placing the horizontal (anteroposterior) position of the center of mass above the knee joint (Campbell and Marcus 1992, Hertel and Campbell 2007). Given these facts, it is not surprising that the femur is a better indicator of body size than the 2 distal leg bones and acts in the discrimination of the 2 groups as such.

The height of the sternal carina was expected to be correlated with flightlessness through the amount of major flight muscles (mm. pectoralis et supracoracoideus). However, it turned out to be one of the least important variables in discriminating the 2 groups (Table 3). Although there was a tendency for flightless anatids to have a relatively small sternal carina compared to body mass, the variation within each group was so large as to obscure the difference between the groups. For example, the sternum and pectoral muscles of flightless Tachyeres species are not quite diminutive because the wings are utilized in other activities than flight, including underwater diving and combating (Livezey and Humphrey 1986). Among volant anatids, the carina is relatively small in B. lobata, perhaps as part of the extreme streamlining of the body related to the specialization for underwater diving (see also Raikow 1970). Cygnus spp., which are also characterized by relatively small carinae, have specialized sterna, including the intrusion of the trachea (e.g., Johnsgard 1961, Humphrey and Clark 1964) and the modification of attachment sites of pectoral muscles (see Woolfenden 1961, Matsuoka and Hasegawa 2007). Such variations due to specializations not directly related to flight ability might affect the morphology of the sternum in Anatidae.

Assessment of Flight Ability in Fossil Anatids

Discriminant rules constructed in this study can be used to infer the presence–absence of flight ability in fossil anatids, with the limitations mentioned above. Cnemiornis calcitrans and C. gracilis from the the Quaternary of New Zealand, Chendytes lawi and C. milleri from the Quaternary of California, and Shiriyanetta hasegawai from the Pleistocene of Japan were assigned to the flightless group (Figure 4 and Table 6). Indeed, Chendytes is characterized by diminutive wing bones (Howard 1947, Livezey 1993c), perhaps more extremely so than any other anatids except Chelychelynechen quassus, although known wing bones are rather fragmentary for the latter species (Olson and James 1991).

Among the Quaternary anatids from Hawaii, Branta rhuax, P. pau, and Thambetochen chauliodous were assigned to the flightless group (Tables 6 and 8). Results are more ambiguous for another Hawaiian anatid B. hylobadistes from Maui Island, which has been described as “at best a weak flyer” (Paxinos et al. 2002:1399). Although it was assigned to the flightless group, discriminant scores for the species mean fall between the ranges of the 2 groups in some models (Table 6 and  Supplemental Material Table S1). It seems likely that Hawaiian Branta populations were “at a stage in the evolution of flightlessness” as suggested by Olson and James (1991:43). More detailed investigation is required for this species, possibly incorporating variations among individuals from various localities (see Olson and James 1991).

The result for Anas marecula from the Holocene of Amsterdam Island was somewhat problematic. Discriminant scores were distributed around the discriminant boundaries (Figure 4 and Table 6). Although Livezey (1993b) and Olson and Jouventin (1996) considered the species to be almost certainly flightless from the skeletal proportion and the morphology of the sternum, the skeletal proportions considered in this study were found to be inconclusive in confirming the presence–absence of flight ability in this species. This result was largely attributable to the fact that the distal wing elements are relatively long compared to the humerus, despite the relatively small sternum and wing elements compared to leg elements. One perplexing factor is that the species is smaller than any modern anatids (Figure 2; Olson and Jouventin 1996); it is not certain whether the discriminant boundary constructed in this study can be extended beyond the observed range of body size in the modern sample.

Chenonetta finschi from the Quaternary of New Zealand has been considered as either possibly flightless (Holdaway et al. 2002) or facultatively flightless (Worthy and Olson 2002), but most authors remained inconclusive about whether the species was flightless. On the basis of measurements collected during the present study from a limited number of specimens from the Earnscleugh Cave, stored at the Natural History Museum (London, UK), the species was assigned to the flightless group. Worthy (1988, 1997) gave measurements from a large number of specimens collected in various localities and demonstrated that the species have undergone a reduction of wing elements during the period from 14,000–20,000 yr BP to 1,000–2,000 yr BP. Results using some measurement data given by Worthy (1988) for this species from Graveyard layer 2 (11,000–14,000 yr BP) and Martinborough Cave (∼1,500 yr BP) indicate that C. finschi was flightless (Figure 5), at least after the reduction of pectoral elements in that period.

Anas chathamica from the Holocene of Chatham Islands has been described as either flightless (Millener 1999, Williams 2015b) or not (Worthy and Holdaway 2002, Mitchell et al. 2014). The discriminant rules and resampling experiments in the present study assigned it to the flightless group (Figure 5 and Table 6). Mergus milleneri from the Chatham Islands was first stated as being flightless by Millener (1999), but Williams et al. (2014) concluded that the species was volant. Based on measurements given by Williams et al. (2014) and discriminant rules constructed in our study, there seems to be little evidence suggesting the flightlessness of this species (Table 6).

Bambolinetta lignitifila, an anatine of uncertain affinity known from an incomplete skeleton from the Miocene of Italy, was redescribed by Mayr and Pavia (2014). They considered the relatively shortened forearm of the species to be indicative of reduced flight ability, and the flattened wing bones as a possible specialization for wing-propelled diving. From the estimated lengths of wing bones given by them, the species was assigned to the volant group in ad hoc LDA models with dimensions of wing bones (Table 7). At least, it can be concluded that the proportion of wing skeleton alone provides only ambiguous evidence for the reduced flight ability of B. lignitifila.

Cayaoa bruneti, known from isolated bones from the Miocene of Argentina, was considered to be flightless and adapted to foot-propelled diving by Noriega et al. (2008). The loss of flight was considered an important factor that enabled the species to attain extremely thick bone walls (De Mendoza and Tambussi 2015). In the present study, the flightless condition of this species was confirmed by an ad hoc model that includes carpometacarpus and femur length (Table 7). The resampling experiment showed that such an extreme proportion was unlikely to be obtained from a volant anatid (Figure 5D).

The “Annaka Short-winged Swan,” a large unnamed anatid of uncertain affinity from the Miocene of Japan, has quite peculiar skeletal proportions and thick bone walls (Matsuoka et al. 2001, 2004). Analysis of the measurements given by Matsuoka et al. (2004) confirmed the flightlessness of that species.

Conclusion

Discriminant rules to assess the presence–absence of flight ability from skeletal measurements in fossil anatids were constructed and were applied to selected fossil anatids. Variable selection based on HAIC gave high support to models including ulna length and femur length, suggesting the usefulness of these variables in discriminating the 2 groups. Apart from the merit of avoiding overfitting, variable selection procedure appears to be useful in reducing the undesired effect of intraspecific variation in inferring ecological traits from limb proportions when only a small number of isolated bones are available for measurements, as is true for many fossil birds.

The discriminant rules of LDA described here can easily be applied to a future observation by calculating its discriminant score—that is, the sum of natural-log-transformed measurements multiplied by corresponding discriminant coefficients (plus the constant term) as shown in Tables 2 and 5 and in  Supplemental Material Tables S2 and  S4. Adjustment for prior probabilities can be made by adding the log-transformed value of the ratio of prior probabilities for volant over flightless groups (3.10 in the present study). A positive score indicates classification to the volant group, and a negative one indicates classification to the flightless group. The use of a (size-corrected) 4-variable model is recommended, although ad hoc models may be required when some of the variables are lacking. Custom-made analyses can be conducted with R codes and data available in the  Supplemental Material.

Although a simple dichotomy between volancy and flightlessness might seem too simplistic given various transitional cases toward flightlessness in Anatidae (Humphrey and Livezey 1982, Worthy 1988, Guillemette and Ouellet 2005a, 2005b, Fulton et al. 2012), the classification rules given in the present study would at least form an objective basis for ecological and evolutionary studies on fossil anatids. Nevertheless, it should be remembered that inferences based on rules constructed from modern samples are generally constrained by the availability of modern samples and that one should be cautious in extrapolating the rules to fossil species that show departures from morphological variations observed in modern species. At present, it is not certain whether the discriminant rules constructed from modern anatids can be applied to other avian groups such as Rallidae. It might be possible that various patterns of skeletal proportions are realized to perform a similar functional performance of flight (many-to-one mapping; Wainwright et al. 2005). Patterns of correspondence between form and function in the avian skeleton, and to what extent phylogenetic signals may affect them, are under active research (Stoessel et al. 2013, Wang and Clarke 2014). A few methods have been proposed to infer functional signals in morphology while incorporating phylogenetic signals (Motani and Schmitz 2011, Cooper et al. 2014), but much remains to be explored in regard to the theoretical and practical aspects. In order to further explore morphological aspects of avian flightlessness and other ecological traits in general, more empirical and theoretical investigations into biomechanics and evolutionary patterns are required.