DOI: 10.5061/dryad.gmsbcc2x4
title: Data from: A classic key innovation constrains oral jaw functional diversification in fishes

description: Modifications to the pharyngeal jaws – a prey processing system located
 posterior to the mouth cavity – are widely considered a key innovation
 that enhanced diversification within several prominent fish clades. Seen
 in cichlids, damselfishes, wrasses, and a few other lineages, these
 musculoskeletal alterations are believed to increase the evolutionary
 independence, and thus the diversification of the oral and pharyngeal jaw
 systems. To test this classic hypothesis, we conducted comparative
 phylogenetic analyses to assess the effect of the pharyngeal novelty on
 the diversification of feeding morphology and kinematics across a
 taxonomically diverse sample of spiny-rayed fishes. We quantified
 movements of the oral jaws and other craniofacial structures from 689
 suction feeding strikes using high speed videos collected from 228 species
 with and without the pharyngeal jaw novelty. Contradicting long-held
 predictions, we find significantly greater disparity across all traits and
 faster rates of oral jaw functional evolution in fishes without the
 specialized prey processing system. The modified pharyngeal jaw is
 undoubtedly a functional innovation as it enhances the strength of the
 prey processing system, facilitating exceptional transition rates to
 feeding on hard and tough prey. However, it also restricts the
 diversification of the feeding system, revealing that the impact of
 pharyngognathy is more nuanced than previously thought. In light of these
 and other recent findings, a reinterpretation of the macroevolutionary
 consequences of the pharyngeal jaw novelty is needed.
