Integration drives rapid phenotypic evolution in flatfishes.

Evolutionary innovations are scattered throughout the tree of life, and have allowed the organisms that possess them to occupy novel adaptive zones. While the impacts of these innovations are well documented, much less is known about how these innovations arise in the first place. Patterns of covariation among traits across macroevolutionary time can offer insights into the generation of innovation. However, to date, there is no consensus on the role that trait covariation plays in this process. The evolution of cranial asymmetry in flatfishes (Pleuronectiformes) from within Carangaria was a rapid evolutionary innovation that preceded the colonization of benthic aquatic habitats by this clade, and resulted in one of the most bizarre body plans observed among extant vertebrates. Here, we use three-dimensional geometric morphometrics and a phylogenetic comparative toolkit to reconstruct the evolution of skull shape in carangarians, and quantify patterns of integration and modularity across the skull. We find that the evolution of asymmetry in flatfishes was a rapid process, resulting in the colonization of novel trait space, that was aided by strong integration that coordinated shape changes across the skull. Our findings suggest that integration plays a major role in the evolution of innovation by synchronizing responses to selective pressures across the organism.

Evolutionary innovations are adaptations or exaptations that result in the colonization of novel regions of trait space, and shift the dynamic of ecosystem interactions in their respective environments. These innovations may involve a dramatic restructuring of an ancestral body plan, allow organisms novel access to ecological resources, and in (the case of key innovations) spur increases in lineage diversification (1–6).

Traditionally, much of the study of evolutionary innovation has focused on extrinsic drivers, such as ecological change or environmental stimuli, which have frequently been viewed as the principal triggers for evolutionary novelty (6–9). Recently, focus has shifted toward the effects of intrinsic (e.g., developmental) processes that may structure patterns of trait diversification (10–12). Among these processes, integration and modularity have emerged as important sources of insight as a growing consensus has found that patterns of trait covariation can constrain and facilitate responses to selection and strongly influence patterns of trait diversification at both contemporary and macroevolutionary timescales (13–20).

Modularity refers to the pattern whereby traits form complexes (i.e., modules) that exhibit a high degree of covariation within themselves, whereas the strength of covariation is far lower between these trait complexes (21, 22). This can result in the compartmentalization of trait complexes across an organism and has been hypothesized to allow modules to respond semiautonomously to different selective pressures (23–27). On macroevolutionary timescales, this modularization of different traits has been shown to result in mosaic patterns of evolution across organisms and is thought to facilitate morphological diversification as different traits are able to fine-tune responses to different selective pressures (27–29). Conversely, integration refers to a pattern whereby different traits exhibit a high degree of covariation (21, 30). Patterns of integration may be the result of pleiotropy or functional coupling (28, 30–33). There is less of a consensus on the macroevolutionary implications of phenotypic integration. Strong integration between traits has traditionally been hypothesized to constrain patterns of trait diversification along specific directions such that traits may exhibit more variation in some directions than others (30, 31). Furthermore, theoretical studies have found that evolutionary flexibility is negatively correlated with the magnitude of integration (27, 28, 34). However, more recent studies have found that the relationship between evolutionary integration and trait diversification is more complex than originally thought. A recent simulation study (28) found that trait integration can actually promote sizeable responses to selection, along lines that are parallel with the direction of trait covariation. This work also indicates that while a modular system has the freedom to explore a wider range of morphospace, this system is less likely to evolve maximally disparate phenotypes as compared to an integrated one that funnels variation along a more restricted trajectory. Indeed, several empirical studies have found rapid rates of shape evolution and high degrees of morphological disparity in clades that also exhibit tight evolutionary integration across the traits of interest (29, 35, 36). These studies suggest that evolutionary integration may play an important role in the evolution of innovation as it allows for clades to rapidly explore novel regions of morphospace.

The majority of our understanding of patterns of vertebrate trait diversification come from tetrapod systems, and primarily studies involving birds and mammals (14, 37–40). While these studies have been deeply informative, they help to explain less than half of the vertebrate story of diversification as these studies frequently exclude ray-finned fishes (Actinopterygii) (24). Ray-finned fishes comprise over half of the species diversity in vertebrates. They also exhibit a diverse array of evolutionary innovations that have allowed them to become arguably the most successful vertebrate radiation on the planet.

