Study specimens and anatomical measurements

We examined 85 Neotropical cichlid species (Family Cichlidae) and 30 North American sunfish species (Family Centrarchidae) with an average of 5 specimens per species. Cichlidae, being one of seven families where modified pharyngeal jaws are an independently derived synapomorphy for all known species, represents fishes with MPJs. Centrarchids were used to represent fishes lacking the MPJ, as all species exhibit a generalized pharyngeal jaw system that lacks all traits defining an MPJ (Fig. 1b). We used published morphological data for 218 cichlid specimens (Burress et al. 2020; Supplemental Table 1) and obtained an additional 170 and 129 adult cichlid and centrarchid specimens, respectively, from museum collections, collection trips, and the pet trade (Supplemental Table 1) in an effort to maximize the morphological diversity and generic coverage in both families. Within our dataset, we have 50% and 100% sampling of genera for Neotropical cichlids and centrarchids, respectively. Though Neotropical and African cichlids both exhibit trophic morphologies and functions unique to their radiations, these two cichlid sister groups show extensive convergent evolution and shared ecomorphs (Burress 2014; Burress et al. 2017). Thus, we restricted our study to the continental radiation of Central and South American river and lake-dwelling cichlid species and we aimed to include species that would showcase the full breadth of ecomorphological variation exhibited by this spectacularly diverse clade (Burress 2014; Arbour et al. 2020). By excluding African cichlids, we also reduce the number of confounding factors that contribute to this group’s ecomorphological diversity, including lake effects (Seehausen 2006) and hybridization (Joyce et al. 2011; Meier et al. 2017, 2019). All specimens were fixed in formalin and stored in 70% ethanol prior to this study. Each specimen was then cleared and stained using procedures adapted from Dingerkus and Uhler (1977) and Taylor and Van Dyke (1985), and stained specimens were stored in 90% glycerin. Upper and lower pharyngeal jaws were dissected from each stained specimen and photographed. The left side of each stained specimen was also photographed in a relaxed posture where oral jaws were mostly protracted, according to supplemental procedures outlined in Roberts et al. 2021.

Using digital calipers and NIH ImageJ software (Schneider et al. 2012), we measured anatomical traits to characterize the size, shape, and functionality of the oral jaw and pharyngeal jaw systems. Traits were separated into oral and pharyngeal jaw systems based on their contributions to prey capture or processing functions, respectively (Alexander 1967; Lauder 1983; Wainwright 2005). We measured and computed 13 oral jaw traits including premaxillary dentigerous arm length, premaxillary ascending process length, gape, oral jaw protrusion, buccal cavity length, the horizontal and vertical positions of the lower jaw joint, the horizontal and vertical positions of the maxilla-nasal joint, lower oral jaw length, opening and closing mechanical advantage, and kinesis (Figs. 1c and 3b; Supplemental Tables 1 and 2). Kinesis, calculated according to procedures in Burress et al. 2020, describes the amount of shape change in the oral, four-bar linkage following a fixed rotation of the mandible (Martinez et al. 2018). Finally, we measured seven pharyngeal jaw traits including lower pharyngeal jaw length, width, depth, and tooth width (averaged across the three largest teeth), as well as upper pharyngeal jaw length, width, and depth (Figs. 1c and 3a; Supplemental Tables 1 and 2). Pharyngeal jaws from especially small species (e.g., Apistogramma and Enneacanthus) were measured with calipers under a dissecting microscope. All pharyngeal jaw measurements for the 218 cichlid specimens from Burress et al. 2020, excluding average tooth width, were retaken, as original trait data collected were not comparable to fishes with generalized pharyngeal jaws.

 Morphological measurements for each specimen were converted to shape ratios where each trait value was divided by the cubed root of (head length x head width x head depth). By size correcting each specimen using the geometric mean of these three main size dimensions, we account for the complexity of body size and our data capture shape variation due to allometry (Claude 2013; Klingenberg 2016; Price et al. 2019). Shape ratios were averaged across specimens, log transformed to achieve normal distributions, and used in all subsequent analyses in R v. 4.2.0 (R Core Team 2022). For select multivariate analyses, each species’ log shape ratios were converted to standard normal deviates (i.e., z-scores) using the scale function (Becker et al. 1988) to ensure that all 20 measurements were in commensurate units and scales (Huttegger and Mitteroecker 2011; Adams and Collyer 2019b). These standardized values were used for analyses where noted.

Statistical analyses

