PHYLOGENY AND FORELIMB DISPARITY IN WATERBIRDS

Xia Wang1 and Julia A. Clarke1

1

Department of Geological Sciences, University of Texas at Austin, Austin, TX 78712, Email: julia_clarke@ jsg.utexas.edu

Previous work has shown that the relative proportions of wing components (i.e., humerus, ulna, carpometacarpus) in birds are related to function and ecology, but these have rarely been investigated in a phylogenetic context. Waterbirds including “Pelecaniformes”, Ciconiiformes, Procellariiformes, Sphenisciformes and Gaviiformes form a highly supported clade and developed a great diversity of wing forms and foraging ecologies. In this study, forelimb disparity in the waterbird clade was assessed in a phylogenetic context. Phylogenetic signal was assessed via Pagel’s lambda, Blomberg’s K and permutation tests. We find that different waterbird clades are clearly separated based on forelimb component proportions, which are significantly correlated with phylogeny but not with flight style. Most of

This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process, which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1111/evo.12486.

This article is protected by copyright. All rights reserved.

1

the traditional contents of “Pelecaniformes” (e.g., pelicans, cormorants and boobies) cluster with Ciconiiformes (herons and storks) and occupy a reduced morphospace. These taxa are closely related phylogenetically but exhibit a wide range of ecologies and flight styles. Procellariiformes (e.g. petrels, albatross, and shearwaters) occupy a wide range of morphospace, characterized primarily by variation in the relative length of carpometacarpus and ulna. Gaviiformes (loons) surprisingly occupy a wing morphospace closest to diving petrels and penguins. Whether this result may reflect wing proportions plesiomorphic for the waterbird clade or a functional signal is unclear. A Bayesian approach detecting significant rate shifts across phylogeny recovered two such shifts. At the base of the two sister clades Sphenisciformes + Procellariiformes, a shift to an increase evolutionary rate of change is inferred for the ulna and carpometacarpus. Thus, changes in wing shape begin prior to the loss of flight in the wing-propelled diving clade. Several shifts to slower rate of change are recovered within stem penguins. KEY WORDS: forelimb evolution, disparity, phylogeny, waterbirds, wing-propelled diving.

This article is protected by copyright. All rights reserved.

2

Procellariiformes (petrels, albatross, shearwaters), Sphenisciformes (penguins), Gaviiformes (loons), Ciconiiformes (herons and storks) and the traditional contents of “Pelecaniformes” (frigate birds, pelicans, cormorants, gannets and boobies), together form a highly supported waterbird clade (Ericson et al. 2006; Livezey and Zusi 2007; Hackett et al. 2008; Morgan-Richards et al. 2008; Pratt et al. 2009) and exhibit significant diversity in flight style (e.g., flapping, flapping/gliding, flapping/soaring, dynamic soaring) and foraging ecology (e.g., feeding on the wing, from the water surface, pursuit underwater; Brewer and Hertel 2007; Hinić-Frlog and Motani 2009; Simons 2010). From soaring taxa with among the longest wingspans observed in living birds (e.g., Diomedea exulans; wingspan 3.5 m) to small bodied taxa (Halocyptena microsoma, wingspan 0.3 m) with highly manoeuvrable flight (del Hoyo et al. 1992), waterbirds present an ideal target clade for investigating patterns of forelimb modification. This clade also includes flighted and flightless taxa that utilize the wing in aquatic locomotion. Gannets and shearwaters use the wing in ascending or descending underwater after plunge-diving from the air, and penguins and diving-petrels habitually use their wings in pursuit of underwater prey (Adams and Walter 1993; Garthe et al. 2000; Burger 2001). In water, the wing must produce propulsive forces in a medium that is over 800 times denser and ~70 times more viscous than air (Denny 1993; Hamilton 2006). Sustained wing-propelled diving (WPD) has been hypothesized as a key locomotor transition in birds (e.g., Simpson 1946; Middleton and Gatesy 2000; Hamilton 2006; Thewissen et al. 2007; This article is protected by copyright. All rights reserved.

3

Clarke et al. 2010). Despite known variation, relatively few explicit hypotheses have been proposed concerning forelimb evolution in waterbirds or the tempo of forelimb proportional change associated with the origin of wing-propelled diving. Penguins were only noted as outliers in the assessment of forelimb disparity of Middleton and Gatesy (2000), and shortening of the wing has repeatedly been noted for flightless WPD taxa (e.g., Simpson 1946; Olson 1977; Livezey 1988, 1989; McCall et al. 1998; Middleton and Gatesy 2000). The avian forelimb is comprised of three main segments: humerus, ulna/radius, and carpometacarpus. The relative proportions of these forelimb segments both exhibit taxonomic variation and show a relationship with flight style (e.g., Padian and Chiappe 1998; Middleton and Gatesy 2000; Rayner et al. 2002; Nudds et al. 2004, 2007). Taxonomic variation in the relative proportions of the proximal wing (brachial index: the ratio between humeral length and ulnar length) of birds was noted by early anatomists (e.g., Beddard 1898; Steiner 1917; Böker 1927; Marples 1930; Verheyen 1961) and recovered more recently with larger datasets (Rayner et al. 2002; Nudds et al. 2004, 2007). This index has also been used as a discrete character in phylogenetic analyses of the relationships of both Mesozoic taxa (Chiappe 1995; Forster et al. 1998; Ji et al. 1998; Chiappe 2002; O’Connor and Zhou 2012) and extant clades of birds (Mayr and Clarke 2003; Cracraft et al. 2004). Middleton and Gatesy (2000) incorporated the length of the carpometacarpus (cmc) into quantitative analyses of forelimb disparity in theropod dinosaurs utilizing a proportion morphospace, an approach that has This article is protected by copyright. All rights reserved.