Flatfishes in particular represent a diverse clade of bottom-dwelling, teleost fishes that possess a striking evolutionary innovation: a degree of cranial asymmetry exceeding that of any other vertebrate lineage (41). Most flatfishes are completely blind on one side of their body, and instead feature both eyes on the same side of the head. This perplexing flatfish phenotype is achieved during the early developmental stages where one eye of a symmetrical larva gradually begins to migrate to the other side of its body, rendering one side “eyed” and the other “blind” (42). Developmental studies indicate that this orbital migration is driven by thyroid hormone expression and is paired with changes in swimming behavior (43), as well as asymmetrical visceral organ rearrangement (44). These developmental patterns recapitulate paleontological trends, with fossils indicating that flatfishes gradually became more asymmetrical from symmetrical ancestors (45). In addition to their cranial asymmetry, flatfishes possess a suite of adaptations that further allow them to exploit their benthic habitats; most nobly among these is their derived form of locomotion that combines strong, whole-body undulations (for fast escapes or burial) with finer-scale, dexterous undulations of their dorsal and anal fins (for slower cruising behaviors) (46–48). The evolution of this derived form of locomotion is coupled with additional changes in the roofing bones of the skull as they evolved projections to support the extended dorsal fin along the supraoccipital crest and the frontal bones (49, 50).

The flatfish body plan is thought to have rapidly arisen shortly after the Cretaceous–Paleogene boundary (66 Mya) (49). Comparative analyses show rapid bursts of body-shape diversification in flatfishes and their close relatives within the broader clade Carangaria, a diverse radiation including disparate lineages like billfishes and remoras in addition to flatfishes (49, 51, 52) (Fig. 1). This episode of phenotypic innovation is hypothesized to reflect the filling of newly available ecological roles in the early Cenozoic, matching patterns reported for other groups.

Fig. 1.

Phylogeny of 102 carangarian species included in the analyses of skull shape evolution. Insets depict representative skull shapes for each clade. Phylogeny based on Ribeiro et al. (52).

The striking and rapid evolutionary dynamics within Carangaria make this clade a tantalizing target for investigating the roles that integration and modularity play in trait diversification and, particularly, evolutionary innovation. Here we use three-dimensional geometric morphometrics and a cutting-edge phylogenetic comparative toolkit to study the evolution of the neurocranium and the evolution of cranial asymmetry across 102 carangarian species. We quantify shifts in the rate of skull shape evolution between species and across the neurocranium as a whole, while also quantifying patterns of integration and modularity between flatfishes and their relatives to test for the effect of integration and modularity on the evolution of innovation. We hypothesize that flatfishes underwent a rapid shift in their rates of skull shape evolution as a result of their orbital migration. We additionally hypothesize that flatfishes will exhibit higher levels of integration compared to their carangarian relatives as a result of their asymmetrical larval metamorphosis that involves coordinated changes across the body, and adaptations associated with their derived locomotory mode.

Results

Skull Shape Evolution in Carangaria.