description: <strong><em>Feeding videos & morphometric
 data</em></strong> Totaling 689 feeding
 strikes, this dataset spans a wide swath of the diversity of suction
 feeding motions in spiny-rayed fishes (Acanthomorpha; Supplemental Table
 1). Specimens of 133 and 95 species with and without modified pharyngeal
 jaws, respectively (Supplemental Table 2), were obtained from commercial
 sources and housed in laboratory aquaria at room temperature where they
 were filmed feeding on mobile prey. We sampled a total of 43 acanthomorph
 families, including wrasses and parrotfishes (Labridae including
 Scaridae), damselfishes (Pomacentridae), the only false scorpionfish
 (Centrogenyidae), and cichlids (Cichlidae) to represent four of the seven
 clades in which all known species exhibit a modified pharyngeal jaw
 (Figure 1B). Though our sampling of MPJ fishes does not
 include surfperch (Embiotocidae), halfbeaks (Hemiramphidae), or flying
 fishes (Exocoetidae), previous work exploring the ecomorphological
 diversity in six taxon-rich MPJ families found that these three clades
 show reduced body shape diversity and, on average, slower rates of body
 shape evolution compared to their non-MPJ sister clades. Further, these
 three MPJ families overlap significantly in body shape and/or ecological
 space with cichlids, damselfishes, and wrasses. Halfbeaks and flying
 fishes – species typically more elongate in body and head shape – do
 occupy distinct body shape space compared to other MPJ fishes. However,
 authors found (1) complete overlap in diet space, (2) 1.4 to 3.1-fold less
 morphological disparity, and (3) significantly slower rates of evolution
 in seven out of eight body and head shape traits among these two clades
 compared to the four other MPJ clades included in the study (Larouche
 <em>et al.</em>, 2020). Thus, the addition of halfbeaks,
 surfperch, and flying fishes may have affected our estimates of
 multivariate variance in MPJ fishes, but it is unlikely that their
 inclusion would significantly impact our overall findings.  Additionally,
 we note that our sampling of Cichlidae only includes Neotropical and
 African riverine cichlids. Recent examination of the diversification of
 prey capture mechanisms across Neotropical and African Rift Lake cichlids
 found (1) that Neotropical cichlids have the most variance in craniofacial
 morphology, but similar or lower variance for most functional traits
 compared to cichlids from Lake Tanganyika; and (2) a negative relationship
 between age and rates of evolution such that young lake radiations display
 faster diversification than the older Neotropical cichlid radiation
  (Martinez <em>et al.</em>, 2024). Thus, the exclusion of Rift
 Lake cichlids may reduce our estimates of functional diversity and rates
 of evolution. Still, the Neotropical and African riverine cichlids
 included herein represent much of the ecomorphological and functional
 variation found throughout Cichlidae (Burress, 2014; Burress <em>et
 al.</em>, 2017; Martinez <em>et al.</em>, 2024).
 Further, by excluding Rift Lake cichlids, we avoid conflating the impact
 of the MPJ with the influence of other factors, such as hybridization
 (Joyce <em>et al.</em>, 2011; Meier <em>et
 al.</em>, 2017, 2019, 2023) and lake effects (Seehausen, 2006)which
 have been shown to drive increased evolutionary rates in African Rift Lake
 cichlids. Feeding strikes were filmed with a digital
 Fastec HiSpec 2G Monochrome high-speed camera at 2,000 frames per second
 and in 592 x 474 resolution. To ensure consistency in motion tracking
 across all strike sequences, videos were recorded from a lateral
 perspective and were only used if fishes exhibited a full-effort suction
 feeding attempt. We used landmark morphometrics to track oral jaw and
 craniofacial motions during each suction feeding sequence using 10 frames
 spanning the feeding strike. After identifying the video frames depicting
 the start (onset of mouth opening; designated as frame 1) and end (maximum
 oral jaw and mouth cavity expansion, prior to mouth closure; designated as
 frame 10) of each strike, we selected eight additional frames that were
 equally spaced in time between these two frames (Figure 2). We used either
 tpsDIG2 software (Rohlf, 2015) or the digitizeImage function in
 STEREOMORPH v.1.6.7 (Olsen & Westneat, 2015) in
 <em>R</em> v.4.2.2 (R Core Team, 2022) to digitize 18
 landmarks on feeding-related structures in each of the ten frames to
 capture the oral jaw motions and craniofacial shape change that occurred
 during each feeding sequence (Supplemental Figure 1A). Ten homologous
 fixed landmarks were used to track the motions of several oral jaw bones
 and the skull. Eight sliding semi-landmarks, bound by fixed landmarks on
 the insertion of the pelvic fin and a ventro-anterior point on the
 mandible, captured shape change along the ventral margin of the head due
 to buccal expansion. We used published landmark data for 155 species of
 cichlids (Martinez <em>et al.</em>, 2024) and reef fishes
 (Corn <em>et al.</em>, 2021). Landmark data for the remaining
 73 species were new data collected for this study.
 Using landmark coordinates, we calculated six ‘motion
 components’, including premaxillary protrusion, hyoid depression, oral jaw
 gape, maxillary rotation, lower jaw rotation, and cranial rotation
 (Supplemental Figure 1B). These traits describe the rotational and linear
 motions of oral jaw bones and the skull during a suction feeding sequence.
 To ensure that all motion components were commensurate, the rotational
 variables were converted to linear distances as outlined in a recent study
 (Martinez <em>et al.</em>, 2024). To account for the effect of
 size on feeding mechanics, each component was scaled by the centroid size
 of the fish’s head when in the closed mouth position. All components were
 then log-transformed to achieve normal distributions. Finally, for each of
 the six motion components, we computed average trait values by specimen
 before calculating average trait values for each of our 228 species (on
 average, n=3 strikes per species; Supplemental Figure 3).
 In addition to examining diversification in the components of
 prey capture motions, we explored the evolution of two composite kinematic
 traits and cranial morphology. Using the gpagen function in GEOMORPH
 v.4.0.5 (Collyer & Adams, 2018; Baken <em>et al.</em>,
 2021; Collyer & D. C. Adams., 2021; Adams & Collyer, 2022;
 Adams <em>et al.</em>, 2023), we performed a generalized
 Procrustes analysis (GPA) which aligns morphometric data through an
 iterative process of scaling, rotating, and translating landmarks to
 minimize the Procrustes distances across motion shapes120. Aligned shapes
 create a trajectory through morphospace for each motion (Adams, 2014a;
 Figure 2A), and a trajectory’s length provides a metric for ‘total
 craniofacial kinesis’ – a univariate trait describing the magnitude of
 shape change that occurs due to movement of the oral jaw system’s
 anatomical features during a feeding sequence (Martinez <em>et
 al.</em>, 2018). Kinesis was computed for each prey capture event by
 summing the Procrustes distances between each pair of successive motion
 shapes (Figure 2B; Martinez et al., 2018; Corn et al., 2021). Kinesis skew
 was then calculated by taking the natural logarithm of kinesis across the
 last five motion shapes divided by total craniofacial kinesis (Figure 2B;
 Supplemental Figure 1A), providing information about the distribution of
 movement across a strike. Total craniofacial kinesis and kinesis skew
 values were computed for each of the 689 feeding strikes, before averaging
 by specimen and then by species (Supplementary Table 2). Finally, from
 each species' aligned and average morphometric data, we extracted the
 starting head shape for each strike—the sole morphological trait included
 in this study—which describes craniofacial form when in the closed mouth
 position at the start of a feeding sequence (Supplementary Table
 2) <strong><em>Morphological and functional
 disparity</em></strong> To conduct our
 analyses in a phylogenetic context, we trimmed a time-calibrated molecular
 phylogeny of ray-finned fishes (Rabosky, Chang, Title, Alfaro et al.,
 2018) to the species in our dataset (Figure 1; Supplemental Table 2) using
 APE v.5.7-1 (Paradis & Schliep, 2019). For species in our
 kinematic dataset that were not present on the phylogeny, we used the
 closest related species within the same genus as proxy. If no congeners
 were present, we selected a substitute species from within the taxon’s
 clade at random, excluding any species sister to those in our dataset. In
 total, 30 substitutions were made (Supplemental Table 2).
 To visualize the multivariate data in shape and kinematic space,
 we performed separate principal component analyses (PCAs) of (1) aligned
 morphometric data from all ten frames, showcasing feeding sequences as
 trajectories of shape change; (2) aligned morphometric data from the first
 frame in the feeding sequence, which captures head shape in the starting,
 closed-mouthed posture; and (3) the six, size-scaled motion components.
 Here, PCAs were performed using covariance matrices in the gm.prcomp
 (Revell, 2009; Collyer & Adams, 2021) and prcomp (Mardia
 <em>et al.</em>, 1979; Becker <em>et al.</em>,
 1988; Venables & Ripley, 2002) functions within GEOMORPH and STATS
 v.4.2.2 (R Core Team, 2022) to analyze morphometric and motion component
 data, respectively. To compare the magnitude of interspecific variance
 among MPJ and non-MPJ morphologies and kinematic patterns, we used the
 full-dimensional trait data and the morphol.disparity function (Zelditch
 <em>et al.</em>, 2012; Collyer & Adams,
 2021)<em> </em>in GEOMORPH to compute disparity – the sum of
 the diagonal elements of the variance-covariance matrix divided by the
 number of observations (Zelditch <em>et al.</em>, 2012) – for
 interspecific head shape, the six motion components (both combined in a
 multivariate dataset, and individually as univariate traits), total
 craniofacial kinesis, and kinesis skew. To test whether the MPJ and
 non-MPJ fishes have statistically different mean craniofacial
 morphologies, we performed a phylogenetic multivariate analysis of
 variance (phylogenetic MANOVA) on the full-dimensional starting head shape
 data over 10,000 iterations under a Brownian Motion model with the
 procD.pgls function (Adams, 2014b; Adams & Collyer, 2015, 2016,
 2018; Collyer <em>et al.</em>, 2015) in GEOMORPH. This same
 method was used to perform individual phylogenetic ANOVAs on the
 full-dimensional trait data for each of the six, individual motion
 components, total craniofacial kinesis, and kinesis skew to determine if,
 on average, there are significant differences in kinematic processes
 employed by MPJ and non-MPJ fishes.
 <strong><em>Patterns of trait
 coevolution</em></strong> To understand the
 impact of the modified pharyngeal jaw on patterns of trait coevolution, we
 examined (1) the overall strength of evolutionary integration for cranial
 morphology and jaw function and (2) pairwise trait correlations in MPJ and
 non-MPJ species. To do this, we first assessed the degree of evolutionary
 integration separately for interspecific head shape and motion components
 using the integration.Vrel function (Pavlicev <em>et al.</em>,
 2009; Conaway & Adams, 2022)in GEOMORPH. We then used the
 compare.ZVrel function (Conaway & Adams, 2022) in GEOMORPH,
 calculating covariances for head shape and motion components in each group
 (i.e., their relative eigenvalue variance, <em>Vrel</em>) and
 converting them<em> </em>to a standardized effect size. The
 function then compares these effect sizes using a two-sample test to
 determine if there is a significant difference in the strength of
 integration between MPJ and non-MPJ species. Second, we examined pairwise
 relationships across all functional traits (i.e., kinesis, kinesis skew,
 and all six motion components) to further understand if and how fishes
 with the pharyngeal jaw novelty use different mechanical strategies to
 successfully capture prey. We first produced evolutionary correlation
 matrices by calculating pairwise correlations between phylogenetic
 independent contrasts for each functional trait using the pic and
 corr.test functions in APE and PSYCH v.2.4.3 (Revelle, 2024),
 respectively. We then removed pairwise correlations that were not stronger
 than expected under Brownian motion. Here, we used the fastBM function in
 PHYTOOLS v.1.5-1 (Revell, 2012) to simulate trait data 1,000 times under a
 Brownian Motion model of evolution, bound by each trait’s empirical
 maximum and minimum values. As done with the empirical data, we estimated
 the pairwise relationships between phylogenetic independent contrasts of
 the simulated data, resulting in null distributions of 1,000 simulated
 correlations for each pairwise trait comparison. Empirical correlations
 that did not fall at or above the 95<sup>th</sup> percentile
 of simulated correlations (14 of 56 pairwise relationships) were removed.
 We visualized each group’s evolutionary correlation matrix and the
 absolute value difference between them using the corrplot function
 (Murdoch & Chow, 1996; Friendly, 2002) in CORRPLOT v.0.92 (Wei
 & Simko, 2021). The strength of pairwise correlations in MPJ and
 non-MPJ fishes was statistically compared using the paired.r function in
 PSYCH. <strong><em>Rates of morphological
 and functional evolution</em></strong> We
 examined the impact of the modified pharyngeal jaw system on rates of
 morphological and functional evolution of feeding structures using a
 Bayesian approach. After coding each species for the presence or absence
 of the MPJ, we tested for an effect of this discrete trait on
 interspecific head shape and kinematic traits using a relaxed clock,
 state-dependent model of multivariate Brownian motion evolution
 implemented with the MuSSCRat model (May & Moore, 2020) in
 RevBayes (Hohna <em>et al.</em>, 2016). MuSSCRat uses
 continuous time Markov chain and Markov chain Monte Carlo processes to
 simultaneously estimate evolution of discrete and continuous traits,
 respectively, to determine the impact of the discrete character on
 continuous character evolution. In its computation of the rate parameter
 for multivariate data that represent complex, functionally integrated
 systems and processes, this model accommodates branch-specific rate
 variation (i.e., rate heterogeneity) using a relaxed-clock model. Further,
 this model accounts for variation in rate that is not due to the discrete
 character (referred to as ‘background rate variation’).
 We ran separate MuSSCRat analyses for (1) scores for the first 10
 principal component axes from a PCA of interspecific head shape, (2) the
 six motion components traits, (3) the univariate total craniofacial
 kinesis trait, and (4) the univariate kinesis skew trait (Supplemental
 Figure 2; Supplemental Table 2). We estimated evolutionary rates for
 interspecific head shape using scores from principal component axes that
 account for over 95% of shape variation due to computational limitations,
 as the model could not reach convergence on the morphometric data in its
 full dimensionality. In each of the four analyses, pharyngeal jaw
 configuration was used as the discrete character where the prior (i.e.,
 expected) number of transitions from the unmodified to modified state was
 set to 4 and the magnitude of the associated rate shift was drawn from a
 lognormal distribution. Each MCMC was run for 500,000 generations with 10%
 burn-in; the prior probability of state-dependent continuous character
 evolution was set to 0.5; and the priors for the branch-specific
 background rates were drawn from an uncorrelated lognormal distribution
 (i.e., UCLN relaxed clock). We determined that varying priors had little
 effect on posterior parameter estimates by repeating these analyses under
 different prior specifications for (1) the number of discrete character
 state transitions, (2) lognormal distribution of the discrete character
 rate shift, and (3) the probability of state-dependence (Supplemental
 Table 7). We used Tracer software v.1.7.2 (Rambaut <em>et
 al.</em>, 2018) to assess model convergence and code modified from
 REVGADGETS package v.1.1.0 (Tribble <em>et al.</em>, 2022)to
 visualize the results on our phylogeny. We assessed the
 degree of rate heterogeneity throughout the evolutionary history of each
 trait to determine if rates of functional and morphological
 diversification shift following transitions to a modified pharyngeal jaw.
 We performed local polynomial regressions between log-transformed, overall
 branch rate estimates from MuSSCRat and internal node ages to visualize if
 and how continuous character rates change through time. Additionally, we
 performed node height analyses using the nh.test function in GEIGER
 (Pennell <em>et al.</em>, 2014), which fit linear models
 between log-transformed, absolute values of phylogenetic independent
 contrasts (PICs) and internal node ages. Univariate node height analyses
 were run for kinesis, kinesis skew, each of the six motion components, and
 head shape scores for each of the first 10 PC axes. For multivariate
 analyses of head shape and motion components, we used means of the
 absolute values of PICs of principal component scores and individual
 motion components, respectively. A significant positive relationship would
 indicate increases in rate through time, while a significant negative
 relationship would indicate that rates have slowed and potentially be
 consistent with an early burst in trait evolution (Freckleton &
 Harvey, 2006; Slater & Pennell, 2014).
 <strong><em>Comparing rate estimates between
 maximum-likelihood and Bayesian frameworks
 </em></strong> Although MuSSCRat presents
 multiple advantages for estimating rates of evolution, we also estimated
 rates of continuous character evolution under a variety of Brownian Motion
 (BM) and Ornstein-Uhlenbeck (OU) models using the OUwie function within
 OUWIE v.2.10 (Hansen, 1997; Butler & King, 2004; O’Meara
 <em>et al.</em>, 2006; Thomas <em>et al.</em>,
 2006; Beaulieu <em>et al.</em>, 2012; Ho & C., 2014)
 to determine whether our findings are impacted by the underlying model of
 evolution. We generated a distribution of 100 stochastic character maps
 using PHYTOOLS, where the pharyngeal jaw condition was the discrete
 character. We then estimated evolutionary parameters under five models of
 evolution – single-rate Brownian Motion (BM1), multi-rate Brownian Motion
 (BMS), single rate and single optimum Ornstein-Uhlenbeck (OU1), single
 rate and multi-optima Ornstein-Uhlenbeck (OUM), and multi-rate and
 multi-optima Ornstein-Uhlenbeck (OUMV) – and compared the fit of these
 models based on AICc scores. While we recovered a mix of best-fitting
 models, generally, the OUwie results did not differ from the MuSSCRat
 findings – when MPJ and non-MPJ fishes show significantly different rate
 parameters (8 of 13 traits), non-MPJ fishes typically evolve faster (7 of
 8 traits; Supplemental Table 8). We note that the OUwie and MuSSCRat
 findings are not fully comparable as OUwie is only able to estimate
 univariate rates, whereas rates were estimated in multivariate ((A)
 interspecific head shape and (B) kinematic components; Supplemental Figure
 2) and univariate ((C) kinesis and (D) kinesis skew; Supplemental Figure
 2) formats within MuSSCrat. Because we do not recover significant,
 state-dependent evolution of kinematic components using MuSSCRat, it may
 be that individual motion components vary regarding which group evolves
 fastest, but that this variation averages out to a non-significant
 difference in state-dependent rates overall. Similar to findings in other
 recent studies (Corn <em>et al.</em>, 2021; Larouche
 <em>et al.</em>, 2023), our results suggest that estimates of
 evolutionary patterns (e.g., rate, integration, modularity) within a BM
 framework can still be informative even if a BM model is not the
 best-fitting evolutionary model. This also highlights some of the current
 limitations for evaluating models of multivariate trait
 evolution. In addition to the above analyses, we
 conducted a simulation study to examine the appropriateness of conducting
 maximum likelihood-based model-fitting analyses on this dataset. Here, we
 used the OUwie.sim function in OUWIE to simulate data under five models of
 evolution – BM1, BMS, OU1, OUM, and OUMV – across 100 stochastic character
 maps where the pharyngeal jaw condition was the discrete character. We
 then analyzed the fit of those same BM and OU models to each simulated
 dataset with the expectation that OUwie would return the model that the
 data was simulated under as the best-fitting model. In most cases (4 of 5
 simulated datasets), OUwie<em> </em>did not return the correct
 model and often preferred a more complex model (Supplemental Table 9).
 Thus, we conclude that a Bayesian approach (e.g., the MuSSCRat model) is
 more appropriate for evolutionary rate estimation within this dataset as
 our data does not have the power to accurately distinguish models of
 evolution within a maximum-likelihood framework where rates of discrete
 and continuous character evolution are not estimated jointly.
 <strong><em>Comparing parameter estimates between
 full and subsampled datasets </em></strong>
 To ensure that cichlids (n = 112 of the 133 MPJ species) do not
 drive the patterns recovered, we estimated evolutionary parameters using a
 subset of the cichlid taxa. Here, we randomly selected 13 cichlid species,
 matching the next highest number of species examined from a single MPJ
 family, resulting in a dataset of 95 non-MPJ and 34 MPJ taxa. Using the
 morphol.disparity and compare.evol.rates (Adams, 2014a; Denton &
 Adams, 2015; Adams & Collyer, 2018, 2019) functions in GEOMORPH ,
 we computed disparity and rates of evolution for interspecific head shape
 (all data used for disparity estimates; scores from PC axes 1-10 used for
 rate estimates following our treatment of the full dataset), motion
 components, total craniofacial kinesis, and kinesis skew. This analysis
 was repeated 100 times before comparing average disparity and rate
 estimates to those found using the full dataset of 133 MPJ species. On
 average, we found that estimates of functional and morphological disparity
 increased for all traits excluding kinesis skew which showed no change
 while rate estimates decreased for all traits excluding kinesis which
 decreased. Regardless of how parameter estimates shifted with a subsampled
 MPJ dataset, all patterns comparing MPJ to non-MPJ parameter estimates
 remained unchanged (Supplemental Figure 4). Thus, the results from the
 subsampling analysis align with our main findings where non-MPJ fishes
 show significantly greater disparity for all traits and a significantly
 greater rate of evolution for total craniofacial kinesis.
