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.

References

• Adams, D. C. 2014a. A method for assessing phylogenetic least squares models for shape and other high-dimensional multivariate data. Evolution (N Y) 68:2675–2688.
• Adams, D. C. 2016. Evaluating modularity in morphometric data: Challenges with the RV coefficient and a new test measure. Methods Ecol Evol 7:565–572.
• Adams, D. C. 2014b. Quantifying and comparing phylogenetic evolutionary rates for shape and other high-dimensional phenotypic data. Syst Biol 63:166–177.
• Adams, D. C., and M. L. Collyer. 2019a. Comparing the strength of modular signal, and evaluating alternative modular hypotheses, using covariance ratio effect sizes with morphometric data. Evolution (N Y) 73:2352–2367.
• Adams, D. C., and M. L. Collyer. 2022. Consilience of methods for phylogenetic analysis of variance. Evolution (N Y) 76:1406–1419. Society for the Study of Evolution.
• Adams, D. C., and M. L. Collyer. 2018. Multivariate phylogenetic comparative methods: Evaluations, comparisons, and recommendations. Syst Biol 67:14–31.
• Adams, D. C., and M. L. Collyer. 2016. On the comparison of the strength of morphological integration across morphometric datasets. Evolution (N Y) 70:2623–2631.
• Adams, D. C., and M. L. Collyer. 2015. Permutation tests for phylogenetic comparative analyses of high-dimensional shape data: What you shuffle matters. Evolution (N Y) 69:823–829.
• Adams, D. C., and M. L. Collyer. 2019b. Phylogenetic comparative methods and the evolution of multivariate phenotypes. Annu Rev Ecol Evol Syst 50:405–425.
• Adams, D. C., M. L. Collyer, A. Kaliontzopoulou, and E. K. Balken. 2023. Geomorph: Software for geometric morphometric analyses.
• Adams, D. C., and R. N. Felice. 2014. Assessing trait covariation and morphological integration on phylogenies using evolutionary covariance matrices. PLoS One 9:e94335.
• Alexander, R. McN. 1967. The functions and mechanisms of the protrusible upper jaws of some acanthopterygian fish. J. Zool 151:43–64.
• Arbour, J. H., C. G. Montaña, K. O. Winemiller, A. A. Pease, M. Soria-Barreto, J. L. Cochran-Biederman, and H. López-Fernández. 2020. Macroevolutionary analyses indicate that repeated adaptive shifts towards predatory diets affect functional diversity in Neotropical cichlids. Biological Journal of the Linnean Society 129:844–861.
• Baken, E. K., M. L. Collyer, A. Kaliontzopoulou, and D. C. Adams. 2021. geomorph v4.0 and gmShiny: Enhanced analytics and a new graphical interface for a comprehensive morphometric experience. Methods Ecol Evol 12:2355–2363.
• Becker, R. A., J. M. Chambers, and A. R. Wilks. 1988. The New s Language: A programming environment for data analysis and graphics. Wadsworth & Brooks/Cole.
• Burress, E. D. 2014. Cichlid fishes as models of ecological diversification: patterns, mechanisms, and consequences. Hydrobiologia 748:7–27. Kluwer Academic Publishers.
• Burress, E. D., C. M. Martinez, and P. C. Wainwright. 2020. Decoupled jaws promote trophic diversity in cichlid fishes. Evolution (N Y) 74:950–961.
• Burress, E. D., L. Piálek, J. R. Casciotta, A. Almirón, M. Tan, J. W. Armbruster, and O. Říčan. 2017. Island- and lake-like parallel adaptive radiations replicated in rivers. Proceedings of the Royal Society B 285:20171762. Royal Society Publishing.
• Claude, M. 2013. Log-shape ratios, Procrustes superimposition, elliptic Fourier analysis: Three worked examples in R. Hystrix 24:94–102.
• Collyer, M. L., and D. C. Adams. 2021. Phylogenetically aligned component analysis. Methods Ecol Evol 12:359–372.
• Collyer, M. L., and D. C. Adams. 2018. RRPP: An R package for fitting linear models to high-dimensional data using residual randomization. Methods Ecol Evol 9:1772–1779.
• Collyer, M. L., D. J. Sekora, and D. C. Adams. 2015. A method for analysis of phenotypic change for phenotypes described by high-dimensional data. Heredity (Edinb) 115:357–365. Nature Publishing Group.
• Conaway, M. A., and D. C. Adams. 2022. An effect size for comparing the strength of morphological integration across studies. Evolution (N Y) 76:2244–2259.
• Denton, J. S. S., and D. C. Adams. 2015. A new phylogenetic test for comparing multiple high-dimensional evolutionary rates suggests interplay of evolutionary rates and modularity in lanternfishes (Myctophiformes; Myctophidae). Evolution (N Y) 69:2425–2440.
• Dingerkus, G., and L. D. Uhler. 1977. Enzyme clearing of alcian blue stained whole small vertebrates for demonstration of cartilage. Stain Technol 52:229–232.
• Friendly, M. 2002. Corrgrams: Exploratory displays for correlation matrices. Am Stat 56:316–324.
• Huttegger, S. M., and P. Mitteroecker. 2011. Invariance and meaningfulness in phenotype spaces. Evol Biol 38:335–351.
• Joyce, D. A., D. H. Lunt, M. J. Genner, G. F. Turner, R. Bills, and O. Seehausen. 2011. Repeated colonization and hybridization in Lake Malawi cichlids. Current Biology 21:R108–R109. Elsevier.
• Kendall, M. G. 1938. A new measure of rank correlation. Biometrika 30:81–93.
• Kendall, M. G. 1945. The treatment of ties in ranking problems. Biometrika 33:239–251.
• Klingenberg, C. P. 2016. Size, shape, and form: concepts of allometry in geometric morphometrics. Dev Genes Evol 226:113–137. Springer Verlag.
• Lauder, G. V. 1983. Functional design and evolution of the pharyngeal jaw apparatus in euteleostean fishes. Zool J Linn Soc 77:1–38.
• Mardia, K. V., J. T. Kent, and J. M. Bibby. 1979. Multivariate Analysis. Academic Press.
• Martinez, C. M., M. D. McGee, S. R. Borstein, and P. C. Wainwright. 2018. Feeding ecology underlies the evolution of cichlid jaw mobility. Evolution (N Y) 72:1645–1655.
• Meier, J. I., D. A. Marques, S. Mwaiko, C. E. Wagner, L. Excoffier, and O. Seehausen. 2017. Ancient hybridization fuels rapid cichlid fish adaptive radiations. Nat Commun 8. Nature Publishing Group.
• Meier, J. I., R. B. Stelkens, D. A. Joyce, S. Mwaiko, N. Phiri, U. K. Schliewen, O. M. Selz, C. E. Wagner, C. Katongo, and O. Seehausen. 2019. The coincidence of ecological opportunity with hybridization explains rapid adaptive radiation in Lake Mweru cichlid fishes. Nat Commun 10. Nature Research.
• Murdoch, D. J., and E. D. Chow. 1996. A graphical display of large correlation matrices. Am Stat 50:178–180.
• Pavlicev, M., J. M. Cheverud, and G. P. Wagner. 2009. Measuring morphological integration using eigenvalue variance. Evol Biol 36:157–170.
• Price, S. A., S. T. Friedman, K. A. Corn, C. M. Martinez, O. Larouche, and P. C. Wainwright. 2019. Building a Body Shape Morphospace of Teleostean Fishes. Integr Comp Biol 59:716–730. Oxford University Press.
• Rabosky, D. L., J. Chang, P. O. Title, P. F. Cowman, L. Sallan, M. Friedman, K. Kaschner, C. Garilao, T. J. Near, M. Coll, and M. E. Alfaro. 2018. An inverse latitudinal gradient in speciation rate for marine fishes. Nature 559:392–395. Springer US.
• Roberts, A. S., J. R. Hodge, P. Chakrabarty, and P. C. Wainwright. 2021. Anatomical basis of diverse jaw protrusion directionality in ponyfishes (Family Leiognathidae). J Morphol 282:427–437.
• Schneider, C. A., W. S. Rasband, and K. W. Eliceiri. 2012. NIH Image to ImageJ: 25 years of image analysis. Nat Methods 9:671–675. Nature Publishing Group.
• Seehausen, O. 2006. African cichlid fish: A model system in adaptive radiation research. Proceedings of the Royal Society B 273:1987–1998.
• Taylor, W. R., and G. C. Van Dyke. 1985. Revised procedure for staining and clearing small fishes and other vertebrates for bone and cartilage study. Cybium 9:107–121.
• Team, R. C. 2022. R: A language and environment for statistical computing.
• Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics with S. Springer-Verlag.
• Wainwright, P. C. 2005. Functional morphology of the pharyngeal jaw apparatus. Pp. 77–101 in Fish Physiology.
• Wei, T., and V. Simko. 2021. R package “corrplot”: Visualization of a Correlation Matrix.
• Zelditch, M. L., D. L. Swiderski, H. D. Sheets, and W. L. Fink. 2012. Geometric morphometrics for biologists: A primer. Academic Press.

• Version 3 (12/23/24) -- Updated trait data (ComboFin2_new.csv, Supplemental_Table_1_data.xlsx, SR_data_mean.csv)
• 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.