Study system and phylogenetic inference
The Psittaculidae comprise 175-180 species distributed across Africa and Asia through to Australasia and Oceania, and vary moderately in body mass, from about 15-30 g (pygmy parrots Micropsitta spp., budgerigar Melopsittacus undulatus and the blue lorikeet Vini peruviana) to over 500 g in the eclectus parrot Eclectus roratus. We restricted the present study to the Australasian region (Australia, New Zealand, New Guinea, and surrounding islands) which harbors the greatest proportion of taxonomic diversity (39 of 45 genera and 68% of species of Psittaculidae inhabit this region) and has been suggested as ancestral area of origin of all extant parrots [30].
We inferred the phylogenetic relationships among species from the consensus tree assembled by [57]. It consists of a calibrated Matrix Representation with Parsimony (MRP) phylogenetic supertree including all 398 extant species of parrots (order Psittaciformes). Since parrots have a poor fossil record [58], the Burgio supertree relies on well-accepted external fossil calibrations points, which are described in detail in [30-31, 41]. To account for phylogenetic uncertainty in some analyses, we also used a sample of 500 trees obtained from BirdTree (www.birdtree.org) based on the ‘stage 2 parrot’ backbone phylogeny. BirdTree combines relaxed clock molecular trees of well-supported avian clades with a fossil calibrated backbone with representatives from each clade [59].
Morphological and ecological data
We first characterized the overall morphology of each of the 117 species from average measurements on body length (size), tarsus length, tail length, culmen length, and the Hand Wing Index (HWI). These measurements were obtained from the Handbook of Australian, New Zealand and Antarctic birds [60] and Parrots of the World [53], except for HWI, which were retrieved from [61]. These traits are inherently linked to ecological and behavioral characteristics such as foraging niche and diet preferences. For instance, tarsus length has strong influence on the foraging mode and provisioning behaviour, whereas the length of the tail affects maneuverability and stability and thus, confers aerodynamic properties [62]. The HWI describes the wing aspect ratio and, consequently, has been frequently used as a proxy of dispersal ability and flight efficiency (e.g., [63]). Specifically, HWI represents the distance between the tip of the first secondary feather and the tip of the longest primary feather (Kipp’s distance) corrected for wing size [61]. The first four variables (body size, culmen length, tarsus length and tail length) were strongly correlated with each other after log-transformation. To avoid this and remove the confounding effects of body size while accounting for phylogenetic relationships among species, we performed phylogenetic regressions in ‘phytools’ using the BM method to obtain the correlation structure [64]. The obtained residuals for the three remaining traits (i.e., size-corrected variables) were used as input in a principal component analysis (PCA) in order to reduce the dimensionality of our trait dataset. The PCA yielded a phylogenetic principal component axis (PCa1) that accounted for 70% of the total variance. Tarsus and culmen length were positively loaded (0.52 and 0.61, respectively), whereas tail length showed the opposite trend (factor loading: -0.59).
We then focused in more detail on the beak morphology. We took advantage of the global dataset of morphological measurements compiled by J. Tobias and collaborators from live-caught individuals and preserved museum skins (e.g., [65]). This dataset includes five traits describing variation in beak shape [(1) beak length measured from fore-edge of cere to tip along the culmen; (2) beak length measured from the tip to the anterior edge of the nares; (3) beak depth measured at the anterior nares; (4) beak width measured at the anterior nares; and (5) gape width measured at the join of the mandibles] from a sample of adult individuals from each species (average number of measured individuals per species: 5.4). Since sexual size dimorphism is absent or minor in this taxonomic group [53], measurements on males and females were pooled together when computing the average values for each species. From these five linear measurements we conducted another PCA in a similar way to that previously indicated once these were log-transformed and size-corrected. The first axis (PCb1) captured the dominant proportion of variation in the beak dataset (89% of the total variance) All factor loadings were negative and >0.40. Lower values of PCb1 indicate species with relatively larger and less downward inclined beaks, the great-billed parrot Tanygnathus megalorynchos being the species with the lowest score. Higher PCb1 values indicate species with relatively smaller and more curved beaks (and high bite force) like the Papuan lorikeet Charmosyna papou or the superb parrot Polytelis swainsonii. Species’ scores on these two morphological axes (PCa1 -appendices- and PCb1 -beak shape-), body size and HWI values were used as input in subsequent analyses.
We defined five main trophic groups (1. omnivore; 2. granivore; 3. frugivore; 4. herbivore; 5. nectarivore) according to dietary information obtained from the literature [57, 66-67]. The diet of Australasian parrots is reasonably well documented and we found no discrepancies among sources. Subsequently, we defined a broader category in which we discerned between species that include nectar in their diet and those that do not (Fig. 1).
Lineage diversification
We generated a lineage-through-time (LTT) plot to test whether cladogenesis in this group fits pure-birth (constant-rate speciation), early-burst, or recent speciation models. The LTT plot includes the observed LTT line and the 95% confidence intervals (CI) from our sample of 500 trees from BirdTree. In conjunction, we estimated the Pybus & Harvey’s γ statistic, for which values significantly different from zero indicate strong departures from the constant-rate null model (i.e., diversification is a function of time and space) [68]. A deceleration in diversification rate is inferred when γ≤ -1.645 [69], whereas positive gamma values are associated to LTT plots with slower initial diversification followed by a rapid accumulation of recent lineages [70]. We computed this metric using both the MRP tree and our sample of 500 trees. Lastly, we also determined whether the rate of diversification remained constant over time following a pure-birth (Yule) or a birth-death process (pb and bd, respectively), changed over time according to a two- or three-rate Yule process (yule2 and yule3, respectively), or decreased in accordance with either of the two diversity-dependent models (DDL = linear decrease, DDX = exponential decrease). Maximum likelihood models of diversification were fitted using the MRP tree and our sample of 500 trees from BirdTree. We also ran these analyses only including the clade of lories (Loriinae), which comprises most nectarivorous species. Diversification analyses were implemented using the R packages ‘ape’ [71], ‘geiger’ [72], and ‘laser’ [73].
We tested for heterogeneity in diversification rates across the phylogeny using Bayesian analysis of macroevolutionary mixtures (BAMM) [74-76]. This program is oriented entirely towards detecting and quantifying heterogeneity in evolutionary rates, allowing to infer mixtures of time-dependent and clade-specific lineage or phenotypic rate regimes on phylogenetic trees. BAMM uses a reversible jump Markov chain Monte Carlo algorithm to explore the universe of candidate cladogenesis models. Shifts in evolutionary rates are detected automatically, with no a priori designations, and can occur at nodes or along branches. After computing tree-appropriate rate priors using the R package BAMMtools [76], we ran BAMM for 100 million generations on the MRP tree, sampling every 20,000 generations. To account for incomplete sampling, we specified the percentage of species that have been sampled by setting the ‘globalSamplingFraction’ parameter in our control file. We used the ‘coda’ package [77] to explore the Markov chain Monte Carlo (MCMC) output and check for convergence (50% generations as burn-in).
Macroevolutionary modeling
(a) Disparity
We analyzed the temporal dynamics of traits’ disparity using disparity-through-time (DTT) plots. This method computes the phenotypic disparity of each subclade relative to the phenotypic disparity of all taxa, and plots the change in average relative subclade disparity through time [72]. The observed DTT trajectory was compared to that expected according to pure Brownian motion (2,500 simulations) and deviations were evaluated using the Morphological Disparity Index (MDI), which represents the area between the curves illustrating the observed and expected patterns of morphological divergence. Disparity values significantly greater than the null suggest that subclades overlap, and all contain a significant proportion of variation found throughout the entire group at a given time. Negative disparity indicates that morphological variation is partitioned among subclades, indicating that each subclade occupies a distinct region of the morphospace. In adaptive radiations, disparity is highly partitioned among clades early in the radiation, with each lineage holding little of the total disparity, and significant negative MDI values are expected. Trait disparity-through-time (DTT) analyses were conducted using the R package ‘geiger’ [78] and the rank envelope test developed by [79]. This test computes a range of p‐values that encompasses the most liberal and most conservative p‐values, respectively, and shows better type I and II error rate than the original pointwise envelope test (see [79] for more details).
(b) Trait model fitting
We then assessed several evolutionary models to find the best fit to explain the evolution of both overall morphology (body size, PCa1 and HWI) and beak morphology (PCb1) in Australasian parrots. Firstly, we fitted a single-rate Brownian motion (BM) model. This model assumes a random walk (time-homogeneous) process in which individual trait values fluctuate at random and increase uniformly as a function of time and thus, indicates undirected and unconstrained stochastic change [80]. That is, traits change as a result of genetic drift (neutral evolutionary change). The dynamics of this model are governed by a single parameter (σ2) that determines the rate at which independent lineages diffuse away from an ancestral condition. A non-uniform multi-regime BM model (BMS) that allows σ2 to differ among states, and where the rubber-band parameter α is set to zero (diffusive model of evolution; [81]) and a directional random walk subject to decreasing rates through time called ‘early-burst’ EB model (also known as ‘ACDC’ model) [34] were also fitted. In addition, we incorporated two modifications of the uniform BM process that capture other conceivable ways in which traits might evolve. Specifically, a model (delta) with an added parameter (δ) to model a rate shift indicating a late or early-burst fashion (i.e., whether trait evolution has accelerated or slowed down over time), and a model that evaluates whether species traits have evolved according to phyletic gradualisms or punctuated evolution through history (kappa model) [82].
Increasing in complexity, the Ornstein-Uhlenbeck (OU) process is a modification of BM that incorporates attraction (α) to a trait ‘optimum’ (θ). OU models describe the evolution of a particular trait towards or around a stationary peak or optimum value, at a given evolutionary rate (i.e., constrained morphological evolution) [14]. We expected that the rate at which a trait changes through time (σ2) would be higher (i.e., faster evolution) in more specialized niches compared to broader generalist niches [83]. We also examined multipeak OU models in which these parameters vary with dietary regime allowing us to estimate adaptation and strength of selection towards separate phenotypic optima. We considered an OU model with attraction to five distinct selective optima corresponding to different dietary groups (OU5 = OUdiet) and another two-regime model with different optima for nectarivorous and non-nectarivorous species (OU2 = OUnectar). In these models, different peaks can only be discovered by certain lineages (which are inferred to be under the same regime), and do not affect the evolutionary dynamics of other regions of the tree. As an alternative, we employed the Fokker-Planck-Kolmogorov (FPK) model recently devised by, which can approximate adaptive landscapes with more complex shapes than those assumed by OU models, thereby accommodating scenarios of disruptive and directional selection. According to the FPK model, the trait evolves to random diffusion and is also subject to an ‘evolutionary potential’ (V) that creates a force that pulls the trait towards specific regions of the trait interval. That is, the trait evolves according to bounded random diffusion (or biased random-walk). Hence, while OU models with several optima might be better at describing transitions (i.e., situations in which a lineage shifts to another adaptive zone [15], multi-peak FPK might represent more genuine diversifying selection towards alternative phenotypic optima [19]. In other words, the OU model is specifically devised to model adaptation in the form of stabilizing selection around an optimal value [13]. In contrast, the FPK model and extended versions (Bounded Brownian Motion, BBMV) are useful to describe a scenario in which traits can drift freely over time in phenotypic space, except that they are confined between two bounds [84]. Discerning between both cases is of utmost importance since several studies have shown that neutral evolution between bounds produces patterns closely resembling the ones obtained under an OU process [85].
Finally, to assess the possible role of pulsed processes in trait evolution, we fitted Lévy process models. Specifically, a pure normal inverse Gaussian process (NIG), which captures constant rapid adaptation [86], and two evolutionary jumps processes that model stasis followed by sudden shifts in trait values between adaptive zones. The evolutionary jumps models were implemented following the approach described in [87], which assumes that jump effects are normally distributed (pure jump-normal model; JN), and a variant of this approach (JEM) that uses an alternative optimization that does not require matrix inversions to sample jump configurations (see [32] for further details). These pulsed models are intended to evaluate the punctuated component of evolutionary processes and were not assessed in conjunction with the conventional Gaussian models (and its modifications) described above (BM, EB and OU), which depict evolutionary patterns (gradual change, gradual stasis, and adaptive radiation) that differ in form and mode (see e.g., [52]).
Model fitting was carried out using the packages ‘geiger’ [72], ‘mvMORPH’ [81], ‘BBMV’ [54] and ‘pulsR’ [86] in R 4.0.2, as well as with the software ‘levolution’ (https://bitbucket.org/wegmannlab/levolution, commit d9898f4, [32]). The relative fit of models was assessed using the Akaike’s Information Criterion corrected for small sample sizes (AICc) and the AICc metric. To consider uncertainty in the ancestral dietary regime, all multiregime models were run on 500 stochastically mapped trees generated using stochastic character mapping (SIMMAP), an empirical Bayes procedure whereby character substitutions are mapped onto the nodes and branches of a phylogenetic tree. This MCMC approach whereby character substitutions are mapped onto the nodes and branches of a phylogenetic tree was implemented using the make.simmap function in ‘phytools’ [88]. The best-fit model of character evolution for the evolution of nectarivory (2-regimes) and diet (5-regimes) was determined by fitting an equal-rates (“ER”) model, a symmetric (“SYM”) model, an all-rates-different (“ARD”) model to the data set using the function fitDiscrete in ‘geiger’. The ER and the SYM transition-rate model were employed to estimate the ancestral states of nectarivory and diet, respectively, as these models yielded the lowest AICc [89]. Complementarily, we used simulation-based phylogenetic ANOVA [90] to test for differences in morphological traits among diet categories. We also used ‘OUwie’ (Beaulieu & O’Meara, 2015) to estimate trait optima and rate of evolution for each regime.
Phenotypic diversification rates
To estimate rates of evolution for continuous traits (body size, HWI, PCa1 and PCb1), we used BAMM in a similar way to that previously described. We performed these analyses across the entire phylogeny and for lories separately in order to determine whether Loriinae, the clade in which nectarivory is the most frequent feeding strategy, underwent an increase in diversification rate relative to other psittaculid lineages. Evolutionary rate dynamics were visualized from BAMM output using phylorate and rate-through-time plots generated in BAMMtools (Rabosky et al., 2014).
In addition, we explored the statistical patterns of the trait data using the reversible-jump algorithm implemented in ‘bayou’ [20]. ‘bayou’ employs Bayesian rjMCMC to identify regimes without a priori hypotheses in phylogenetic comparative data. This Bayesian approach alleviate many of the challenges inherent to fitting OU models, due to ridges in the likelihood space. By examining the posterior distribution, rather than point estimates, this method provides a more realistic interpretation of evolutionary models. Convergence of optima in ‘bayou’ was assessed by visual comparison of data within the phenogram.