The evolutionary origin of avian facial bristles and the likely role of rictal bristles in feeding ecology

Samples

We examined 1,022 avian species (~ 10% recorded species) in this study, representing 418 genera, from 91 families (37% recorded families) and 29 orders (73% of all orders). Specimens were from the skin collection of the World Museum Liverpool, Tring Natural History Museum, Manchester Museum and Wollaton Hall Museum, all situated in the United Kingdom. All work was carried out in accordance with ethical regulations at Manchester Metropolitan University and with the permission of all aforementioned museums. Only the best-preserved adult specimens (no signs of cut off feathers or holes in the skin near the beak) were chosen for this study to ensure accurate measurements of bristle length, shape and presence, which should not be affected by the process of skin removal and specimen conservation. Species were randomly chosen, without targeting our sampling towards species known a priori to have bristles. Where possible, two specimens per species were measured (occurring in 82% of all species examined). Specimens of each sex were measured when present; however, this was not always possible since labelling was often inaccurate or missing. In total, the sample included 508 males, 412 females and 374 individuals of unknown sex. Both sexes were examined in 274 species and there was no difference whatsoever between the presence of bristles on male or female species (n = 97 with bristles present and n = 180 with bristles absent for both males and females). Length (Mann–Whitney U test, W = 37,962, N = 552, P = 0.94) and shape (Chi-square test, χ² = 0, N = 552, df = 3, P = 1) of rictal bristles also did not significantly differ between males and females. Therefore, rictal bristles are likely to be sexually monomorphic and data for males and females was pooled for further analyses. Overall, rictal bristles were absent in 64% of species examined (n = 656) and just over a third of species (n = 366) had bristles present.

Bristle descriptions

Facial bristles were initially identified by sight and touch in each specimen. Bristles were recorded as either present or absent from the upper rictal, lorial, lower rictal, narial and interramal regions (Fig. 1a). We use the term ‘rictal bristle’ here for bristles on both the upper rictal and/or the lorial region, since there was no clear differentiation and morphological differences between the bristles found in these regions forming a continuum of bristles above the edge of the beak. When present, rictal bristle shape was recorded as: (i) unbranched rictal bristles, (ii) rictal bristles with barbs only at the base (“Base”) and (iii) branched rictal bristles (“Branched”), i.e. barbs and barbules present along the bristle rachis (Fig. 1b). The three longest rictal bristles were measured on both sides of the head of each specimen using digital callipers, and these lengths were averaged to provide a mean length of rictal bristles per species. In species lacking rictal bristles, a length of “0” and a shape category of “Absent” was recorded.

Ancestral reconstruction of facial bristle presence