Carangarians exhibit a broad diversity of skull shapes (Fig. 2). An analysis of phylogenetic signal finds significant, but weak phylogenetic structure in neurocranium shape across this clade (P = 0.001; Blomberg’s K = 0.178). The phylogenetic signal is most apparent along the primary axis of shape variation (principal component 1, PC1), which corresponds overwhelmingly to variation in the shape and position of the right (migrating) orbit. Interestingly, aside from the shape divergence toward the root of the phylogeny, there appears to be rampant convergence among the flatfishes, particularly among species with elongated neurocrania (e.g., Pleuronectes quadrituberculatus). It’s possible that this convergence is driving the low phylogenetic signal observed in this clade. The results of the phylogenetically aligned component analysis reveal a similar pattern of dispersion (). In both analyses, a strong asymmetrical pattern of shape variation can be observed across the neurocranium with the highest amount of variance concentrated at the right, frontal-sphenotic margin (Fig. 2). Using the phylogenetically aligned components analysis, we found that this variation in orbital position accounts for over 93% of the total shape variation for carangarians as a whole, suggesting that the early-diverging flatfish asymmetry is the primary contributor to phylogenetic signal in the shape data. In our phylomorphospace analysis, along the PC1 axis, flatfishes occupy a distinct region of morphospace with respect to their “nonflatfish” relatives, with the Alaska plaice (P. quadrituberculatus) exhibiting the lowest score and the common remora (Remora remora) exhibiting the highest score. The second PC axis (PC2) corresponds to variation in the skull depth and relative skull length (i.e., heterocephaly) with species like the roughscale tonguesole (Cynoglossus lida) and the moonfish (Mene maculata) exhibiting the lowest scores while more dorsoventrally compressed species, like the Petrale sole (Eopsetta jordani) and the Pacific barracuda (Sphyraena argentea), exhibit higher scores.

Fig. 2.

Skull shape diversity and evolution across 102 carangarian species. PC analysis showing the primary axes and regions of variation in neurocranium shape across Carangaria (A–F). (G) Phylomorphospace analysis showing the diversity of neurocranial shapes across carangarian fishes. Insets depict representative extremes for each PC axis.

Tempo and Mode of Skull Shape Evolution.

Carangarians exhibit variable rates of skull shape evolution. We recover strong model support for a “variable-rate” model of trait evolution as opposed to a “single-rate” model (). A disparity through time analysis indicates a rapid decline in subclade disparity that exceeds the Brownian motion expectation early in carangarian history (though on average, does not differ significantly across the entire history of the clade; mean disparity index = −0.012, P = 0.36). This rapid decline in subclade disparity is consistent with adaptively radiating lineages and suggests that skull shapes partitioned early in evolutionary history. After the initial rapid decline, subclade disparity remains fairly consistent over time, with a recent rapid spike less than 5 Mya coinciding with rapid rates of shape evolution occurring primarily within flatfishes specifically: Cynoglossus, Gymnachirus, and Lepidopsetta (Fig. 3). We found additional evidence for rapid shape divergences at the base of Carangaria using a BayesTraits analysis, which recovers a strong rate shift at the base of Carangaria. Additionally, we found that flatfishes underwent a rapid burst in the rate of skull shape diversification that coincided with the evolution of cranial asymmetry. Within flatfishes, rates of shape evolution generally decline with subsequent shifts in achirids and cynoglossids, coinciding with variation related to orbital position and relative skull length. Outside of flatfishes, we found the largest rate shifts among remoras (Echeneidae), bill-fishes (Istiophoridae), and lookdown jacks (Carangidae). These rate shifts coincide primarily with heterocephaly (variation between elongated and foreshortened skull shapes), dorsoventral compression, and the height of the supraoccipital crest. Comparisons of mean evolutionary rates between flatfishes and “nonflatfishes” found significantly higher rates of skull shape evolution among flatfishes (P < 0.001; rate ratio = 2.83). These results suggest that carangarians rapidly morphologically diversified early in their evolutionary history and that flatfishes continued to rapidly diversify in smaller subclades after initially evolving cranial asymmetry.

Fig. 3.

Tempo and mode of skull shape diversification in Carangaria. Results from BayesTraits and disparity though time analyses showing shifts in the rate of skull shape evolution across the carangarian phylogeny. Dashed line on DTT indicate Brownian motion expectation and shaded region indicates 95% confidence interval.

Evolutionary Modularity.