To conduct our analyses in a phylogenetic context, we used a time-calibrated molecular phylogeny of ray-finned fishes (Rabosky, Chang, Title, and Alfaro et al. 2018). This phylogeny, which includes a dense sampling of species in both focal groups, was trimmed to include fishes represented by molecular sequence data and then to match the species included our study. For species that were not present on the phylogeny, we used the closest related species within the same genus as proxy. Four substitutions were made, including Apistogramma hoignei for Apistogramma cacatuoides, Cichlasoma dimerus for Cichlasoma orientale, Pterophyllum leopoldi for Pterophyllum scalare, and Amatitlania nigrofasciata for Amatitlania siquia.

To visualize the multivariate morphological data in shape space, we performed Principal Component Analyses (PCA) on the correlation matrix of unstandardized, log shape ratios using the prcomp function (Mardia et al. 1979; Becker et al. 1988; Venables and Ripley 2002) in STATS v. 4.2.0 (R Core Team 2022). Analyses were run on all 20 morphological traits for cichlids and centrarchids separately, as well as in one combined analysis. We ran additional PCAs on the correlation matrices of 13 oral jaw traits and the seven pharyngeal jaw traits separately to visualize how cichlids and centrarchids fall in shape space based on each jaw unit. To test whether the cichlids and centrarchids have statistically different oral and pharyngeal jaw morphology, we performed a phylogenetic multivariate analysis of variance (phylogenetic MANOVA) on the z-scores of all 20 morphological traits over 10,000 iterations under a Brownian Motion model. This same method, implemented with the procD.pgls function (Adams 2014a; Adams and Collyer 2015, 2016, 2018; Collyer et al. 2015) in GEOMORPH v. 4.0.4 (Collyer and Adams 2018; Baken et al. 2021; Adams and Collyer 2022; Adams et al. 2023), was used to perform individual phylogenetic ANOVAs on each of the 20 traits in our dataset.

To quantify and compare the magnitude of interspecific variation among cichlid and centrarchid morphologies, we used morphol.disparity (Zelditch et al. 2012; Collyer and Adams 2021) in GEOMORPH. Using unstandardized, log shape ratios, we estimated overall disparity for all 20 morphological traits, two subsets of traits representing each jaw unit, and for each individual trait. We also estimated and contrasted the Brownian multivariate rate parameter for all 20 traits,  two subsets of traits representing each jaw unit, and for each individual trait using compare.evol.rates (Adams 2014b; Denton and Adams 2015; Adams and Collyer 2018, 2019a) in GEOMORPH.

We used the phylo.integration function (Adams and Felice 2014; Collyer et al. 2015; Adams and Collyer 2016, 2018, 2019a) in GEOMORPH to quantify the degree of morphological covariance between oral jaw and pharyngeal jaw traits in an evolutionary context over 10,000 iterations. As this was done separately on the z-scores from cichlids and centrarchids, we then compared the strength of between-jaw-unit integration between families with the compare.pls function (Collyer et al. 2015; Adams and Collyer 2016) in GEOMORPH. For each family, we used standard normal deviates to calculate the degree of evolutionary integration within each jaw system using the integration.Vrel function (Pavlicev et al. 2009; Conaway and Adams 2022) and compared the results using the compare.ZVrel (Conaway and Adams 2022) in GEOMORPH. Further, using standard normal deviates and the phylo.modularity and compare.CR functions (Adams and Felice 2014; Adams 2016; Adams and Collyer 2019a) in GEOMORPH, we quantified and contrasted the degree of evolutionary modularity in each family. Finally, we produced the evolutionary correlation matrices for each family using code modified from integration.Vrel and the cor function (Kendall 1938, 1945; Becker et al. 1988) in STATS. We then used the corrplot function (Murdoch and Chow 1996; Friendly 2002) in CORRPLOT v. 0.92 (Wei and Simko 2021) to visualize these matrices.

• Version 2 (7/21/23) -- Supplemental materials now include an R markdown file that outlines all analyses completed in R studios as well as two dataframes and a tree file to be uploaded when running analyses. This is in addition to the three supplemental files that were previously published. Please note that .csv files include species mean data for each trait. Also note that Supplemental Table 2 now includes information on which data columns in each .csv file correspond to which traits. Please see the README file for additional details about each published file.
• Version 1 (6/21/23) -- Supplemental files include trait data for each specimen, definitions for all measurements taken and calculated, and PCA loadings.