4

been utilized in an array of subsequent studies (e.g., Gatesy and Middleton 2007; Dyke and Nudds 2009; Nudds et al. 2013; Benson and Choinere 2013). In this study, forelimb disparity in the waterbird clade is assessed in a phylogenetic context to investigate evolution of the pectoral limb, potential relationships with flight style and to characterize potential shifts associated with evolution of wing-propelled diving. Our methods expand on those utilized in previous studies with dense taxonomic sampling including extant as well as extinct waterbird taxa (clade H; Hackett et al. 2008). Previous studies evaluated waterbird taxa as parts of broader analyses of variation in the forelimb in extant birds (e.g., Gatesy and Middleton 2000; Nudds et al. 2004, 2007), focused on one subclade (“Pelecaniformes”; Simons 2010) or assessed the appendicular morphometrics of aquatic taxa across an array of extant avian clades and outside the avian crown clade (Hinić-Frlog and Motani 2009). In all of these cases the primary aim of these studies was to investigate the relationship between proportions and function (e.g., flight style; aquatic foraging mode). While some of our questions and approaches are shared with these previous studies, initial assessment of waterbird forelimb proportions (Fig. 1) led to a focus on exploring phylogenetic trends in forelimb proportions in marked contrast to functional focus of these previous studies. Here, we assess metrics of phylogenetic signal, utilize ancestral trait reconstruction to test for rate shifts in forelimb proportions and describe patterns in the evolution of forelimb in the clade revealed in ternary space and phylomorphospace. This article is protected by copyright. All rights reserved.

5

Material and Methods DATASET AND VISUALIZATION Lengths of the three main forelimb segments (i.e., humerus, hu; ulna/radius, ul/ra; and carpometacarpus, cmc) for 335 specimens representing 146 species of extant and 7 extinct species were assembled from direct measurements of specimens and the literature (see electronic supplementary material, Table S1). The living bird taxonomic sampling includes all species of Gaviiformes, all genera of and approximately one half of extant species diversity within Sphenisciformes. Within Procellariiformes, all named families, more than one-half (15/26) of included genera were sampled. All genera and more than 70% (49/64) of extant species were sampled within traditional “Pelecaniformes”. Within Ciconiiformes (inc. storks, herons, hammerkop, shoebill and spoonbills), all named families, more than one-half (25/42) of genera were sampled. Data for all 334 waterbird specimens were used for the ternary plots (Fig. 1B). A larger data set comprising 447 extant birds (from Dyke and Nudds 2009) were used to show waterbirds relative to the total avian morphospace in ternary space (Fig. 1A). Mean values for 54 exemplar species were used for phylogenetic analyses (Figs. 1C, 2, 3). Ternary diagrams (Fig.1A, B) were first used to visualize limb proportion diversity in waterbirds. The length of the humerus, radius (or ulna), and carpometacarpus were summed to obtain total forelimb length. Each of the three segments was then divided by total forelimb This article is protected by copyright. All rights reserved.

6

length to calculate its percentage of the total. The percentage of each segment for each specimen was plotted on a ternary diagram (see Middleton and Gatesy 1997, 2000 for a more detailed explanation of the application of these diagrams). To quantitatively compare the forelimb morphospace occupied by different clades, the ternary diagram was subdivided into 400 triangular cells and the number of occupied triangles (cells) was counted, following the methods of Gatesy and Middleton (1997). The amount of variation in each forelimb segment was also summarized by comparing the standard deviation (SD) associated with mean values. ANOVA, including a post-hoc test and a t-test, was also used to evaluate if the standard deviation values were significantly different among wing segments for the major waterbird subclades considered here (e.g., Procellariformes, ‘Pelecaniformes’+Ciconiformes, Sphenisciformes).