Following Felice et al.¹⁹, a single consensus phylogenetic tree was generated from the Hackett posterior distribution of trees from Birdtree.org²⁰ with a sample size of 10,000 post burn-in, using the TreeAnnotator utility in BEAST software²¹ with a burn-in of 0. Maximum Clade Credibility (MCC) with the option “-heights ca” was selected as the method of reconstruction. The common ancestor trees option (-heights ca) builds a consensus tree by summarising clade ages across all posterior trees. Both the consensus tree and posterior distribution of 10,000 trees were imported into RStudio v. 1.2.5 for R^22,23 and pruned so that only species present in the dataset of this study remained in the phylogeny. Taxon names were modified where necessary to match those from the Birdtree.org (http://birdtree.org) species record. Negative terminal branches in our consensus tree were slightly lengthened to be positive using ‘edge.length[tree$edge.length < 0] = 1e−6’. Then, we coerced the tree to be ultrametric using the method ‘nnls.tree’ (phangorn package)²⁴ in the function force.ultrametric (package phytools v0.7–70)²⁵. The root age of the basal node of the consensus tree was calculated using the function tree.max from the FossilSim package²⁶.

The ancestral state reconstruction (ASR) of the presence/absence character was conducted using the function make.simmap (package phytools v0.7-70)²⁵, which simulated stochastic character mapping by using the binary character presence/absence of rictal bristles on the consensus tree (nsim = 10,000), with the results summarised by using the function describe.simmap (phytools). Three commonly-used evolutionary transition rate models—equal-rates (ER), symmetrical (SYM) and all-rates-different (ARD)—were evaluated across the posterior distribution of 10,000 trees using the ace function in ape v5.4-1 package²⁷. For the ER model, all transitions rates were governed by a single parameter; for the SYM model, transitions 0 → 1 and 1 → 0 occur at the same rate but may differ from the transition rate between state 1 and 2; and in ARD model, each rate referred to a different parameter. Model fits were evaluated using the fitDiscrete function in the R package geiger v2.07²⁸. Model fits were determined using the AIC values (Akaike Information Criterion) and AIC weight (AIC.w)^29,30,31,32. The comparison between different transition rate models revealed that the equal rate (ER) model was rejected in favour of a more parametrised all rates different (ARD) model for rictal, lower rictal and interramal bristle presence while the ER model was selected for the narial bristle presence (Table 1). The densitymap function (phytools) was used to plot the consensus tree on to which the posterior density of rictal bristle presence/absence was mapped. The mapped value represented the probability of having rictal bristles present. Adobe Photoshop CC 2018 software (Adobe Systems Incorporated, San Jose, California) was used to customise the resulting radial tree.

Sensitivity analysis

To check the robustness of our analyses to the trait used (bristle presence) and examine the uncertainty of our ancestral reconstruction for the presence of each facial bristle, we calculated the phylogenetic signals of each trait and evaluated both taxonomic sampling bias and tree.

Phylogenetic signals

We examined the phylogenetic signals to determine whether the distribution of each trait across our phylogeny follows a Brownian motion evolution hypothesis or shows divergent evolutionary trends suggesting potential selective pressures acting on our trait. We used Pagel’s λ³³ for rictal bristle presence character (binary traits), using the phylosig function (package phytools)²⁵, to evaluate its phylogenetic signal; a lambda value (λ) close to 1 and P < 0.05 signifies that the trait evolution would be predicted by Brownian motion model and a lambda value (λ) of 0 would signifies that the character has no phylogenetic signal.

Sampling bias

The robustness of the ancestral state reconstruction findings for facial bristle presence was evaluated by exploring the effects of sampling bias of the taxa and trees within the datasets for each make.simmap analysis. To evaluate taxonomic sampling bias, we calculated a probability of oversampled or undersampled families (number of representative of a family in the dataset divided by the total number of species known to exist in the family) to weight the likelihood of a given bird species appearing in a downsampled dataset (e.g.³⁴). The probability was subsequently used in the weighted analysis by downsampling the datasets (full dataset, and presence only subset) to 70%, 80% and 90% of the species, preferentially removing families that are overrepresented. The stochastic mapping procedure was then repeated 10,000 times, with the mean number of transitions, gains and losses calculated, as well as the average time the character spent in each state, for each of the subsamples (70, 80, 90%) (Table 2).

Tree topology uncertainty

We explored the uncertainty of tree topology and branch length in the ancestral state reconstruction analysis by randomly sampling 100 trees from the posterior distribution, in addition to the consensus tree. For each of these 100 trees, 100 character histories were randomly sampled, generating 10,000 character maps to account for the uncertainty associated with tree topology, branch length and timing of the transitions between morphological states (Table 2). We also tested if alternative phylogenetic topologies for Palaeognathae and Caprimulgimorphae, which would match recent studies (e.g.^13,35), would have an effect on the ancestral state reconstruction analysis for facial bristle presence (Table S2), and obtained similar results to the analysis with our original consensus tree (Table S1, S2). 

Ecological traits

Species-specific ecological traits were added to the dataset using birdsoftheworld.org¹⁴, and included the following trait variables: (i) period of activity, (ii) habitat type, (iii) foraging method, (iv) foraging height, and (v) diet (Table 3). Instances where some species belonged to more than one dietary guild, diet categories were based on a combination of maximum two dietary guilds—the first guild (a; Table 3) corresponding to the main food type and the second guild (b, c, d; Table 3) corresponding to the secondary food type of a species (e.g. a–b, a–c, a–d, etc.). A total of 22 combinations were generated in our dataset for diet, e.g. Invertivore–Granivore, Invertivore–Vertivore and Invertivore–Frugivore, or Vertivore–Invertivore, Frugivore–Granivore or Frugivore–Herbivore, etc. Similarly, in instances where species exhibited more than one foraging method, the foraging method category was based on a combination of the two main foraging method. A total of 4 combinations were generated for the foraging method category, i.e. Gleaning–Hawking or Gleaning–Sallying, Hawking–Sallying, and Sallying–Plunge-diving.

Model construction for bristle association with ecological traits

The same dataset of 1022 avian species and the consensus phylogenetic tree was used for this analysis. However, for this investigation, the rictal bristle length was calculated per individual (the average length of the rictal bristles measured per individual), rather than using a mean value per species.

To determine the relationship of both rictal bristle presence and length with species-specific ecological traits, a phylogenetically controlled Markov chain Monte Carlo generalised linear mixed model (MCMCglmm) was conducted, using the R package “MCMCglmm”³⁶ in RStudio²³. For rictal bristle presence, a binomial ‘threshold’ model was used to account for the binary response variable, whereas a “Gaussian” model was used for the rictal bristle length continuous response variable. In both models, the period of activity, the habitat type, the foraging method, the foraging height and diet were included as fixed effects. In both models, phylogeny and individual ID were included as independent random effects.

Following Hadfield³⁶, a weak informative inverse-Gamma prior was used in the models, with variance (V) set to 1, and the belief parameter (nu) set to 0.002 for both the random effects structure (G-structure) and residual structure (R-structure). Residual variance was fixed in the absence of this information for the rictal bristle presence model since this used binary data (as per Hadfield, 2010³⁶). Other parameter combinations were systematically explored but the models did not converge with them. The model was run for 800,000 iterations, with a burn-in period of 80,000, and a thinning of 40, which were determined using diagnostics in the coda package³⁷. Three independent MCMC chains were run per model to check for model convergence using Gelman-Rubin diagnostics, with model convergence confirmed when the potential scale reduction factor required value was < 1.1³⁸. Effective sample sizes (> 200) and autocorrelation (P < 0.05) values between successive iterations were also examined. Non-significant fixed effects (pMCMC > 0.10) were permanently excluded from the model formula if, in doing so, the fit of the model improved. Using ggplot2 package³⁹, we constructed caterpillar plots representing the mean parameter estimates and the 95% credible intervals (CI) for each model. If the credible intervals were found to exclude zero, the parameter was considered significant with the model P-value given by the pMCMC value.

It was not possible to obtain a converging model for bristle shape since this was a categorical variable that exhibited a large range in the number of species in each category i.e. unbranched bristle shape was found in 292 species, while the branched shape was only found in 47 species, and branched at the base found to be present in 26 species. While visual inspection of MCMC chains suggested convergence after 12.8 million iterations, convergence was not supported by the Gelman–Rubin statistic; thus, a model for bristle shape was not considered further.

Following model construction and validation, a suitable ‘reference category’ was selected for pairwise comparisons. These reference categories are compared to all others categories within their ecological traits, and tested for significant differences⁴⁰. Since all Anatidae species recorded in our dataset did not have rictal bristles and shared the same ecological categories, these categories were selected as reference categories for the rictal bristle presence and length models. Therefore, the reference categories were diurnal, open, dabbling, low and aquatic herbivore.

Source: Ecology - nature.com