The upper and lower jaws of some wrasses (Eupercaria: Labridae) possess teeth that have been coalesced into a strong durable beak that they use to graze on hard coral skeletons, hard-shelled prey, and algae, allowing many of these species to function as important ecosystem engineers in their respective marine habitats. While the ecological impact of the beak is well-understood, questions remain about its evolutionary history and the effects of this innovation on the downstream patterns of morphological evolution.  Here we analyze 3D cranial shape data in a phylogenetic comparative framework and use paleoclimate modeling to reconstruct the evolution of the labrid beak across 205 species. We find that wrasses evolved beaks three times independently, once within odacines, and twice within parrotfishes in the Pacific and Atlantic Oceans. We find an increase in the rate of shape evolution in the Scarus+Chlorurus+Hipposcarus (SCH) clade of parrotfishes likely driven by the evolution of the intramandibular joint. Paleoclimate modeling shows that the SCH clade of parrotfishes rapidly morphologically diversified during the middle Miocene. We hypothesize that possession of a beak in the SCH clade coupled with favorable environmental conditions allowed these species to rapidly morphologically diversify.

Shape Analyses

For the analysis of skull shapes, we used three-dimensional geometric morphometrics to quantify shape variation and diversity among our labrid dataset. We digitized the left side of each specimen with 79 landmarks and 118 semi-landmarks (Supplementary Table 1; Supplementary Figure 4) following the approach described in Larouche et al. (2022). Landmarks encompassed the skull, jaws, hyoid region, and pharyngeal jaws. The teleost fish skull exhibits immense biomechanical complexity, with highly kinetic, articulating elements (Westneat 2004; Hulsey et al. 2005). This kinesis poses a challenge to studies of shape change across the skull because preservational artifacts related to the relative positions of these individual elements can strongly bias any downstream analyses (Vidal‐García et al. 2018; Evans et al. 2019b). To account for this rotation and translation of mobile elements, we performed a local superimposition to standardize the position of the different skull elements (Rhoda et al. 2021a; Rhoda et al. 2021b). Specimen positions were standardized to the specimen that was closest to the mean shape (Ophthalmolepis lineolata) using the findMeanSpec function in geomorph. To study the shape evolution of the individual beak elements, we subset our larger skull shape dataset into smaller premaxilla and dentary datasets, the raw coordinates of which were individually superimposed. After the local superimposition, a principal component analysis (PCA) was performed to assess the primary axes of shape variation across all three of our wrasse datasets. All geometric morphometric analyses were performed in the R-package geomorph version 4.0.3 (Adams et al. 2016).

Morphological Sampling

Skull shape was sampled across 234 wrasse specimens representing 204 different species of wrasse (25% taxon sampling). For most species, sampling was limited to a single adult individual per species due to limitations associated with the overall size of many of the large parrotfish species and the time associated with the data collection process (Supplementary Table 2). For our shape analysis we excluded a sub-adult Bolbometopon muricatum specimen. The lack of intraspecific sampling in our dataset prohibits the ability to estimate morphological variation within species and can inflate the estimation of rates of morphological evolution between species. However, studies using geometric morphometric data have shown that across datasets where sample variance is high (as is the case with our dataset) shape relationships between species tend to be accurately estimated with even a low number of individuals per-species (Cardini and Elton 2007). To further demonstrate this, we performed an intraspecific analysis using five individuals of Chlorurus spilurus that were all collected from the same place (Moorea, French Polynesia) around the same time (within three weeks of each other, June 2021) and show that the specimens of this species have significantly (p=0.006) less Procrustes variance (0.00015), than the other species in our dataset (0.00078) (Supplementary Figure 5). We used micro-computed tomography scanning (micro-CT scanning) to collect three-dimensional osteological data for each specimen. Scans were conducted at Rice University, the University of Minnesota, the Uniersity of Chicago and the University of Washington Friday Harbor Labs in conjunction with the oVert initiative. Scans were segmented in Amira™ to isolate the skull, and exported as three-dimensional mesh files. Mesh files were then imported into Stratovan Checkpoint ™ where they were digitized.

Phylogenetic Comparative Methods

We analyzed patterns of skull and jaw shape variation and shape diversification using a phylomorphospace analysis (Sidlauskas 2008), which displays the principal component scores of each species with an underlying phylogeny and estimates the skull shapes at ancestral nodes.

Beaks are typically considered to be discrete characters and defined by the fusion or coalescence of teeth in the upper and lower jaws. However, for our analyses, we were also interested in the overall shape of the beak among labrid fishes. To determine whether shape could predict the presence or absence of beak morphology, we tested for the relationship between beak presence and skull shape using a Procrustes analysis of variance (ANOVA) and a Procrustes phylogenetic generalized least squares (PGLS) (procD. pgls function in geomorph) analysis to account for phylogenetic non-independence of our shape data (Adams and Collyer 2018) . We performed these analyses on the entire skull shape configuration, as well as the individual dentary and premaxillary datasets. We recover significant relationships across all Procrustes ANOVAs but not for PGLS analyses suggesting that shape is strong predictor of beak presence, but that this relationship is phylogenetically restricted, which implies that some beaked species may have evolved beaks that differ substantially in shape from other beaked species in different clades (Supplementary Tables 3:8).