PHYLOGENETIC TREE Because no single analysis included all of our target waterbird taxa, we utilized a supertree with the phylogeny of Hackett et al. (2008) as the backbone. The relationships among waterbirds recovered by Hackett et al. (2008) are the same as those in more recent published phylogenies (Jetz et al. 2012; Yuri et al. 2013) for taxa in common among these analyses. Turacos (Musophagiformes) were included as sister taxon to the waterbird taxa as recovered in these studies (Hackett et al. 2008; Jetz et al. 2012; Yuri et al. 2013). We used This article is protected by copyright. All rights reserved.

7

the following references for the relationships within families represented by more than two species in the analysis: Sphenisciformes (Baker et al. 2006; Ksepka and Clarke 2010; Clarke et al. 2010), Phalacrocoracidae (Kennedy et al. 2009, 2013; Holland et al. 2010); Procellariidae (Gangloff et al. 2012); Sulidae (Patterson et al. 2011); Fregatidae (Smith 2010). Because ancestral state reconstruction and inferred patterns of character evolution are sensitive to the branch lengths used (Hutcheon and Garland 2004; Garland et al. 2005; Hinić-Frlog and Motani 2009; Kilbourne and Makovicky 2010; Laurin 2010), the ideal situation would be to know the true branch lengths in terms of absolute time (or in number of generations) since common ancestry (Blomberg et al.2003; Garland et al. 2005). However, divergence time estimation for higher level avian phylogeny remains controversial (e.g., Ericson et al. 2006; van Tuinen et al. 2006; Brown et al. 2008; Ksepka et al. 2011), and no comprehensive analysis exists for the taxonomic sample we consider. In such cases, many comparative studies employ “arbitrary” branch lengths by either setting all branches to unit length, or aligning contemporary species and varying the depths of internal nodes (e.g., Grafen 1989; Pagel 1992; S. Nee cited in Purvis 1995). Here, we utilize molecular branch lengths where available and scaled wholly extinct branches. We used the exact branch lengths from an ultrametric time tree of the Fig. 3 topology of Hackett et al. (2008) for all available nodes. This tree was used in Braun et al. (2011) and This article is protected by copyright. All rights reserved.

8

Han et al. (2011) in which non-parametric rate smoothing was employed. For all other branches, we first set all branch lengths to unit 1, then applied a semi-parametric penalized likelihood approach (Sanderson 2002) to render these branches ultrametric by using function chronos in R package Ape (Paradis et al. 2013). The relative lengths of the recovered branch lengths were similar in relative lengths to those in ultrametric tree from Hackett et al. (2008) for comparable internodes but with a max tip to base length= 1. We then multiplied these branch lengths by 93.54 to match the max tree tip to base length (93.54 Ma) used in ultrametric transformation of the Hackett et al. (2008) tree (Braun et al., 2011). For those nodes with fossils, we further subtracted the fossil ages from the total branch length. This transformation (hereafter “primary”) while not ideal, approximates what is known about branch lengths in the clade better than unit branch lengths or simple ultrametric tranformations of arbitrary sets of branch lengths.

ASSESSMENT OF PHYLOGENETIC SIGNAL We compared three approaches for two different sets of branch lengths (primary tree and unit=1) to investigate the degree to which patterns of phylogenetic signal in the forelimb proportion data might be sensitive to the statistic methods and branch length used. First we performed a permutation test (Klingenberg and Gidaszewski 2010) in Mesquite. In this test, the phylogeny is held constant and the wing proportion data for each taxon are repeatedly This article is protected by copyright. All rights reserved.

9

and randomly swapped across the tree 100,000 times. The null hypothesis of no phylogenetic signal is rejected if less than 5% of permutations result in a tree length (calculated using squared-change parsimony) that is shorter than or equal to the value obtained from the original data (Klingenberg and Gidaszewski 2010). Second, we tested Pagel’s lambda (Pagel 1999) using the R package GEIGER (Harmon et al. 2008). The influence of the phylogeny increases with lambda from 0 (no phylogenetic signal) and 1 (strong phylogenetic signal). When lambda = 0, a star phylogeny results with all tips radiating from a basal node, describing a model where traits evolve independent of the phylogeny. When lambda = 1, trait evolution followed a Brownian motion model, where branch lengths would be proportional to divergence. To determine whether lambda is significantly different from zero, we used a likelihood ratio test in R (Harmon et al. 2008). Finally, we estimated Blomberg’s K (Blomberg et al. 2003) using the R package PHYTOOLS (function phylosig, Revell 2012). Blomberg's K statistic (Blomberg et al. 2003) is used to test whether the observed distribution of traits exhibits more or less divergence than expected for traits evolving under Brownian motion. Values of K close to 1 indicate trait similarity is proportional to divergence and a Brownian motion model of evolution fits the data. K > 1 indicates that close relatives are more similar than expected, and K < 1 indicates more divergence between taxa than expected under a Brownian model. The metrics of Pagel’s lambda and Blomberg’s K are quite different from one another. λ is a scaling parameter for the covariances between species, relative to the This article is protected by copyright. All rights reserved.