We found strong patterns of modularity across the carangarian neurocranium. We recovered strong model support for the six-module “bone” hypothesis encompassing the basicranium, ethmoid, orbit, parasphenoid, and supraoccipital regions (Table 1). Support for this model was retained after rerunning the analysis on a smaller subsample of landmarks (). Support for the bone hypothesis was also recovered using the covariance ratio (CR) coefficient (CR = 0.943; P < 0.001). We subsequently divided the Carangaria into two subclades consisting of flatfishes (Pleuronectiformes) and nonflatfishes (Carangiformes, Istiophoriformes, Centropomidae, Leptobramidae, Menidae, Polynemidae, and Toxotidae) and reanalyzed patterns of phylogenetic integration and modularity. We found that flatfishes exhibit a strong signal of integration (partial-least squares [pls] = 0.968) while still exhibiting a weak but significant degree of modularity (CR = 0.96; P < 0.001). Meanwhile, we found that nonflatfishes exhibit significantly (P < 0.001) lower levels of integration (pls = 0.730) and significantly (P < 0.001) higher levels of modularity (CR = 0.782; P < 0.001) across the six regions of the neurocranium (Fig. 4 and ).

Table 1.

Results for the evaluation of modularity hypotheses (using phylo.EMMLI) for 102 carangarian species

Model	MaxL	K	n	AICc	dAICc	Model_L	Post_Pob
No.modules. default	−29182.5	2	17956	58369.09	22557.38	0	0
Bone modules
Same within-module p + same between-module p	−19009.5	3	17956	38024.93	2213.216	0	0
Separate within-module p + same between-module p	−18374.7	8	17956	36765.47	953.755	7.85E-208	7.85E-208
Same within module p + separate between module p	−18518.6	17	17956	37071.16	1259.446	3.27E-274	3.27E-274
Separate within module p + separate between module p	−17883.8	22	17956	35811.72	0	1	1
Functional modules
Same within-module p + same between-module p	−25411.4	3	17956	50828.89	15017.18	0	0
Separate within-module p + same between-module p	−22287.6	5	17956	44585.19	8773.472	0	0
Same within module p + separate between module p	−25274.9	5	17956	50559.88	14748.17	0	0
Separate within module p + separate between module p	−22151.1	7	17956	44316.18	8504.465	0	0
Developmental modules
Same within-module p + same between-module p	−22261.3	3	17956	44528.56	8716.839	0	0
Separate within-module p + same between-module p	−20675.5	6	17956	41363.03	5551.312	0	0
Same within module p + separate between module p	−22059.1	8	17956	44134.13	8322.416	0	0
Separate within module p + separate between module p	−20473.3	11	17956	40968.61	5156.892	0	0

Bold text indicates optimal model.

Fig. 4.

Phylogenetic integration across six modules of the carangarian skull (A). Violin plots showing results from a phylogenetic two-block partial least-squares for flatfishes and nonflatfishes showing significant differences in mean partial least squares values for each group (B). Network graphs showing the magnitude of integration between modules for flatfishes and nonflatfishes (C and D).

Mosaic Evolution.

Among carangarians, we found that the fastest rates of shape evolution are concentrated in the ethmoid, orbital, and ventral basicranial regions of the neurocranium (all regions that undergo dramatic changes during flatfish metamorphosis), while the parasphenoid and parietal regions exhibit slower rates of shape evolution (Fig. 5 and ). Rate analyses (using σ2 mult) also recovered significant rate differences between modules (). Mosaic patterns of evolution appear to be driven primarily by the slower-than-average rates of shape evolution in the anterior region of the parasphenoid, which differs significantly from every other region. Excluding the parasphenoid, no other significant rate differences were recovered between modules.

Fig. 5.

Mosaic evolution across the carangarian skull. Per landmark rates of shape evolution for 102 carangarian species (A–C). Ridgeline plot comparing per landmark rates across the six modules of the neurocranium (D).

Discussion

Evolutionary innovations allow the organisms that possess them to colonize novel regions of trait space, and grant them access to new and different adaptive zones. In this study, we examined a major evolutionary innovation in flatfish cranial asymmetry. We found that cranial asymmetry evolved within the Carangaria shortly after the K/Pg boundary and resulted in a rapid shift in rates of skull shape evolution at the base of Pleuronectiformes. This innovation allowed flatfishes to rapidly colonize a completely novel region of morphospace as they invaded benthic habitats. Among carangarians, we recovered asymmetrical patterns of variation across the skull driven by the orbital migration that flatfishes experienced at the macroevolutionary timescale. This degree of cranial asymmetry is not seen anywhere else in the vertebrate tree of life, even among other teleost fishes that inhabit the benthos (but see ref. 53). Typically, benthic fishes exhibit dorsoventrally flattened bodies to allow for partial burial in substrates (e.g., lizard fishes) (51, 54). Dorsoventrally compressed phenotypes are even present within carangarians (e.g., remoras), making the cranial asymmetry of flatfishes even more puzzling since some of their relatives were clearly able to evolve flattened skulls in a more typical fashion.