To reconstruct the evolution of beaks across space and time in labrids and specifically parrotfishes, we estimated the ancestral parrotfish skull shape by warping a mesh of the specimen closest to the shape mean (Ophthalmolepis lineolata) to the ancestral parrotfish node using geomorph.

Ancestral State Reconstruction

The ancestral state of beaks for 205 wrasse species was estimated using the R-package phytools version 1.0-1(Revell 2012). Species were designated as beaked or non-beaked following visual inspections via micro-CT scans. We used stochastic character mapping (Huelsenbeck et al. 2003; Bollback 2005) to reconstruct the evolutionary history of the presence or absence of a beak. Stochastic character mapping was performed using the make.simmap function in the R package phytools (Revell 2012). The maximum likelihood Q-matrix was estimated from the data, and the prior for the root state was determined using its stationary distribution, conditional on the Q-matrix. We used the fitDiscrete function from geiger to determine the best-fitting model of character evolution between and equal rates (ER), symmetrical (SYM)and all-rates-different (ARD) model. To account for phylogenetic uncertainty, we fit different models of discrete trait evolution across 100 randomly sampled trees from the posterior distribution of the BEAST run used to build the consensus phylogeny. We also used these randomly sampled trees for our stochastic character mapping analysis and estimated the ancestral state at each node using posterior probabilities summarized across 1000 simulations for each of the 100 phylogenies. For the evolutionary history of the presence or absence of a beak, the lowest AIC values were obtained for the ER and SYM models (as transition rate matrices are identical for binary characters), however these models were not found to be substantially better fitting compared to an ARD model with a ΔAIC of only 0.93. We conservatively chose to perform the stochastic character mapping for the beak presence/absence using an ER model of trait evolution.

Assessing Convergent Evolution Hypotheses

After recovering support for multiple independent origins of beaked dentition among labrids, we tested for whether these beaked species evolved similarly-shape skulls and beaks via convergence, or if the shapes of the skull and beak are entirely distinct among the different taxa. We used the distanced-based metrics of Stayton (2015) (C1-C4) to test for convergence between beaked wrasses. We also used these metrics to quantify the degree of convergence between these species. We ran our convergence analysis on the entire skull configuration as well as on the premaxilla and dentary bones separately.

Quantifying Rates of Shape Evolution

To test the effect of beaks on rates of morphological diversification, we used the compare.evol.rates function in geomorph which calculates the rate of shape evolution between groups of specimens under a Brownian motion model (Denton and Adams 2015). To account for phylogenetic uncertainty in our rate estimates, we ran this analysis over 100 randomly sampled phylogenies from the posterior distribution of the BEAST run and used a paired t-test to test for differences in the mean evolutionary rate between beaked and non-beaked labrids following the approach of Evans et al. (2019a). We also estimated the branch-specific rates of skull shape evolution for 204 labrid species. Rates were estimated using a variable rates model implementation in the BayesTraitsV4 program (Venditti et al. 2011). This Bayesian method uses a reversible-jump MCMC chain approach to estimate the probability of rate shifts in trait data across a phylogeny revealing clade and species-specific rate shifts in trait data. To reduce the dimensionality of our data, we used the first 23 principal components (PCs) because they accounted for 85% of our total shape variance. While PC axes are mathematically orthogonal, and thus uncorrelated, trait variation can still be evolutionarily correlated. To account for this, we ran our analyses using the “TestCorrel” function which constrains the correlation between trait axes to zero. We used uniform, uninformative priors and ran four independent chains each for 200,000,000 generations discarding the first 60,000,000 as burn-in. The chain was sampled every 1,400,000 generations after convergence using a stepping-stone. Model convergence was evaluated for each model by running the analysis a second time and visually assessing the trace of the marginal likelihoods using Tracer (Rambaut et al. 2018). We evaluated a “variable rates” model that allows for rate heterogeneity and identifies regions of the tree where evolutionary rates differ across different branches and internal nodes (Venditti et al. 2011). The resulting output of the variable-rates analysis is a set of phylogenies where each branch is scaled by its Brownian motion rate of evolution. The BayesTraits analyses were run for the entire skull shape configuration. Additionally, due to the fact that many of the species in our analysis are represented by a single specimen, we ran the BayesTraits analysis using a lambda transformation of 0.754 estimated from the trait data which lengthens the branches at the tips of a phylogeny, accounting for measurement error, and within-species variation (Baken et al. 2021; Goswami et al. 2022). We also performed a variable rates regression in BayeTraits following the approach of Kubo et al. (2019) and Baker and Venditti (2019) using the same parameters mentioned above. Briefly, this method simultaneously estimates shifts in the rate of trait evolution using the variable rates model while also estimating the parameters of a phylogenetic regression. This method has the advantage of accounting for background rates of trait evolution, which has been shown to be important for accurately estimating the effect of discreet characters on rates of trait evolution (May and Moore 2020).

Analyses were run in R version 4.1.3 and BayesTraits Versions 3 and 4.