DOI: 10.5061/dryad.qrfj6q5s7
title: Feather evolution following flight loss in crown group birds: relaxed selection and developmental constraints

description: Feathers are complex structures exhibiting structural/functional disparity
 across species and plumage. Flight was lost in &gt;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.
description: <em>Taxon Sampling.</em> 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 (<em>Spheniscus</em>,
 <em>Eudyptes</em>, <em>Aptenodytes</em> [Ksepka
 <em>et al.</em> 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 <em>et al</em>. 2017): 1)
 <em>Struthio camelus</em>, 2) <em>Rhea americana, 3)
 Dromaius novaehollandiae</em>/<em>Casuarius
 unappendiculatus</em>, and 4) <em>Apteryx
 australis</em>. Six of our flightless taxa are
 recently extinct (<em>Porzana sandwichensis</em>,
 <em>Porzana palmeri</em>, <em>Gallinula
 nesiotis</em>, <em>Xenicus lyalli</em>,
 <em>Podilymbus gigas</em>, <em>Pinguinus
 impennis</em>). One recently extinct species was
 poor-flying/‘incipiently flightless’ (<em>Mergus
 australis</em>). 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
 (<em>Oceanites oceanicus</em>, <em>Pelecanoides
 urinatrix</em>, <em>Diomedea
 immutabilis</em>). We expanded sampling of crown
 birds to include at least one species from each commonly recognized order
 (Reddy <em>et al.</em> 2017) All these broadly sampled species
 are volant, except for <em>Mesitornis unicolor</em>, whose
 flight capabilities are unknown. <em>Data
 Acquisition.</em> 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 <em>et al</em>.
 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. <em>Equipment.</em> 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). <em>Feathers sampled.</em> 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. <em>Observations of each feather.</em> 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 <em>et al.</em> 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 <em>et al</em>. 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 <em>et al</em>. 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 <em>et al</em>. 2003; Prum 2005; Lin <em>et al</em>. 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.         <em>Analysis.</em> 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 <em>et al.</em> 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: Δ = Metric<sub>Flightless</sub> – Metric<sub>Volant</sub>. For penguins, volant values were averaged from three tubenoses (<em>Oceanites oceanicus</em>, <em>Pelecanoides urinatrix</em>, <em>Diomedea immutabilis</em>). 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. <em>Phylogenetic comparative methods (PCM).</em> A series of phylogenetically controlled comparisons tested how rapidly traits changed after flight loss (Supplementary Data S4). Using the phylogeny from Jetz <em>et al</em>. (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 <em>et al</em>. 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 <em>et al</em>. (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 <em>Fulica </em>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 <em>et al.</em> 2012; Bürkner 2017, 2018, 2021; Fuentes-G <em>et al</em>. 2020). These two-slope, two-rate phylogenetic mixture model regressions (MM) were fit using the R package RJAGS (Fuentes-G <em>et al</em>. 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). <em>Sampling error</em>. 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 <em>Fulica</em>, 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 <em>et al</em>. (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).
description: # Feather evolution following flight loss in crown group birds: relaxed
 selection and developmental constraints
 [https://doi.org/10.5061/dryad.qrfj6q5s7](https://doi.org/10.5061/dryad.qrfj6q5s7) ## Description of the data and file structure Supplemental Material Sections 1-10.pdf and Supplemental Data S1-S5.zip can be found here. They relate to the publication "FEATHER EVOLUTION FOLLOWING FLIGHT LOSS IN CROWN GROUP BIRDS: RELAXED SELECTION AND DEVELOPMENTAL CONSTRAINTS". ### Files and variables #### File: Supplementary\_Data\_S1-S5.zip **Description:**  S1_Dataset folder: Dataset.csv Dataset.xlsx The main dataset containing morphological and ecological data from flightless and volant birds. All variables are fully written out in the first row. Can be opened with software such as Microsoft Excel. Empty cells are NA. S2_Main R code folder: Dorsal contours subfolder: Dorsal Contour_Asymmetry.R Dorsal Contour_Density.R Dorsal Contour_Middle_Long branches_Lengths.R Dorsal Contour_Middle_No Long branches_Lengths.R R scripts used to analyze evolutionary shifts in the morphology of dorsal contour feathers. Can be opened with software such as RStudio. See S1_Dataset for variable descriptions. Primaries subfolder: Primary_Asymmetry.R Primary_Density_SCALED ATTEMPT.R Primary_Density.R Primary_Middle_Long branches_Lengths.R Primary_Middle_No Long branches_Lengths.R R scripts used to analyze evolutionary shifts in the morphology of primary feathers. Can be opened with software such as RStudio. See S1_Dataset for variable descriptions. Rectrices subfolder: Rectrix_Asymmetry.R Rectrix_Density.R Rectrix_Middle_Long branches_Lengths.R Rectrix_Middle_No Long branches_Lengths.R R scripts used to analyze evolutionary shifts in the morphology of retrix feathers. Can be opened with software such as RStudio. See S1_Dataset for variable descriptions. Tertials subfolder: Tertial_Asymmetry.R Tertial_Density.R Tertial_Middle_Long branches_Lengths.R Tertial_Middle_No Long branches_Lengths.R R scripts used to analyze evolutionary shifts in the morphology of tertial feathers. Can be opened with software such as RStudio. See S1_Dataset for variable descriptions. Ventral contours subfolder: Ventral Contour_Asymmetry.R Ventral Contour_Density.R Ventral Contour_Middle_Long branches_Lengths.R Ventral Contour_Middle_No Long branches_Lengths.R R scripts used to analyze evolutionary shifts in the morphology of ventral contour feathers. Can be opened with software such as RStudio. See S1_Dataset for variable descriptions. Whole plumage subfolder: ALL_Middle_Long branches_Lengths.R ALL_Middle_No Long branches_Lengths.R R scripts used to analyze the morphology of all feathers in the dataset. Can be opened with software such as RStudio. See S1_Dataset for variable descriptions. S3_Asymmetry or Rachis Legnth vs Mass folder: Asymmetry vs Mass code.R R script used to analyze the correlation between either feather asymmetry or rachis length versus body mass. Can be opened with software such as RStudio. See S1_Dataset for variable descriptions. S4_PCM Github files folder: see Github ([https://github.com/paleomitchelljs/FeatherEvo](https://github.com/paleomitchelljs/FeatherEvo)). See S1_Dataset for variable descriptions. data subfolder: fulica_intraspecific.csv primary.csv sisterpairs.csv Weeks_etal_2020_Data.csv These data files are the main inputs to the code contained in the /scripts/ folder. Empty cells in all files indicate inapplicable data (values of NA). primary.csv contains the same values as in S1_Dataset, although they have been slightly renamed. "swim" is the swimming score value from S1_Dataset "body_len" is the body length (beak tip to end of tail in mm) value from S1_Dataset "mass" is the Best single body value (g) from S1_Dataset "wing" is the tip of longest primary to wrist (mmm) from S1_Dataset "tail" is the tail fan length (mm) from S1_Dataset "tarsus" is the tarsus length (mm) from S1_Dataset "barb_angle_m_l" is the Primary: Angle barb LM from S1_Dataset "barb_angle_m_t" is the Primary: Angle barb TM from S1_Dataset "rachis_width_m" is the Primary: middle rachis width (mm) from S1_Dataset "rachis_length_m" is the Primary: middle rachis length (mm) from S1_Dataset "barb_length_m_l" is the Primary: Length barb LM (mm) from S1_Dataset "barb_length_m_t" is the Primary: Length barb TM (mm) from S1_Dataset "barb_density_m_l" is the Primary: Middle leading barb density (#/mm) from S1_Dataset "barb_density_m_t" is the Primary: Middle trailing barb density (#/mm) from S1_Dataset "barbule_length_m_t" is the Primary: Distal barbule length on barb TM (mm) from S1_Dataset "barbule_density_m_distal" is the Primary: Distal barbule density on barb TM (#/mm) from S1_Dataset "barbule_density_m_proximal" is the Primary: Proximal barbule density on barb TM (#/mm) from S1_Dataset fulica_intraspecific.csv contains the same data, although for multiple specimens of species within the genus Fulica. sisterpairs.csv contains the sister-taxon pairing guide for species within our dataset. The first column (Taxon) contains the species name, the second column (flighted) contains whether the species is volant (1) or flightless (0), the third column (clade) contains the name of the larger clade to which the species belongs, the fourth column (max_age) contains the maximum divergence time (if known) of the species from its sister-taxon in this dataset, and the fifth column (pairs) contains a code denoting which sister-pair a given taxon belongs to. Weeks_etal_2020_Data.csv contains the Weeks et al. (2020) dataset used in the intraspecific comparison analyses. Data used for phylogenetic comparative analyses. Can be opened in software such as Microsoft Excel. See S1_Dataset for variable descriptions. Empty cells are NA. phy subfolder: crane_tree_renamed_alda.tre crane_tree_renamed.tre fulltree.tre garcia_gruiformes.newick prum_etal_2015_vegavis.tre PrumSpecies_hand.csv simtrees.tre wholetree.tre Phylogenetic trees used for phylogenetic comparative analyses. Can be opened in software such as Mesquite or R. PrumSpecies_hand.csv is a reference file converting the tip labels in the Prum et al (2015) trees to species names that match those in our dataset and contain the location (column 3, Bind) in the tree where an extinct or subspecific tip should be placed using bind.tip() in R. fulltree.tre is the primary phylogeny file used in the analyses described in the paper. The other files were used to construct this file, and are described below. The script xx00_trees.R uses garcia_gruiformes.newick to generate crane_tree_renamed.tre and crane_tree_renamed_alda.tre. The script xx00_reconcileTrees.R uses prum_etal_2015_vegavis.tre, PrumSpecies_hand.csv, and crane_tree_renamed_alda.tre to generate wholetree.tre The script xx04_simmaps.R generates the stochastic maps for flightlessness and saves them as simtrees.tre scripts subfolder: multi.txt spindlePlot.R spreadlabels.R xx00_readFormatDat.R xx00_reconcileTrees.R xx00_trees.R xx01_sisterPairs.R xx02_logisticRegression.R xx03_phyloLogReg.R xx04_simmaps.R xx05_PCMs.R xx06_table2_pcmRatePlots.R xx07_pPCA.R xx08_LogRegPlots.R xx09_rateTest.R xx10_hypTests.R xx11_fig5.R xx12_tableBayesReg.R R scripts used for phylogenetic comparative analyses. Can be opened in software such as RStudio. All scripts contain setwd() calls that reference other subfolders in the S4 folder, such as /data/ or /phy/, and these working directories will need to be changed to match the local computer to run them. spindePlot.R and spreadlabels.R are generic functions required to make figures look nicer, and are used by xx11_fig5.R to generate Figure 5 in the main text. xx00_reconcileTrees.R and xx00_trees.R help generate the tree file, and their inputs and outputs are described in the section above. xx00_readFormatDat.R reformats the data files from the /data/ subfolder into the format needed for further analysis (e.g., mean-standardizing and log-transforming where needed). xx01_sisterPairs.R uses the sisterpairs.csv file as input to produce subsets of the main datafile that correspond to the traits for flightless and volant species. xx04_simmaps.R generates the Stochastic Mappings used by subsequent analyses. xx09_rateTest.R contains the script to run the rJAGS model of Fuentes-G et al. (2019) for all of our data across all of our trees. xx11_fig5.R contains the script to generate Figure 5 in the main text. xx12_tableBayesReg.R contains the script to generate Table 2 in the main text, containing the resulting estimates from fitting the Fuentes-G et al (2019) regression model to our data. S5_Sampling error folder: Strix comparison.xlsx Measurement error between two researchers.xlsx Fulica intraspecific comparison.xlsx Data used to better judge sampling/measurement error in the dataset. Can be opened in software such as Microsoft Excel. See S1_Dataset for variable descriptions. Feo et al 2015 subfolder: Feo Code.R Feo-2015_PRSB_PaperData.csv Feo-2015_PRSB_PaperData.xlsx New Feo Dataframe.csv Data from Feo et al. (2015) which we also use to better judge sampling/measurement error in the dataset along with the R script used for this analysis. Can be opened in software such as Microsoft Excel or RStudio. See S1_Dataset for variable descriptions. #### File: Supplementary Material Sections 1-10.pdf A detailed appendix. Can be opened in software such as Adobe Acrobat. ## Sharing/Access information Links to other publicly accessible locations of the data: * [https://github.com/paleomitchelljs/FeatherEvo](https://github.com/paleomitchelljs/FeatherEvo)