The role that integration and modularity play in the evolution of innovation is a hotly debated topic (23, 55, 56). Traditional hypotheses suggest that phenotypic diversification and novelty can be facilitated by the modularization of traits to allow them to behave semiautonomously over macroevolutionary timescales (23, 27, 56). In contrast to modularity, integration between traits across macroevolutionary timescales is thought to constrain patterns of trait diversification (29, 32, 57, 58; but see ref. 28). Here we found that the evolution of cranial asymmetry in flatfishes was an integrated process, with coordinated patterns of shape change across the neurocranium. This result may seem counterintuitive given the superficially modular nature of the orbital migration typified in flatfishes. However, the benthopelagic transition that flatfishes experience at the developmental level reveals a more integrated process across the skull and body that includes changes in bilateral pigmentation, hindbrain morphology, inner-ear morphology, swimming behavior, relative muscle masses between eyed and blind sides, and the relative placements and sizes of internal organs (59–61). Additional evidence for the integrated nature of the benthopelagic transition can be found in flatfish locomotion, which utilizes both the dorsal, and anal fins as paired appendages for benthic “walking” (46). This walking behavior is accomplished by the more anterior placement of the dorsal fin (relative to nonflatfishes), which interdigitates with the supraoccipital crest and other dermal-roofing bones in many species (e.g., Achirus). Similarly, widespread burial behavior exhibited across most flatfish species also involves coordinated sinusoidal undulations across the entire body and fins to cover the animal with sediment. These behavioral novelties almost certainly required integrated modifications to both axial and cranial morphologies.

