Feather evolution following flight loss in crown group birds: relaxed selection and developmental constraints

Feathers are complex structures exhibiting structural/functional disparity across species and plumage. Flight was lost in >30 extant lineages from ~79.58 Ma–15 Ka. Effects of flight loss on senses, neuroanatomy, and skeletomusculature are known. To study how flightlessness affects feathers, we measured 11 feather metrics across the plumage of 30 flightless taxa and their phylogenetically closest volant taxa, with broader sampling of primaries across all orders of crown birds. Our sample includes 27 independent flight losses, representing nearly half of extant flight­less species. Feather asymmetry measured by barb angle differences between trailing and leading vanes decreases in flightless lineages, most prominently in flight feathers and weakest in contour feathers. Greatest changes in feather anatomy occur in older flightless lineages (penguins, ratites). Comparative methods show that many microscopic feather traits are not dramatically modified after flightlessness compared to body mass increase and relative wing and tail fan reduction. Changes involved with greater vane symmetry show stronger shifts, however. Relaxing selection for flight does not rapidly modify feather flight adaptations, apart from asymmetry. Developmental constraints and relaxed selection for novel feather morphologies may explain some observed changes. Macroscopic changes to flight apparati (skeletomusculature, airfoil size) are more evident in recently flightless taxa and could more reliably detect flightlessness in fossils, with increased feather symmetry as a potential microscopic signal. We observed apical modification in later stages of feather development (asymmetric displacement of barb loci), while mor­phologies arising during early developmental stages are only altered after millions of years of flightlessness.

Methods
Taxon Sampling.

There are ~64 species of extant flightless birds, ~28% of which are penguins (Sphenisciformes). Ecologically, ~59% of flightless species are terrestrial and ~41% are semiaquatic (Roots 2006). We examined bird skins at the Field Museum of Natural History (FMNH) and American Museum of Natural History (AMNH). Our sample (83 total taxa, 307 feathers) contained 30 taxa (nine orders) fully flightless as adults (Palaeognathae, Gruiformes, Psittaciformes, Charadriiformes, Anseriformes, Podicipediformes, Suliformes, Sphenisciformes, Passeriformes), nearly half of which are rare/threatened (IUCN 2022). Our sample is ~57% terrestrial and ~43% semiaquatic, resembling the ecological breakdown among all living flightless species.

We reduced sampling of flightless species within non-volant lineages to avoid redundancy (Supplementary Material, section 1): our dataset is only 10% penguins, while broadly covering major penguin lineages (Spheniscus, Eudyptes, Aptenodytes [Ksepka et al. 2006]). Only one species from each of the five living genera (13 total species) of ratites (i.e., flightless Palaeognathae) is included because they represent only four independent losses of flight (Yonezawa et al. 2017): 1) Struthio camelus, 2) Rhea americana, 3) Dromaius novaehollandiae/Casuarius unappendiculatus, and 4) Apteryx australis.

Six of our flightless taxa are recently extinct (Porzana sandwichensis, Porzana palmeri, Gallinula nesiotis, Xenicus lyalli, Podilymbus gigas, Pinguinus impennis). One recently extinct species was poor-flying/‘incipiently flightless’ (Mergus australis).

For each flightless/‘incipiently’ flightless taxon, we included its phylogenetically closest volant taxon for sister-taxon comparisons. For penguins, we included three volant species of tubenoses (Procellariiformes), representing major lineages within that clade (Oceanites oceanicus, Pelecanoides urinatrix, Diomedea immutabilis).

We expanded sampling of crown birds to include at least one species from each commonly recognized order (Reddy et al. 2017) All these broadly sampled species are volant, except for Mesitornis unicolor, whose flight capabilities are unknown.

Data Acquisition.

Macro- and microscopic measurements were made on skins of size/proportions, densities, and positions of various skeletal, scale, and feather traits (Supplementary Data S1), supplemented with values from the literature (e.g., species body mass was frequently taken from Dunning [2008]). Chosen variables (Supplementary Material, section 3) were considered relevant from prior studies (e.g., Cubo & Arthur 2001; Feo et al. 2015), sufficient to describe anatomy, and practical to measure. Macroscopic traits included body mass, wing length from the wrist to tip of the longest primary, and tail fan length from the posterior end of pygostyle to tip of the longest rectrix. Data, scripts, and archiving are available at https://github.com/paleomitchelljs/FeatherEvo and in Supplementary Data S4.

Equipment. Macroscopic traits were measured with analog calipers, digital calipers, and measuring tape. Microscopic observations of feathers utilized a Dino-Lite Edge AM73915MZTL (10X~140X; 5.0MP; USB 3.0) Digital Microscope with calibrated scalebar (supplemented with a physical scale bar) with a maneuverable Dino-Lite MS36A-C2 adjustable precision mount and clamp. Measurements were taken from micrographs using Adobe Photoshop (version 20.0.6).