10

covariances expected under Brownian evolution (Pagel 1999). K is a scaled ratio of the variance among species over the contrasts variance and can inform us about trait variation that is more similar than expected under Brownian motion (Blomberg et al. 2003). These metrics vary in sensitivity to branch length and taxon sampling, with Pagel's lambda least sensitive to taxon number and Blomberg's K more sensitive to the absence of branch length information. The upper limit of Pagel’s lambda is close to one, while Blomberg’s K can yield higher values indicating much stronger phylogenetic signal (Blomberg et al. 2003; Cavender-Bares et al. 2006; Mϋnkemϋller et al.2012).

ANCESTRAL STATE RECONSTRUCTION AND ESTIMATION OF EVOLUTIONARY RATE SHIFTS We used Mesquite (version 2.75; Maddison and Maddison 2011) to map the relative forelimb proportions onto the reference phylogenetic tree. Each character was traced onto the tree using the “reconstruct ancestral state” module of Mesquite with weighted squared-change parsimony (Maddison 1991). Given the tree and observed character distribution, this method finds the ancestral states that minimize the number of steps of character change. A 3-D phylomorphospace (Sidlauskas 2008) plot was also obtained from Mesquite (version 2.75, build 564; Maddison and Maddison 2011; Fig. 1C).

This article is protected by copyright. All rights reserved.

11

We also investigated shifts in rate of evolution for each forelimb segment using both Pagel’s delta (Pagel 1999) and a Bayesian approach. Pagel's delta transformation is a model of increasing or decreasing rates of trait change through time (Pagel 1999). Values of delta < 1 transform branch lengths to be increasingly shorter towards the tips, describing a model in which trait change occurs rapidly early in the history of a clade and then slows through time. Models with delta > 1 describe an increasing rate of character evolution through time, while delta = 0 is identical to a Brownian motion model. We tested Pagel’s delta for two different sets of branch lengths, those from the primary tree and unit=1. The Bayesian approach, implemented in the R package AUTEUR version 0.12.10 (Eastman et al. 2011) uses reversible-jump Markov Chain Monte Carlo (MCMC) sampling to draw from a distribution of multi-rate models in inverse proportion to their poorness of fit. Available models range from one in which every branch segment of the tree has its own rate parameter, to one in which a single rate characterizes the whole tree. This process results in model-averaged posterior rate distributions and estimates of the posterior probability of rate shifts occurring at particular nodes on the tree (Eastman et al. 2011). The set of evolutionary rates scalars bears on the stochastic variance of trait values expected to accumulate within the tree. We implemented this approach over 1,000,000 generations ten times for each parameter utilizing the branch lengths of the primary tree. The first 250,000 generations were discarded as burn in. The results of Bayesian rate estimation on relative forelimb length are This article is protected by copyright. All rights reserved.

12

shown in Figs. 2 and 3. The significance of observed rate shifts was assessed using randomization tests (Eastman et al. 2011).

PHYLOGENETIC ANOVA In addition to testing for phylogenetic signal in the forelimb proportion data, we were also interested to see if forelimb proportions are influenced by flight styles. Phylogenetic ANOVA including a post-hoc test (using Holm–Bonferroni method) on means was conducted in R using the package PHYTOOLS (Revell 2012) to see if forelimb proportions significantly different between any two flight styles. Four flight styles for living birds were used: ‘continuous flapping’ (CF) (e.g., loons, herons, cormorants and diving petrels); ‘flapping/soaring’ (FS) (e.g., pelicans, storks, frigatebirds, gannets and anhingas); ‘dynamic soaring’ (DS) (e.g., albatrosses, shearwaters); ‘flapping/ gliding’ (FG) (e.g., boobies, storm petrels) according to Bruderer et al. (2010) and Simons (2010). Flightless WPD species were coded as a distinct class. Results Distinct waterbird subclades are clearly separated based on wing component proportions (Fig. 1B). Most of the traditional contents of “Pelecaniformes” (e.g., pelicans, cormorants and boobies) cluster with Ciconiiformes (herons and storks) and occupy a limited morphospace (3 of 400 cells=0.75%), exhibiting significantly longer ul% (t=18.79, P

Phylogeny and forelimb disparity in waterbirds.

Previous work has shown that the relative proportions of wing components (i.e., humerus, ulna, carpometacarpus) in birds are related to function and e...
984KB Sizes 2 Downloads 4 Views