description: # Code & Data from: A classic key innovation constrains oral jaw
 functional diversification in fishes
 Article: [https://doi.org/10.1093/evlett/qrae046](https://doi.org/10.1093/evlett/qrae046) Data & Code: [https://doi.org/10.5061/dryad.gmsbcc2x4](https://doi.org/10.5061/dryad.gmsbcc2x4) ## Data Collection & Trait Computation ***Feeding videos & morphometric data*** Totaling 689 feeding strikes, this dataset spans a wide swath of the diversity of suction feeding motions in spiny-rayed fishes (Acanthomorpha; Supplemental Table 1). Specimens of 133 and 95 species with and without modified pharyngeal jaws, respectively (Supplemental Table 2), were obtained from commercial sources and housed in laboratory aquaria at room temperature where they were filmed feeding on mobile prey. We sampled a total of 43 acanthomorph families, including wrasses and parrotfishes (Labridae including Scaridae), damselfishes (Pomacentridae), the only false scorpionfish (Centrogenyidae), and cichlids (Cichlidae) to represent four of the seven clades in which all known species exhibit a modified pharyngeal jaw (Figure 1B). Though our sampling of MPJ fishes does not include surfperch (Embiotocidae), halfbeaks (Hemiramphidae), or flying fishes (Exocoetidae), previous work exploring the ecomorphological diversity in six taxon-rich MPJ families found that these three clades show reduced body shape diversity and, on average, slower rates of body shape evolution compared to their non-MPJ sister clades. Further, these three MPJ families overlap significantly in body shape and/or ecological space with cichlids, damselfishes, and wrasses. Halfbeaks and flying fishes – species typically more elongate in body and head shape – do occupy distinct body shape space compared to other MPJ fishes. However, authors found (1) complete overlap in diet space, (2) 1.4 to 3.1-fold less morphological disparity, and (3) significantly slower rates of evolution in seven out of eight body and head shape traits among these two clades compared to the four other MPJ clades included in the study (Larouche *et al.*, 2020). Thus, the addition of halfbeaks, surfperch, and flying fishes may have affected our estimates of multivariate variance in MPJ fishes, but it is unlikely that their inclusion would significantly impact our overall findings.  Additionally, we note that our sampling of Cichlidae only includes Neotropical and African riverine cichlids. Recent examination of the diversification of prey capture mechanisms across Neotropical and African Rift Lake cichlids found (1) that Neotropical cichlids have the most variance in craniofacial morphology, but similar or lower variance for most functional traits compared to cichlids from Lake Tanganyika; and (2) a negative relationship between age and rates of evolution such that young lake radiations display faster diversification than the older Neotropical cichlid radiation  (Martinez *et al.*, 2024). Thus, the exclusion of Rift Lake cichlids may reduce our estimates of functional diversity and rates of evolution. Still, the Neotropical and African riverine cichlids included herein represent much of the ecomorphological and functional variation found throughout Cichlidae (Burress, 2014; Burress *et al.*, 2017; Martinez *et al.*, 2024). Further, by excluding Rift Lake cichlids, we avoid conflating the impact of the MPJ with the influence of other factors, such as hybridization (Joyce *et al.*, 2011; Meier *et al.*, 2017, 2019, 2023) and lake effects (Seehausen, 2006)which have been shown to drive increased evolutionary rates in African Rift Lake cichlids. Feeding strikes were filmed with a digital Fastec HiSpec 2G Monochrome high-speed camera at 2,000 frames per second and in 592 x 474 resolution. To ensure consistency in motion tracking across all strike sequences, videos were recorded from a lateral perspective and were only used if fishes exhibited a full-effort suction feeding attempt. We used landmark morphometrics to track oral jaw and craniofacial motions during each suction feeding sequence using 10 frames spanning the feeding strike. After identifying the video frames depicting the start (onset of mouth opening; designated as frame 1) and end (maximum oral jaw and mouth cavity expansion, prior to mouth closure; designated as frame 10) of each strike, we selected eight additional frames that were equally spaced in time between these two frames (Figure 2). We used either tpsDIG2 software (Rohlf, 2015) or the digitizeImage function in STEREOMORPH v.1.6.7 (Olsen & Westneat, 2015) in *R* v.4.2.2 (R Core Team, 2022) to digitize 18 landmarks on feeding-related structures in each of the ten frames to capture the oral jaw motions and craniofacial shape change that occurred during each feeding sequence (Supplemental Figure 1A). Ten homologous fixed landmarks were used to track the motions of several oral jaw bones and the skull. Eight sliding semi-landmarks, bound by fixed landmarks on the insertion of the pelvic fin and a ventro-anterior point on the mandible, captured shape change along the ventral margin of the head due to buccal expansion. We used published landmark data for 155 species of cichlids (Martinez *et al.*, 2024) and reef fishes (Corn *et al.*, 2021). Landmark data for the remaining 73 species were new data collected for this study. Using landmark coordinates, we calculated six ‘motion components’, including premaxillary protrusion, hyoid depression, oral jaw gape, maxillary rotation, lower jaw rotation, and cranial rotation (Supplemental Figure 1B). These traits describe the rotational and linear motions of oral jaw bones and the skull during a suction feeding sequence. To ensure that all motion components were commensurate, the rotational variables were converted to linear distances as outlined in a recent study (Martinez *et al.*, 2024). To account for the effect of size on feeding mechanics, each component was scaled by the centroid size of the fish’s head when in the closed mouth position. All components were then log-transformed to achieve normal distributions. Finally, for each of the six motion components, we computed average trait values by specimen before calculating average trait values for each of our 228 species (on average, n=3 strikes per species; Supplemental Figure 3). In addition to examining diversification in the components of prey capture motions, we explored the evolution of two composite kinematic traits and cranial morphology. Using the gpagen function in GEOMORPH v.4.0.5 (Collyer & Adams, 2018; Baken *et al.*, 2021; Collyer & D. C. Adams., 2021; Adams & Collyer, 2022; Adams *et al.*, 2023), we performed a generalized Procrustes analysis (GPA) which aligns morphometric data through an iterative process of scaling, rotating, and translating landmarks to minimize the Procrustes distances across motion shapes120. Aligned shapes create a trajectory through morphospace for each motion (Adams, 2014a; Figure 2A), and a trajectory’s length provides a metric for ‘total craniofacial kinesis’ – a univariate trait describing the magnitude of shape change that occurs due to movement of the oral jaw system’s anatomical features during a feeding sequence (Martinez *et al.*, 2018). Kinesis was computed for each prey capture event by summing the Procrustes distances between each pair of successive motion shapes (Figure 2B; Martinez et al., 2018; Corn et al., 2021). Kinesis skew was then calculated by taking the natural logarithm of kinesis across the last five motion shapes divided by total craniofacial kinesis (Figure 2B; Supplemental Figure 1A), providing information about the distribution of movement across a strike. Finally, we extracted starting head shape for each strike – the sole morphological trait included in this study – which describes craniofacial form when in the closed mouth position at the start of a feeding sequence. Here, we performed a separate GPA to align the morphometric data extracted from the first frame of each feeding sequence (Supplemental Table 1). Total craniofacial kinesis, kinesis skew, and starting head shape values were computed for each of the 689 feeding strikes, averaged by specimen, and then by species (Supplemental Table 2). Additional information on data analyses can be found in the methods of the manuscript.  ### Code/Software, Data, & Supplemental Files * **RCode_EVL3-23-0152**: This markdown file includes all analyses reported in the manuscript and the supplemental material, excluding the MuSSCRat analyses which were run in RevBayes. There are notes and comments throughout the file to guide others in running the code and/or similar analyses. The code will import the following data files: * **PJTree_final_trim.tre**: Phylogenetic tree with 228 tips representing all study species.  * **SuppTable1_scaledLM.csv**: Data file that contains landmark morphometric data for each frame of the feeding sequence for all 689 strikes; all data is scaled but unaligned. X-Y positions are reported for each of the 18 landmarks (descriptions in Supplemental Figure 1). We include additional information on the family classification; pharyngeal jaw configuration; dataset that the morphometric data originates from (ARH – original data collected for this study; CMM – Martinez, Corn et al., 2024; KAC – Corn et al., 2021); and the video frame that each row of data originates from, where 1 represents the start of the feeding sequence. This entire data file matches the data within Supplemental Table 1.  * **FamFeed_data.csv**:** **Data file that contains family classification for fishes with a modified pharyngeal jaw.  * **sliding_semi.txt**: Data file used for the alignment of morphometric data.  * **phyl_KinData.csv**: Data file that contains all species average trait data including *Total Craniofacial Kinesis*, *Kinesis Skew*\ (log value), *Premaxillary Protrusion*, *Maxillary Rotation*, *Lower Jaw Rotation*, *Cranial Rotation*, *Hyoid Depression*, *Oral Jaw Gape*, *Starting Head Shape* (average landmark data from the first feeding sequence and the scores from the first 10 PC axes). We include additional information on the family and family classification; pharyngeal jaw configuration; dataset that the morphometric data originates from (ARH – original data collected for this study; CMM – Martinez, Corn et al., 2024; KAC – Corn et al., 2021); This data matches the data within Supplemental Table 2, excluding some other pieces of information regarding species proxies.  * **MuSSCRat output file**: Including *output_UCLN_KinesisFin_params.log*, *output_UCLN_KinSkewFin_params.log*, *output_UCLN_KinCompsFin_params.log*, and *output_UCLN_HeadFin_params.log*, these files include all of the posterior parameters from the MuSSCRat analyses run in RevBayes.  * **PJstateTrimmed_rb.nex**: Data file containing the pharyngeal jaw configuration for each species. The species order matches that in other .nex data files.  * **MuSSCRat Code Files**: Including *MuSSCRat_UCLN_script_KinesisTrim.txt*, *MuSSCRat_UCLN_script_KinSkewTrim.txt*, *MuSSCRat_UCLN_script_KinCompsTrim.txt*, and *MuSSCRat_UCLN_script_HeadshapeTrim.txt*, these script files outline how to run the MuSSCRat analyses RevBayes. These script files will import the following data files: * **PJTree_final_trim.tre**: Phylogenetic tree with 228 tips representing all study species.  * **Data files**: KinesisTrimmed_rb.nex, KinSkewTrimmed_rb.nex, KinCompsTrimmed_rb.nex, HeadshapePCTrimmed_rb.nex. All data matches that in **phyl_KinData.csv**, but has been placed in the proper format for MuSSCRat analysis.  * **PJstateTrimmed_rb.nex**: Data file containing the pharyngeal jaw configuration for each species. * **How to run MuSSCRat in RevBayes.docx**: Document with instructions and notes on how to run MuSSCRat analyses in RevBayes, analyze the quality of the output, and visualize the results. While MuSSCRat analyses are run in RevBayes, researchers will primarily use Tracer (software) and R to analyze the quality of the output and visualize the results.  * **Supplemental_Tables_1-9_final.xlsx**: Supplemental Tables 1 and 2 contain data used to complete analyses. Supplemental Tables 3-9 include results from the study. Definitions and descriptions of all traits can be found in Supplemental Figure 1. * **Supplemental Table 1**  contains landmark morphometric data for each frame of the feeding sequence for all 689 strikes; all data is scaled but unaligned. X-Y positions are reported for each of the 18 landmarks (descriptions in Supplemental Figure 1). We include additional information on the family classification; pharyngeal jaw configuration; dataset that the morphometric data originates from (ARH – original data collected for this study; CMM – Martinez, Corn et al., 2024; KAC – Corn et al., 2021); and the video frame that each row of data originates from, where 1 represents the start of the feeding sequence. * **Supplemental Table 2** contains average species data for all traits including total craniofacial kinesis; kinesis skew; six size-scaled motion components; and the aligned, morphometric data from the first frame of the feeding sequence representing interspecific head shape. Principal component (PC) scores for each species from the first 10 PC axes of interspecific head shape are also included as these data represent 95% of the total variation of head shape and were used to compute state-dependent rates of continuous character evolution in MuSSCRat. We include additional information regarding each species' accepted scientific name (Fishbase, July 2023) if it differs from the species name on the tree, the species used as proxy for taxa in the dataset that are not present in the phylogenetic tree (n = 30), family classification (Fishbase, July 2023), pharyngeal jaw configuration, dataset that the land- mark morphometric data originates from (ARH – original data collected for this study; CMM – Martinez, Corn et al., 2024; KAC – Corn et al., 2021), and the number of observations that comprise each species’ average trait values. ## Access information Data was derived from the following sources: * All derivations from source code are noted within the R markdown file