Feathers sampled. For clades containing flightless taxa, five feathers across the body were measured: longest primary remex, dorsal-most tertial remex, dorsal-most rectrix, and dorsal and ventral contours from roughly the middle of the torso. For orders without flightless taxa, only the longest primary remex was measured, assuming that effects of flight loss are strongest in them (a hypothesis ultimately supported). For internal consistency, we define ‘trailing’ vanes of contour feathers as those more medial and ‘leading’ vanes as those more lateral, even though these terms are inapplicable to feathers that do not generate lift.

Each feather was manipulated by hand for maximum exposure from the plumage. Paper was placed underneath for contrast. Barbs were exposed using a metal pick and sewing pins.

Observations of each feather. Observations were made at apical, middle, and basal regions along the exposed extent of the rachis (Supplementary Data S1). We focused on the middle region because the apex is often damaged in museum specimens and exhibits high asymmetry between comparable feather positions on left and right sides of an individual. Basal regions would sometimes include plumulaceous regions, depending on how much the feather could be exposed.

Feather metrics of greatest interest are shown in Fig. 1A (for kiwi, a barb pair with at least one well-developed barb was measured [Fig. 1B; Supplementary Material, section 2]). Barbules and mid-barb width were measured only on the trailing barb, since the trailing vane has a greater area than the leading vane on flight feathers, relating to the wing’s lift surface. To control for overall feather size and position along the rachis’ length, linear feather metrics (rachis width, mid-barb width on trailing barb, leading barb length, trailing barb length, and distal barbule length on trailing barb) were divided by the length from the apical most barb node to the measured barbs. Although some groups (e.g., penguins, ratites) may contain outliers, body mass and feather length positively correlate in flying species (Sullivan et al. 2017). In our dataset, feather size (i.e., rachis length at middle region) shows greater correlation to taxon body mass than does feather asymmetry (i.e., difference in trailing and leading barb angle at middle region), suggesting that body size does not confound our conclusions (Supplementary Material, section 9; Supplementary Data S3). Proximal barbule length was not examined due to difficulties in manually unfurling tightly overlapping barbules in many taxa.

Error in barb orientation can result from curvature at the base of barbs or manipulation (Feo et al. 2016). Barb angle differences were used to measure feather asymmetry at a given barb pair along the rachis, controlling for distortion at that position along the rachis or for variation in measurement. Feather asymmetry equaled the angle of the trailing barb subtracted by the angle of its paired leading barb – consistent with developmental models of feather asymmetry (Feo & Prum 2014). Lower values indicate greater symmetry. Trailing barb angles are typically greater than leading barb angles in volant birds (Feo et al. 2016), and greater asymmetry is indicated by larger positive differences in barb angles. Other studies often measure feather asymmetry using the width of the leading versus trailing vanes perpendicular to the rachis. As a fractal structure, feather asymmetry is also reflected in other metrics (e.g., barb densities or barb lengths on leading versus trailing vanes).

Filament densities were measured as the number of filaments per unit distance (i.e., number of barbs or barbules per mm along rachis or barb, respectively). This approach does not consider scaling (i.e., similarly shaped feathers can vary in density due to size alone), but is appropriate with feather development (Prum 1999). It assumes larger feathers do not develop isometrically, but rather allometrically (consistent with models by Maderson & Alibardi 2000; Prum & Brush 2002; Sawyer et al. 2003; Prum 2005; Lin et al. 2006; Chernova & Kiladze 2019). Unscaled filament density values reflect the rate of branching during development. To account for size confounding, we also analyzed scaled densities of barbs and barbules in primaries (Supplementary Material, section 8), where density was multiplied by rachis width at that position, converting filaments/mm into filament counts.        

Analysis.

We triangulated (Munafò & Smith 2018) using simple and sophisticated statistical analyses with increasing phylogenetic control, rerun under varying conditions (Supplementary Material, sections 4-7; Supplementary Data S2-S4). We compare relative effect sizes, rather than solely relying on arbitrary thresholds of statistical significance (Amrhein et al. 2019). Given available sample sizes, we made broad comparisons (e.g., volant versus flightless, recent versus ancient divergences, inter-order differences, terrestrial versus semiaquatic). Semiaquatic ecologies were coded 0 to 5 inspired by Fish (2016): 0, terrestrial; 1, adaptive surface hindfeet paddling (beyond wading or ability to swim); 2, surface wing paddling (i.e., ‘steaming’) with surface hindfeet paddling; 3, underwater hindfeet rowing; 4, asymmetrical underwater flight; 5, symmetrical underwater flight. However, fine parsing leads to small sample sizes per category, so separate analyses lumped categories 1–5 into a single semiaquatic ecological designation.

