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.

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).

Morphological and functional disparity

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 et al., 1979; Becker et al., 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 et al., 2012; Collyer & Adams, 2021) in GEOMORPH to compute disparity – the sum of the diagonal elements of the variance-covariance matrix divided by the number of observations (Zelditch et al., 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 et al., 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.

Patterns of trait coevolution

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 et al., 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, Vrel) and converting them 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^th 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.

Rates of morphological and functional evolution

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 et al., 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 et al., 2018) to assess model convergence and code modified from REVGADGETS package v.1.1.0 (Tribble et al., 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 et al., 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).

Comparing rate estimates between maximum-likelihood and Bayesian frameworks

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 et al., 2006; Thomas et al., 2006; Beaulieu et al., 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 et al., 2021; Larouche et al., 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 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.

Comparing parameter estimates between full and subsampled datasets 

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.