DOI: 10.5061/dryad.prr4xgxrz
title: Data from: Beaks promote rapid morphological diversification along distinct evolutionary trajectories in labrid fishes (Eupercaria: Labridae)

description: 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.
description: 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).
description: Analyses were run in R version 4.1.3 and BayesTraits Versions 3
 and 4.