Analyses were performed using R coding language (RStudio version 1.3.1056). We started with sister taxon comparisons between volant taxa and their phylogenetically closest flightless relative by subtracting the value of the volant taxa (e.g., a density) from that of their sister flightless taxa: Δ = MetricFlightless – MetricVolant. For penguins, volant values were averaged from three tubenoses (Oceanites oceanicus, Pelecanoides urinatrix, Diomedea immutabilis). Data were normalized during principal component analysis (PCA).

Data were not log-transformed when feathers without barbules/barbs would yield undefined values (i.e., PCA of linear metrics of rachis, barbs, and barbules, or shifts in density and barb angle symmetry), allowing for inclusion of taxa lacking barbs/barbules (i.e., some ratites [Table 1]), although the data was still normalized during PCA. PCA here reveals positions of taxa relative to each other within the morphospace, not absolute distances between taxa.

Phylogenetic comparative methods (PCM). A series of phylogenetically controlled comparisons tested how rapidly traits changed after flight loss (Supplementary Data S4). Using the phylogeny from Jetz et al. (2012) combined with the Gruiformes phylogeny from Garcia-Ramirez & Matzke (2021), we added recently extinct taxa and dates to their appropriate branches at a stochastically determined time. We repeated this procedure 100 times, creating a set of trees for 82 living (cassowary dropped because 0 values of various primary metrics cannot be log-transformed) and recently extinct taxa. We then employed stochastic mapping (Revell 2012) to generate timings for flight loss along the tree. Transition rates were set to be unidirectional (flight could be lost but not re-evolved). The ancestral state for crown birds was fixed as volant.

Using this set of calibrated trees, we compared rates of change across a suite of macroscopic traits (supplemented from literature) and microscopic traits at the middle region of the exposed primary remex.

We compared the overall distributions of these macroscopic and microscopic traits between flightless and volant taxa without phylogenetic correction. We then compared 14 volant and flightless sister pairs in a strict sense. We looked for differences in rates between macro and microscopic traits by fitting one- and two-rate Brownian motion models using the R package OUwie (Beaulieu et al. 2016) to each of our 100 stochastically mapped phylogenies for each mean-centered, deviation-standardized trait. We accounted for measurement error for each trait by using maximum species standard deviation for each macroscopic trait from Weeks et al. (2020). We used maximum within-species standard error for each microscopic trait, estimated by measuring these traits on five flightless and five volant specimens of Fulica to assess the uncertainty around each trait.

Finally, we used a set of phylogenetic Bayesian regressions to evaluate changes in traits upon onset of flightlessness (Kruschke 2011; de Villemereuil et al. 2012; Bürkner 2017, 2018, 2021; Fuentes-G et al. 2020). These two-slope, two-rate phylogenetic mixture model regressions (MM) were fit using the R package RJAGS (Fuentes-G et al. 2020). These mixtures allowed lambda (phylogenetic signal) to vary, along with the evolutionary rates, intercepts, and slopes for flightless and volant species. We set priors following Fuentes-G et al. (2020), and all #Beta parameters had priors centered at 0, allowing us to test if the data support different regression structures (slopes & intercepts) as well as different rates for volant and flightless taxa. These rate estimates are separate from the rate estimates using OUwie described above; the MM rate estimates are on the direct trait values rather than standardized ratios and do not incorporate measurement error as our other analyses did. We performed these mixture-model (MM) regressions for the relationship between wing length as a function of body mass, as well as for the leading barb angles, lengths, and densities being a function of the trailing values. We compared estimates of the intercepts, slopes, and evolutionary rates between volant and flightless lineages using the highest density intervals of the Bayesian posteriors via the R package HDInterval (Meredith & Kruschke 2022).

Sampling error. Sampling error and inter-individual variation (often modeled phylogenetically as error) were examined to the extent possible given our dataset (Supplementary Data S5). We evaluated inter-individual variance in our estimates of phylogenetic rate across traits by using the greatest measured standard error in each feather trait from 10 specimens of Fulica, increasing uncertainty and support for the simpler one-rate Brownian motion model over the more complex two-rate model. Measurement error itself was examined by comparing the same measurements between observers, within and between feathers of a single individual, and between individuals of a single species.

Limited multi-individual sampling due to difficulty of measuring fine-scale feather traits from species, especially rare species, is common. Considering feather asymmetry, Feo et al. (2015) used a single specimen for each taxon, while Kiat & O’Connor (2024) measured three individuals, or less individuals for rarer taxa (e.g., flightless taxa).