While our analyses recover strong patterns of evolutionary integration across extant flatfishes, the fossil record suggests that trait integration may have gradually increased over evolutionary time. Studies of stem flatfishes indicate that the evolution of asymmetry was a gradual, and stepwise process [contrary to the saltatory hypothesis proposed by Goldschmidt (45, 62, 63)]. While stem flatfishes exhibited some characteristic patterns of cranial asymmetry observed in extant flatfishes (primarily asymmetrical orbits), they retained a number of primitive and generalized morphologies as well [including a “generalized perciform caudal skeleton” (45), which suggests that cranial asymmetry may have predated the adaptations associated with flatfish locomotion (however, more complete fossils will be needed to confirm this). Many of these characters can still be observed in Psettodes, an early-branching, phenotypically conserved flatfish. Members of this genus exhibit incomplete orbital migration, pigmentation on both sides of the body, and left-handed and right-handed morphs that are present at equal frequencies in the population (considered a precursor to directional asymmetry) (42). This fluctuating asymmetry, may be indicative of a less-canalized developmental program, and could possibly suggest that the asymmetry in flatfishes arose originally from developmental instability in the cranial growth fields that may have subsequently been reorganized in later clades (64).

Taken together, our findings suggest that integration played a key role in coordinating changes across different traits during the benthopelagic transition in flatfishes. We further hypothesize that integration likely plays a key role in the evolution of innovation by distributing variation across multiple traits and allowing organisms to mount a coordinated response to selective pressures.

Materials and Methods

Morphological Sampling and General Shape Analyses.

Skull shape was quantified across 123 specimens representing 102 carangarian species (61 flatfishes, 41 nonflatfishes, ∼10% of total clade diversity, and 24 of the 27 carangarian families) (). Specimens were micro-CT scanned at the University of Michigan and the University of Washington using a Nikon XT H 225 ST micro-CT scanner and a Bruker Skyscan 1173 micro-CT scanner, respectively, in conjunction with oVert and ScanAllFishes scanning initiatives. Additional scans were also downloaded from http://morphosource.org. Scans were segmented using Amira and the neurocranium was isolated from each specimen. Segmented neurocrania were converted to three-dimensional meshes and exported as .ply files.

Flatfishes are known to exhibit both sinistral and dextral forms. These drastic differences in asymmetry have the ability to confound geometric morphometric analyses. To account for these differences in asymmetry, we inverted the meshes for all sinistral (left-eyed) specimens in MeshLab turning them into dextral specimens. The inversion of the meshes effectively constrained the asymmetrical variation present in flatfish neurocrania to the right side for all specimens. While useful for standardizing our specimens, it bears noting that our approach may systematically underestimate shape variation and disparity in the flatfish skull.

After standardization, meshes were digitized with 134 (34 fixed; 100 sliding) three-dimensional landmarks using the Checkpoint software (). Bone nomenclature follows Hilton et al. (65). Landmarks were not treated as bilaterally symmetrical due to the high degree of asymmetry present in flatfishes, which comprise over half of the shape data. Landmark coordinates were then imported into R and analyzed using the geomorph package (66). For specimens with more than one representative per species, landmark averages were calculated. Specimen coordinates were subjected to a generalized procrustes analysis to remove the effect of orientation and scaling between specimens (67). Semilandmarks were slid along their tangent directions while minimizing bending energy (68).

Allometric scaling can play an important role in structuring patterns of shape change (69, 70). To quantify the effect of allometric scaling on our skull shape data, we performed a phylogenetic generalized least-squares regression of shape vs. log(centroid-size) using the proc.pgls function in geomorph (71, 72). We found a small (r2 = 0.07) but significant (P = 0.0015) effect of size on shape (). Due to its relatively small effect on shape variation and its general importance as relevant biological information, we elected not to remove the effect of allometric scaling from our shape data following the approach of Evans et al. (16) and Felice and Goswami (37).

Phylogenetic Reconstruction.

To study the evolution of shape change across Carangaria, we used the phylogeny of Ribeiro et al. (52). This time-calibrated phylogeny consists of 508 tips and comprises a phylogenetic backbone of Rabosky et al. (73), with additional carangarian taxa grafted onto it. We pruned the Ribeiro et al. (52) phylogeny in ape (74) to only include the taxa present in our study. The resulting pruned tree consisted of 102 carangarian species.

Phylogenetic Visualization.

To visualize major axes of shape variation and evolution, we utilized a phylomorphospace approach using our first two PCs of shape variation (accounting for 56% of total shape variance) and projected the Ribeiro et al. (52) phylogeny onto them (75). In order to compare phylogenetic trends in our data, we coupled our PC analysis with a phylogenetically aligned components analysis (76), which uses an ordinary least-squares approach to align the trait data to the axis of greatest phylogenetic signal and projected the first two components onto the aforementioned phylogeny. This approach allowed us to differentiate between stronger and weaker phylogenetic signals localized to a particular proportion of the shape axes. We additionally calculated phylogenetic signal in our shape data using the multivariate implementation of Blomberg’s K (kmult), to test for phylogenetic structure in our skull shape data (77).

Quantifying Clade Rates of Shape Evolution.

Rates of shape evolution across carangarian subclades were quantified using the variable-rates model implementation in the BayesTraitsV3 program (78). This Bayesian method uses a reversible-jump Markov chain Monte Carlo chain approach to estimate the probability of rate shifts in continuous trait data across a phylogeny. This method has the ability to detect clade and species-specific rate shifts in trait data. There are known difficulties with evolutionary model-fitting using data with high dimensionality (79). To reduce the dimensionality of our data, we used the first 20 PCs, which accounted for 95% 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 million generations, discarding the first 1 million as burn-in. The chain was sampled every 2 million generations after convergence using a stepping-stone sampler and convergence was evaluated before the first 1 million generations were removed as burn-in. 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 (80) (). We evaluated two models of trait evolution: a single-rate model that assumes a single Brownian motion rate of trait evolution across the phylogeny, and 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 (78). Model comparisons were performed by calculating Bayes factors from the marginal likelihoods of the single-rate and variable-rates models. 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. To summarize these rate shifts, we built a consensus tree out of the 100 trees produced using the BayesTrees program.

We also quantified the rate of shape evolution between flatfishes and nonflatfishes using the compare.evol.rates function in geomorph (81). Significance was assessed using the phylogenetic simulation approach [see Adams and Collyer (79) for more information].

Disparity through Time.

A disparity through time (DTT) analysis was used to model changes in subclade shape disparity for carangarians under a Brownian motion model (using our entire shape dataset) (82). This method calculates changes in relative subclade disparity through time across every node in the phylogeny. We then compared the observed disparity of carangarians to a Brownian motion model that was simulated 1,000 times across the phylogeny. We used the observed and simulated disparities to calculate a morphological disparity index (MDI), which quantifies the difference in relative subclade disparity under a Brownian motion expectation. Negative MDI values are generally indicative of adaptively radiating clades (82, 83). The DTT analysis was performed in the R-package geiger (84).

Evolutionary modularity across different regions of the carangarian neurocranium was assessed using the phyloEMMLi approach; this method uses maximum-likelihood to compare different hypotheses of modularity across a shape dataset while accounting for phylogenetic nonindependence of the trait data (85). We fit three models of modularity (and a “no modularity” model) to our data (). The first bone model consisted of six modules (vomer, parasphenoid, parietals, “basicranium,” supraoccipital, and “orbit”) roughly corresponding to individual bones; with the exception of the orbital and basicranial regions, which are composed of multiple smaller bones that were not densely sampled enough in the landmarking scheme to constitute separate module designation. The second model consisted of three “functional” modules (ethmoid, optic, and basicranium) corresponding to the location of the individual sensory capsules in the neurocranium (86). The otic region was not given an individual module designation apart from the basicranium due to the frequent lack of ossification in the area, which limited adequate landmark sampling. The third model consisted of four “developmental” modules (face, orbit, basicranium, and skull roof) following Evans et al. (16, 87, 88). The robustness of the EMMLI results were evaluated with a random subsampling down to 90% of the full landmark dataset, and EMMLi ran iteratively 100 times with an Akaike Information Criterion cutoff of 500, using the “subSampleEMMLi” function from the EMMLi v.0.0.3 R package. Mean results from the 100 subsamples were then compared to the analysis of the full dataset following the approach of Bon et al. (89).

Evolutionary modularity was also assessed using the CR coefficient (90), which quantifies the degree of modularity between shape variables while taking phylogeny into account. Significance is determined by comparing the observed CR coefficient with a randomly generated distribution.

Differences in the strength of modularity was assessed between flatfishes and nonflatfishes using the compare.CR function in geomorph (91). This function compares the strength of modularity between shape configurations using the CR effect sizes from each dataset.

Evolutionary Integration.

Evolutionary integration was assessed using the phylo.integration function in geomorph (92). This method uses a two-block partial least-squares analysis to quantify the degree of shape covariation between hypothesized modules while considering phylogenetic nonindependence. This analysis was run separately for flatfishes and nonflatfishes. Our best-supported modular configuration exceeded three modules, so we calculated an r coefficient value for each individual pairwise comparison using the phylo.integration function and tested for differences between the mean partial least-squares r coefficients for each group using a two-sample t test. The r coefficient is sensitive to differences in sample size, and differences in the number of variables between datasets (93). To account for this, we randomly removed 20 flatfish species prior to running the integration analyses to standardize the sample sizes between both groups.

Quantifying Mosaic Evolution.

Mosaic patterns of shape evolution across the neurocranium were assessed and quantified by calculating the per landmark rate of evolution for every landmark in the shape configuration using the Hotdots package in R (56). Landmark rates were then subdivided by module, and rate differences between modules were assessed using a Tukey’s honest significant difference to compare pairwise rate differences between modules while adjusting for multiple comparisons (94). Module rates of evolution were also quantified using the compare.multi.evol.rates function in geomorph which calculates rates of shape evolution from a priori defined modules [see Denton and Adams (81) for more information].