A rich diversity of opercle bone shape among teleost fishes

Abstract

The opercle is a prominent craniofacial bone supporting the gill cover in all bony fish and has been the subject of morphological, developmental, and genetic investigation. We surveyed the shapes of this bone among 110 families spanning the teleost tree and examined its pattern of occupancy in a principal component-based morphospace. Contrasting with expectations from the literature that suggest the local morphospace would be only sparsely occupied, we find primarily dense, broad filling of the morphological landscape, indicating rich diversity. Phylomorphospace plots suggest that dynamic evolution underlies the observed spatial patterning. Evolutionary transits through the morphospaces are sometimes long, and occur in a variety of directions. The trajectories seem to represent both evolutionary divergences and convergences, the latter supported by convevol analysis. We suggest that that this pattern of occupancy reflects the various adaptations of different groups of fishes, seemingly paralleling their diverse marine and freshwater ecologies and life histories. Opercle shape evolution within the acanthomorphs, spiny ray-finned fishes, appears to have been especially dynamic.

About the Authors:

Charles B. Kimmel

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Institute of Neuroscience, University of Oregon, Eugene, Oregon, United States of America

ORCID http://orcid.org/0000-0002-1122-4371

Clayton M. Small

Roles Formal analysis, Writing – review & editing

Affiliation: Institute of Ecology and Evolution, University of Oregon, Eugene, Oregon, United States of America

Matthew L. Knope

Roles Formal analysis, Writing – review & editing

Affiliation: Department of Biology, University of Hawaii at Hilo, Hilo, Hawaii, United States of America

Introduction

Understanding the patterns of morphological diversifications among organisms constitutes a central goal of evolutionary biology. These patterns can be quantified using multivariate approaches such as geometric morphometrics [1,2] and visualized as plots on spaces of reduced dimensionality. Such spaces are known as morphospaces, and the occupancy within the space may be termed the morphospace landscape (i.e., denoted by the collective of datapoints of the scatterplot). Of special interest in our study is that it has been understood from previous theoretical and empirical studies that a morphospace landscape generally is not uniform in appearance, but sparsely occupied and clumpy [3–10]. The clustering is claimed to be independent of scale–“from populations to the highest taxonomic categories” [8], which if true would seem remarkable. What kind of evolutionary process(es) could result in such a broadly shared evolutionary pattern?

Within a morphospace framework, evolution can be thought of as a suite of processes—adaptive or otherwise—that modify landscape features. The change in position between an ancestor and a descendant would be visualized as a trajectory, a specific pathway, through the morphospace. This phylomorphospace approach, that we use here to visualize these trajectories [11], requires not only being able to measure morphologies (e.g. of extant forms or well-preserved fossils) such that the morphologies can be plotted, but also knowing the phylogenetic relationships of the forms in question. With a phylogeny in hand the positions of ancestors, essentially unavailable for direct measurement, even if there is something of a fossil record, can be estimated by computation (e.g. by squared change parsimony; [12]). These computed ancestral positions in the landscape are then compared with the measured positions of the descendent set. The phylomorphospace, showing the trajectories, provides a way to understand trait evolution.

In this study we examine the shapes, and their evolution, of a single skull bone, the opercle (OP), among an extremely broad group of organisms, the teleost fishes. Teleosts, by far the most diverse group of vertebrates, are classified into 66 orders (sensu [13]), more than 500 families, and more than 30,000 species. The morphological differences among the head skeletons of these species can be quite dramatic [14]. Our interest in examining just a single bone is because it is possible, as we have shown for zebrafish, to consider how it is shaped in considerable detail, at the level of cell biology and developmental genetics [15–18]. This means that the evolutionary patterns we observe, here at the macroevolutionary scale, may eventually be understandable at a very refined and mechanistic level through complementary approaches from comparative and developmental biology. The OP, providing the principal support of the gill cover, functions critically in respiration and mouth opening [19,20]. Hence, we can contemplate possible ecological determinants of morphological change. Also motivating our current study was our exploration of evolutionary shape change within threespine stickleback. We found the local morphospace representing the oceanic-freshwater divergence of a number separate populations across the world, to be rather densely occupied (S1 Fig; [21]). Does this microevolutionary pattern have a counterpart in macroevolution?

Here we extend the stickleback work by examining opercle disparities among teleost families, in contrast to the within-species stickleback data. Earlier studies have broadly examined single kinds of bones broadly among fish species, but with different motives than ours. For example, McAllister examined branchiostegal rays in an attempt to use them as systematic characters (i.e., to determine phylogenetic relationships) [22], and Arratia and Schultze focused on homology of a hyoid tendon bone, the urohyal, describing it as a unique feature of teleosts [23]. We choose the family taxonomic level for this analysis because, as we show below, the shape differences among the families are prominent, and as expected, considerably more so than within species, or even within families as revealed by recent work examining Cichlidae and Ariidae [24,25]. To look for directional channeling of evolutionary shape change, especially evolutionary convergences, we use a phylomorphospace approach. We find mostly dense and fairly unclustered morphospace distributions. The phylomorphospace findings suggest the presence of marked evolutionary excursions through space, including both divergences and convergences. Especially considering recent work describing “explosive” responses of acanthomorphs to fill in a region of morphospace left barren by the end Cretaceous mass-extinction event [26,27], our findings lend support to ecology as being a major determinant of morphospace patterning, again matching interpretations of the possible cause of stickleback microevolution.

Materials and methods

Ethics

The study was approved by the University of Oregon IACUC. No living animals were utilized or euthanized specifically for the purpose of this study: The preserved fish (or just their heads) were obtained as by-products from other projects, including gifts from individuals or gifts or loans from scientific collections (US Marine Fishery, NOAA; Oregon State University Ichthyology Collection) and otherwise (e.g. sport fishery, hobby tropical fish suppliers).

Specimen collecting and preparation

Our survey is based on species from 110 families, therefore representing about a fifth of the approximately 500 families of Actinopterygians, ray-finned bony fishes. All but three families are teleosts. The species and the families they represent are broadly distributed across the phylogenetic tree, however, virtually no deep-sea groups are represented. 103 of our families are included in the tree of Betancur-R et al. [13], which we used to construct the phylogeny for our analyses. Dissections, photographic imaging and morphometric analyses were initially carried out with 237 preparations representing 202 species: eventually this dataset was reduced to single samples for each species by averaging shape coordinates in cases where we had multiple individuals from the same species. In the interest of equal taxon sampling, for families for which we had multiple representatives, we culled to a single ‘exemplar’ species for each family; the final list is reported in S1 Table. Because quantitative bone shape studies have not previously been done for the great majority of species studied here, we have no way of knowing whether the chosen examples are truly representative of the groups. Examples of the within-family variation as compared with among family variation for our study set are shown in S2 Fig.

We dissected the OP from the skull, noting positional relationships with neighboring bones and the locations of muscle attachments to aid in eventual identification of landmark positions. Dissection, allowing the individual bones to be examined in detail in their entirety, is by far our preferred method for analyzing bone morphology and landmarking. Dissection was sometimes aided by digestion of soft tissues with trypsin and cleaning by dermestid beetles. We also preferred to avoid formaldehyde (formalin) preserved material, which toughens soft tissue, including ligaments, and makes dissections of the sometimes very delicately mineralized bones more difficult. For this reason, and because dissection is destructive we mostly avoided using museum specimens. We obtained most specimens as fresh frozen preparations. Some specimens were preserved in 70% or 95% ethanol. Smaller species were briefly fixed or postfixed (generally for 2–4 hours) with 2 or 4% neutral formaldehyde, and then double-stained with Alcian blue–Alizarin Red, following a procedure that minimizes mineral loss (and hence preserves Alizarin Red staining) during preparation [28]. Following staining and clearing, the bones were separated by dissection in a solution of 50% glycerol and 0.01% KOH (final concentrations). Whether the bones were stained or not, we did not allow drying throughout the procedures: Wet preparations are preferable for revealing subtleties of bone structure which aids in landmarking the bones (Fig 1).

[Figure omitted. See PDF.]

Fig 1. Landmarking the OP.

The distribution along the bone edges of four landmarks (LM, yellow, red border) and 14 semilandmarks (yellow) used in this study, in a smelt (Osmerus mordax, Osmeridae). Developmental studies in zebrafish, stickleback, cory catfish, and Oncorhynchus species inform the positions of the landmarks, those that appear to correspond among all of the species examined [15,29,30], and unpublished: These landmarks include LM1 (joint-associated), LM2 (socket of the ball-and-socket joint the OP makes with the hyomandibula), LM6 (ventral apex) and LM12 (posterior apex, associated with a prominent strut that Gregory termed a tract of folded trabeculae [14] traversing the bone from the joint region). The opercular dilator muscle connects to the OP at LM1. The opercular levator muscles connect between LM1 and LM12. The prominent curvature at LM9 is also shared among many species, including zebrafish (termed ‘c’ in [15]. The semilandmarks were placed either at constant intervals along the appropriate segment of the edge, or to capture elaborations of the edge in particular species. There is no basis currently for supposing that any of the semlandmarks capture homologous locations among the families, even though some them appear consistent among species within families.

https://doi.org/10.1371/journal.pone.0188888.g001

The OPs, usually from both the left and right sides of individual skulls, and often along with other selected craniofacial bones, were photographed, including both medial and lateral views. Digital images were taken with a Leica stereomicroscope equipped for epifluorescence for the smaller preparations and with a Nikon macro setup for the larger ones. In either case we used an incident light source, offset for the bright field images to provide relief, which aids in morphological study. Comparisons of the left-right morphologies of the bones themselves, as well as the images, provided a check for abnormalities that might arise as developmental accidents or due to injury. In fact, we essentially never encountered major left-right differences. A single OP from each specimen was selected for morphometric analysis.

Morphometrics and data analyses

The OP, supporting the gill cover (operculum) is quite flat, such that 2D rather than 3D morphometric analyses are appropriate for shape analyses. Eighteen landmarks were placed around the edge of the OP for geometric morphometric analyses (Fig 1, [31,32]). The legend to Fig 1 provides detail concerning landmark placement. Landmark positions were marked on the digital images in Photoshop, and then digitized using tps Dig, version 2.12 [33]. The digitized raw datasets were moved to MorphoJ [32] for Procrustes transformation, which removes size and orientation from the bone shape analysis, and the resulting shape variables were ordinated and examined by principal component analysis (PCA). Because the first two principal components explained more than 50% of the total shape variation, along with a precipitous drop in variation explained by the remaining PCs, we based subsequent analyses on the PC1-PC2 morphospace. PC scores were transferred to JMP (SAS Institute Inc., Cary, NC, v. 8) for examining scatter plots including nonparametric density landscapes. The PC scores and time-calibrated phylogeny pruned from a massive phylogenetic tree of the teleosts [13] were then transferred to Mesquite [34] for visualization of phylomorphospace [11]. We also used morphol.disparity and compare.evol.rates functions from the R package geomorph (version 3.0.0) [35] for comparing morphological disparities [31] and evolutionary rates of divergence [36,37] in a phylogenetic context. We tested for phylogenetic signal using “K” metrics [38], determining a multivariate version of K from the Procrustes coordinates [39], as implemented in geomorph. Using the R package GEIGER [40,41], we evaluated the fit to our PC data (the first two principal components) of four models of phenotypic evolution: “early burst” [42,6,43], Brownian motion, white noise, and Ornstein-Uhlenbeck. We also tested for phenotypic convergence in the PC1-PC2 phylomorphospace using the R package convevol [44] (see [45] for a working definition of convergence, as used in this study). Specifically, we tested for convergence in groups of taxa associated with density features of the landscape, including high- (“peaks”) and low- (“valleys”) density regions. Hypothesis tests for convergence, using the metrics described in Stayton [44], were carried out in convevol via 1000 Brownian motion simulations per test.

Results

Diverse opercle morphologies

In agreement with classical work [14], we find prominent shape variation of the OP, the principal bone of the gill cover (operculum) of ray-finned fishes (Actinopterygii) (Fig 2). The bone is usually broad and sheet-like in planar view, but there are evident exceptions, e.g., Fig 2D. Some OPs exhibit prominent protrusions serving as attachment points to musculature or other skull elements, or apparently serving as defensive spines (e.g., Fig 2B). In all of our samples, the OP is flattened within the plane of the operculum. The socket of the ball-and-socket joint the OP makes with the hyomandibula (shown to the upper left for each example in Fig 2) is invariably present and appears conservative in structure. Other than the flattening and joint socket, we see marked diversity of opercle form. We note that shape deformations do not appear to be at all uniform across the bone edges. Spines, if present, generally adorn the posterior and ventral bone edges, sensible locations for protection since they are present over the region of opening of the gill chamber. Interestingly, a straight stretch of the anterior-ventral edge of the OP (between landmarks 2 and 6; Fig 1) seems highly conservative in shape among our samples (see Discussion), with a few exceptions (e.g. Fig 2E).

[Figure omitted. See PDF.]

Fig 2. Opercles (OPs) show exceptional morphological diversity across species and families.

A Salmon, Oncorhynchus gorbuscha, family Salmonidae, B sand bass, Paralabrax nebulifer, Serranidae, C reineta, Brama australis, Bramidae, D batfish, Dibranchus atlanticus, Ogcocephalidae, E toadfish, Opsanus tau, Batrachoididae. All of these examples show teleosts, and all but A are spiny-ray finned neoteleosts, acanthomorphs. Scale bars: 5 mm. F Principal component ordination arranges the OPs in a morphospace according to their shapes. The presentation emphasizes the exceptional variation we encountered in this survey. Here the bones, presented as silhouettes, are stripped of all of their distinguishing features (the location of the joint socket, presence of struts crossing the bone, outgrowth incremental banding pattern, density of Alizarin Red staining) except their shapes. We selected the examples to show the spread and mapping across the PC1 by PC2 morphospace. The white outline surrounding the occupied region of the space, the ‘landscape’, was made from nonparametric density contours. See text and Fig 3).

https://doi.org/10.1371/journal.pone.0188888.g002

Features of morphospace occupancy

A PCA-derived morphospace plot reveals a landscape that mostly looks densely occupied. Although some clustering is evident among the samples, and thus their distribution is not at all even, the landscape is not notably clumpy (Fig 3). To provide a check that this filled in pattern is not due to some technical artifact, we used an independent method of shape analysis of our samples, based on outlines rather than landmarks, and obtained a similar dense occupancy pattern to that in Fig 3, including a largely matching up of the spreads of the individual samples across the landscapes (S2A and S2B Fig). Similarly, lack of prominent clumping does not depend on our examining only a single species for each family. We also examined a morphospace in which some of the families are represented by multiple species (S2C Fig). This multiple species morphospace includes about twice as many data points and accordingly is more filled in (compare S2A and S2C Fig), but its pattern of distribution of points is otherwise reasonably similar to the single species one. The unequal taxon sampling present for the multiple species analysis might be expected to increase clustering, if anything.

[Figure omitted. See PDF.]

Fig 3. Principal component-based morphospace occupancy.

PCA reveals that most of the bone shapes densely occupy a single large region of the space, here constructed from the two leading eigenvectors, PC1 and PC2, derived from landmark analysis. The percentages within the parentheses give variances explained. The accompanying overlay configurations show shape changes associated with the PC axes. PC1 captures details of the outlines: forked OPs toward high positive PC1 values and smooth configurations, often with prominent dorsal regions (upper in the figure) toward negative PC1 values. In contrast PC2 separates bones that are thin along the anterior-posterior axis (positive PC2 or much broader along this axis (negative PC2. The data points of the scatter plot show single exemplar species from each of the families across the phylogeny in Fig 4, with matching colors between the points here and the branches of the tree. The lines included with the scatterplot show nonparametric density contours. Heat map colored contours are according to occupancy density, with warm colors indicating density peaks. Morphospaces using other PCs (e.g. PC3 by PC4) also reveal dense occupancies.

https://doi.org/10.1371/journal.pone.0188888.g003

As pointed out above, the distribution of data points in Fig 3 is not completely uniform, perhaps due in part to taxon sampling, but it also could have biological meaning–i.e., reflecting a regionally rugged morphological landscape, with multiple peaks and valleys. The presentation in Fig 3 conveys these rugged regions by including contour lines that show the nonparametric densities of local space occupancy. We used color coding of the data points to see, at large phylogenetic scale, how families broadly related to one another are distributed within these occupancy landscapes. The colors of the points in Fig 3 correspond to clades across the phylogenetic tree in Fig 4. For example, red indicates the neoteleosts, the largest teleost clade, arising deep within the teleosts. In our dataset, all of the neoteleost families except the basal-most one, Synodontidae, are spiny ray finned fishes, acanthomorphs, and they are mostly located on the right side of the landscape–toward high PC1. In contrast, most nonacanthomorph groups are present to the left, especially the upper left (low PC1, high PC2). Except for the ‘far right’ which is exclusively red, mixing among the groups represented by the various colors is considerable. The three families color-coded green, and representing the most basal teleost clade, are located toward the bottom of the landscape (low PC2), and share this location with acanthomorphs. The acanthomorph lineages outnumber the others, which partly accounts for their broader distribution. However, particularly the exclusive occupancy of the right side of the landscape (Fig 3, high PC1 values) by acanthomorphs strongly suggests that among the sampled taxa, their OP morphologies are rather more diverse than in the nonacanthomorphs examined. In agreement, geomorph tests using the morphol.disparity and compar.evol.rates functions, respectively, suggest that acanthomorph OP morphological disparity, estimated as Procrustes variance [31] is significantly higher than that of the nonacanthomorphs included in this study (1.52 times, p = 0.002; Table 1) and, further, that acanthomorph OP shapes have evolved significantly more rapidly than OPs of the nonacanthomorphs (2.33 times, p = 0.001; Table 1, [36,37]. We measured the tendency for related taxa to resemble one another, known as phylogenetic signal [38], using a multivariate version of Blomberg’s K [39]. We found K to be significant and relatively low, 0.53 (p = 0.0001, determined from permutation analysis). We can interpret low phylogenetic signal to be consistent with our visual impression of scattering among lineages, for scattering would reduce phylogenetic signal.

[Figure omitted. See PDF.]

Fig 4. The time-tree showing families analyzed in this study and their times of origin.

The tree was pruned from an extensive species-level time tree provided to us by Richard Betancur-R, using data from [13]. Branches along the tree indicate families. Color coding shows clades across the tree, and matches Fig 3.

https://doi.org/10.1371/journal.pone.0188888.g004

[Figure omitted. See PDF.]

Table 1. Morphological disparities and morphological divergence rates of the OP in acanthomorphs and nonacanthomorphs.

https://doi.org/10.1371/journal.pone.0188888.t001

Modeling with GEIGER [40, 41] yields a first-pass understanding of evolutionary dynamics that potentially could help to explain features of the morphospace patterning. As shown in Table 2, white noise, essentially a distribution not dependent on phylogeny, was not supported. The analysis also yielded only little support for “early burst” [43], a model characterizing ‘Simpsonian’ adaptive radiation in which rapid evolutionary changes occur approximately synchronously during the first stages of a radiation [42]. There is also little support for the “random walk” model of Brownian motion. Arguably, the Ornstein-Uhlenbeck model, positing adaptive directional evolution, is the best supported, considering PC1 and PC2 collectively.

[Figure omitted. See PDF.]

Table 2. Tests of evolutionary models using the GEIGER platform, based on individual Principal Components.

https://doi.org/10.1371/journal.pone.0188888.t002

Evolutionary convergences occur abundantly

We saw above that samples from different lineages can occupy closely neighboring positions on the morphospace landscape. For a striking example, consider the principal occupancy peak toward low PC1 and low PC2 on the landscape (i.e., to the lower left in Fig 3). Ten families are present above the highest contour line (Polypteridae, Lepisosteidae, Notopteridae, Ariidae, Osmeridae, Gobiidae, Atherinopsidae, Channichthyidae, Nototheniidae and Triglidae). The ten families include eight orders (sensu [13]) and nearly span the tree shown in Fig 4. This cohort, along with other examples, motivate us to entertain evolutionary convergence as a feature of the patterning of the landscape. We use convergence in the broad sense, to mean a pattern as defined by Stayton [41], and not to imply any underlying process such as adaptation.

We used the platform convevol [44] to directly examine whether convergence is a significant determinant of morphospace occupancy (Fig 5). The approach relies on the analysis of phylomorphospaces [11], which use morphological disparities, along with the phylogenetic connectivities and times of divergences of the branches of the phylogenetic tree to portray computed evolutionary pathways through morphospace. The positions of ancestors (branch points, indicated by smaller points in the Fig 5 plot) are calculated by squared-change parsimony, and the positions of the terminal taxa (larger points at the terminal ends of the branches in Fig 5) are empirical, matching the points of the Fig 3 landscape plot. The convevol output metric that we consider most important for our dataset, termed C5, describes the frequency of convergences into user-specified regions of a phylomorphospace, based on the computed positions of ancestors. The example in Fig 5 illustrates our findings. The phylomorphospace shows convergences onto the major landscape occupancy peak described above, here indicated by the lavender ellipse encompassing 13 families–the lower panel shows a zoomed-in view. Convevol detected 12 independent convergences into the region (red arrows), significant at an alpha level of 0.05, as determined from 1,000 simulations based on Brownian motion (P = 0.049).

[Figure omitted. See PDF.]

Fig 5. Convevol frequency-based measurement yields evidence for convergent evolution on the principal landscape occupancy peak (termed peak 1 in Table 3) of PC1 by PC2 phylomorphospace.

The large points show terminal taxa (families) corresponding to points of the landscape plot shown in Fig 3. The smaller points show inferred positions of ancestors. Note the evolutionary trajectories, many of them quite long, that crisscross the landscape in a variety of directions. The ellipse (lower panel shows a zoom in) indicates the occupancy peak of interest, and the red arrows show independent entries (convergences) into this region. See text and Table 3 for further explanation.

https://doi.org/10.1371/journal.pone.0188888.g005

These data are compared with seven similar analyses of the landscape, and examining not only the frequency-based measure C5, but also four similarity-based measures C1-C4 provided by convevol (Table 3; Fig 6 shows the analyzed locations on the landscape). C1, for example, quantifies the ratio of dissimilarity between terminal taxa to the maximal dissimilarity between any taxa in the two lineages; convergence is indicated where the terminal taxa are more similar. Of the eight regions analyzed, we selected four to represent the four prominent occupancy density peaks present on the landscape, and other four to represent off-peak regions (i.e., regions of lower elevations, such as the regions bordering the lowest elevation valleys; the very lowest troughs did not include enough datapoints for analysis). Peak 1, the example used in Fig 5, achieved significance for all of the metrics (C1 -C5), as did peak 2 and peak 4. Convergence onto peak 3 was only weakly supported, showing significance for only C1 and C2. Perhaps surprisingly, three of the four off-peak regions showed significance for all of the metrics, and the other, region 6, achieved significance for all of the metrics except C5. With six of the eight regions analyzed showing statistical support for all the metrics, and the other two showing partial support, convevol thus provides strong evidence for convergences occurring at multiple locations across the landscape.

[Figure omitted. See PDF.]

Fig 6. Key to the peak (in red) and off-peak (black) regions of the morphospace landscape analyzed by convevol, and reported in Table 3.

Fig 3 shows the details of the morphospace.

https://doi.org/10.1371/journal.pone.0188888.g006

[Figure omitted. See PDF.]

Table 3. Frequency (C5) and similarity (C1-C4) based measures and probabilities (p) examining convergences within the indicated focal regions of the morphospace.

Convevol analyses.

https://doi.org/10.1371/journal.pone.0188888.t003

Both evolutionary convergence and divergence appear to contribute to morphospace patterning

In the preceding analyses, we examined particular small regions of the morphospace landscape for evidence of independent convergences among the neighbors occupying these regions. We can turn this approach around, and look at the morphological trajectories of OP shapes along lineages that are neighboring on the phylogenetic tree with respect to the entire morphospace. That is, what are the OP morphological differences among closely related families, and what does the phylomorphospace approach suggest about the nature of the divergences responsible for these differences?

To accomplish this analysis, we pruned small sets of neighboring branches out of the phylogenetic tree, with new color coding to aid in making out the details of the patterns of trajectories within these small sets. Exemplifying this approach, Fig 7A shows the phylomorphospace representing 18 families from the base of the neoteleosts. Fig 7B shows their phylogenetic relationships. Close relatives show loose clustering in the phylomorphospace (e.g., the blue clade encompassing four gadiform families and the sister group Zeidae). In contrast, the trajectories are not seemingly channeled, but scattered with respect to both direction and length. Syngnathidae (Syng) shows a particularly prominent long terminal trajectory on the phylomorphospace. The accompanying landscape plot (Fig 7C) reveals that the eighteen families of this clade are scattered across nearly the entire morphospace in the vertical dimension (representing PC2). Besides the divergences, apparent convergences are also suggested (e.g. Synodontidae (Syno) and Bramidae (Bram) converge toward one another at the top of the landscape).

[Figure omitted. See PDF.]

Fig 7. Phylomorphospace analyses reveals extensive evolutionary divergences, as well as apparent convergences among OP shapes in basal acanthomorphs.

A The families included in the plot represent a contiguous paraphyletic grade of Eurypterygian fishes, as shown in the pruned tree in B. Fishes in the sister group to Scombromorpharia (sensu [13]) have been pruned from the tree. Colors facilitate tracing separate lineages, as explained in the text, across the landscape. C The landscape plot (points representing individual families as in Fig 3) shows the distribution of these families (color matched). For explanation of the abbreviations, see text and S1 Table.

https://doi.org/10.1371/journal.pone.0188888.g007

This pattern, suggesting both divergences and convergences, occur irrespective of what set of branches we analyze across the entire phylogenetic tree. We show another example that examines a perciform clade sampling the ‘bush at the top’ (as the collective of taxa in this terminal region of the phylogeny is frequently called) (Fig 8). As in Fig 7, we see as convergences–Nototheniidae (Nott), Channichthyidae (Chan) and Triglidae (Trig) onto occupancy peak 1, and Serranidae (Serr) and Scorpaenidae (Scor) onto peak 2. And we see, perhaps more prominently than convergences, long trajectories separating sister families, that are suggestive of marked divergences (Cirrhitidae (Cirr)–Centrarchidae (Cent) and Liparidae (Lipa)–Cottidae (Cott)).

[Figure omitted. See PDF.]

Fig 8. Phylomorphospace analysis of the ‘terminal’ (distal-most) clade of acanthomorphs, including Cottoidea and their allies.

Presentation as in Fig 7. For explanation of the abbreviations, see text and S1 Table.

https://doi.org/10.1371/journal.pone.0188888.g008

We suggest from our phylomorphospace analyses that both divergence and convergence need to be considered as prominent mechanisms underlying landscape patterning and its ruggedness. We make this interpretation with caution, however, since using the phylomorphospace approach is not without errors, more substantial in the earlier stages of the lineages [11, 46].

Discussion

Teleost opercle macroevolution has produced a richly occupied morphospace

A principal finding of our study is that the local morphospace for the OP within teleosts is largely filled. Might dense occupancy be due to limited diversity, i.e., that OP shape is not evolving to any large extent? This potential explanation is not supported by our study: Our study set, while only representing about a fifth of teleost families, essentially spans the full phylogenetic tree of this numerous and richly variable group of vertebrates, although some groups such as deep-sea lineages are poorly represented. Visually comparing the bone shapes suggests the variation to be high. Furthermore, measurement of disparity using Procrustes variance suggests high disparity: For a rough comparison, we calculated OP morphological disparities using the Procrustes shapes of a single pair of oceanic and freshwater stickleback populations of three-spine stickleback from the study shown in S1 Fig: Our estimate for the microevolutionary stickleback disparities were 10-fold less than the values for our macroevolutionary dataset. To illustrate the difference directly on the morphospace examined in this study, S2 Fig (panel C) shows two points representing the disparity of the two stickleback morphs, as compared with the spread of the full species dataset.

Our finding of a well-filled morphospace in which clustering is not prominent is surprising, given the general understanding, since the pioneering work of Raup [3,4] and see [47], that morphospace occupancy is sparse [48,49,9]. In this model, occupied regions of the space, the ‘clumps’, are contained within substantial regions of unoccupied space. The lack of occupancy of the latter regions is presumed to be due to some kind of constraint on evolution. For example, lack of occupancy might be due to functional constraint–i.e., the morphologies that might have occupied the vacant regions are maladaptive [9,47]. Alternatively, developmental constraint could make morphospace regions inaccessible by the developmental process; [50,51,5]. The PCA-based empirical method that we examine in this study should have revealed a mostly empty morphospace if prominent constraints indeed were present: Examples of sparsely occupied morphospaces from our own work are the morphospaces describing OP development in zebrafish and stickleback. In both, variation is limited to a restricted zone of the space around nonlinear developmental trajectories, similar in form in both species. The trajectories completely avoid most of the space, including the central regions (Fig 5 in [15]; Fig 5 in [30]). Nevertheless, the unoccupied regions of the space we observe in the current study are not impressively large. In fact, they may result from our incomplete sampling of teleost taxa [9], and are not due to constraints.

An important dimension in the understanding of the evolution of our observed morphospace patterning is the addition of the phylogeny to its visualization–the phylomorphospace. A caveat is that phylomorphospaces are based on the positions of ancestors computed by squared change parsimony, not measurement of the morphologies of known ancestors. For this and other reasons the phylomorphospaces are likely to be only inaccurate portrayals of actual evolutionary changes (see also [11,46]). In spite of such problems, several examples in our study argue for both evolutionary convergences and divergences abundantly occurring within the morphospace. Both the convergences and divergences occur in a variety of directions, and frequently by long trajectories across the space. Both convergence and divergence have seemingly participated in a major way in filling the landscape.

Developmental modularity could provide flexibility in shape evolution

Given the substantial variety of freshwater and marine habitats that different teleosts occupy, it seems likely that varying ecology and Darwinian adaptation have prominently influenced the patterning of the OP morphospace. Development might also play a role. In a study addressing opercle development in zebrafish we discovered that the different edges of the opercle grow at separate times from one another. Edge-specific growth is mediated by selective acquisition of early sp7-expressing osteoblasts [15]. A mutational study showed that the osteoblasts along one of these edges specifically, use a Hedgehog signal to stimulate cell division in an adjacent pre-osteoblast population. The signal thus provides regulation of the size of the pre-osteoblast pool, and indirectly, the rate of local growth of the bone [16]. These observations suggest that patterning morphology of the bone is modular, where developmental modules are defined as semiautonomous units of morphological patterning [52–54]. In this case, we proposed that modularity is controlling time- and position-specific features of local bone formation and growth–and hence, region by region, the shape of the opercle [15,16,55].

We extended the zebrafish studies to sticklebacks, comparing development in ancestral (oceanic) and derived (freshwater) populations, and showed that opercle evolution in this species involves a kind of “dissociated heterochrony” [56] in which at one developmental stage, a single region of the bone develops at different rates in the two forms. Strikingly, examining the opercle’s covariance structure by Klingenberg’s method [57] provided direct evidence for modularity, and found that the evolutionary change in bone shape is along a module boundary, the location of the boundary reasonably well predicted from the zebrafish work [30].

This previous work, with single species, informs macroevolution [55]. Modules, as semiautonomous units, would allow for regional dissociability, providing evolutionary ‘flexibility’ among species and higher taxa, and potentially channeling evolution in different directions in different taxa. We propose that the rich morphospace occupancy and the diverse trajectories through the space that we observe depend in part on modular development. As we better understand the developmental genetic underpinnings of these modules, one hope is that future comparative work leveraging detailed genomic information from a host of teleost species will help explain how module tinkering contributes to biodiversity over evolutionary time.

An evaluation of possible constraints on OP shape

That the OP is invariably a flat bone, and invariably forms the socket for the ball-and-socket joint with the hyomandibula shows constraint on morphology. One can easily imagine that these are functional constraints, given that the bone supports a flattened operculum, the movement of which mediates gill-ventilation. The OP also makes close associations with other bones, notably the subopercle, but in spite of this potentially constraining influence, we observed considerable shape variation of this OP edge, the posterior-ventral edge where it articulates with the subopercle. In contrast, we noted in the Results section the apparently low shape variation of the anterior-ventral edge of the OP. This edge abuts the preopercle, connecting with it via a ligament, and movement across the edge occurs during respiration. Accordingly, one might suppose that shape conservation represents functional constraint. However, from our zebrafish study, this OP edge is seemingly an early module in itself [14]; and one might alternatively suppose that variability of the module and edge shape is constrained developmentally, rather than functionally. There exist several routes toward attempting to find out whether a putative constraint is developmental or otherwise [9, 58]. In what he termed “the logic of monsters”, Alberch [58] points out that one might observe a maladaptive and strange phenotype in artificial situations, and proving that the phenotype is possible developmentally even though not observed in nature (see also [9]). In analyses of zebrafish homozygous for a lethal mutant allele of the transcription factor-encoding gene mef2ca we find that the anterior ventral edge of the OP can re-shape dramatically (Fig 9). By the logic of monsters argument, this finding suggests that developmental constraint is not the best explanation for why this edge is conservative among naturally occurring variants.

[Figure omitted. See PDF.]

Fig 9. In this zebrafish mef2ca mutant a bridge of bone (*) forms between the OP anterior-ventral bone and a branchiostegal ray (BSR).

Early larva stage–mutants die soon after this time (examples from [17]; see also [18]).

https://doi.org/10.1371/journal.pone.0188888.g009

Dynamic evolution of the opercle in acanthomorphs

Adaptation might be the principal driver of opercle morphology. We observed prominent occupancy peaks across the landscapes which could be due to adaption, but we note that equating adaptive and occupancy peaks requires additional supporting evidence [47,59]. Our finding multiple convergences with the convevol platform might also suggest adaptation, but this interpretation also is not compelling, since convergences could be due to constraints, or even to neutral processes [45]. Our GEIGER analysis may be relevant to this issue in that an adaptive model, Ornstein-Uhlenbeck, is supported using each of the leading two principal components of OP shape in our trait analysis.

Evolution of opercle shapes in acanthomorphs may have been especially dynamic. The acanthomorph OPs spread more than other groups across the landscape, and the phylomorphospace plots suggest prominent trajectories occur throughout. Previous studies of body morphology [26,27] see also [60] demonstrate that acanthomorphs have opportunistically and ‘explosively’ occupied morphospace newly vacated by the end-Cretaceous extinction event. Our results seem to complement these recent paleontological and comparative findings. Further, we wonder whether, in addition to the body forms studied by Friedman and associates, our results could mean that explosive morphological evolution within the acanthomorphs extends to head bone shapes. Further study of skull bone morphologies, particularly as related to the ecology of lineages of acanthomorph fish could address this question.

Supporting information

[Figure omitted. See PDF.]

S1 Table. Families and exemplar species studied.

Branch number refers the terminal branches of the pruned phylogenetic tree [13]. 0 indicates the family was not included in the Betancur-R study.

https://doi.org/10.1371/journal.pone.0188888.s001

(DOCX)

S1 Fig. Morphospace showing microevolution in stickleback.

Global microevolutionary OP divergences in threespine stickleback obtained from principal component analysis (PCA) yield a well-filled local morphospace. Red filled circles: oceanic (ancestral form, 8 populations). Black open circles, freshwater (derived forms, 14 populations. Notice that within the single cluster of OP forms, the two morphs are mostly separated along PC1. The configurations (right) show the deformation along PC1 (gray: oceanic, black: freshwater, arrows indicating prominent shape change). The plot is original, the data (available from Dryad) are from [21].

https://doi.org/10.1371/journal.pone.0188888.s002

(TIF)

S2 Fig.

Comparability of OP datasets A-C), and examples of within-family variation (D). Morphospaces as in Fig 3 of the main text. A The same landmark-based data with sample points repeating the presentation Fig 3 of the main text. B shows the same samples, but with a PCA morphospace using an outline method based on Elliptical Fourier transformation. The outlines were captured by digitizing 400 xy points for each sample, moved to PAST [Hammer Ø, Harper DAT, Ryan PD. PAST: Paleontological Statistics Software Package for Education and Data Analysis. Palaeontol. Electronica 2001;4: 9pp], and transformed into 2D shape coordinates including 30 modes (harmonics). Size information is removed by this procedure. Silhouettes in A and B show the sample shapes and colors show phylogenetic relationships, matching Fig 3 and the branches of the tree in Fig 4. The individual OPs generally map to approximately the same regions of both spaces, such that samples near one another on one of the spaces, are also usually also near one another on the other. Whereas the elliptical Fourier and landmarking methods are clearly revealing different aspects of Op shape, there is also some consistency between the two kinds of space-shape mappings. C Including more species in PCA morphospace using landmarking does not yield a substantially more clustered distribution of samples. Rather the space is more densely occupied, as compared with the plot made from only single species per family (A). The two larger points show overlays representing oceanic (filled, black circle) and lake (open circle) stickleback, described in the Discussion. The configuration diagrams indicate shape changes associated with the PC axes (consensus shape in light blue in the overlays). D The same morphospace as C, but with the landmarked species of five families highlighted to show examples of within-family shape variation, as compared with among-family variation represented by the whole plot. Colors match the other panels. Cichlidae (Cich) shows the most extreme within-family variation encountered in this study. The others are more typical: Characidae (Char), Salmonidae (Salm), Labridae (Labr), Scombridae (Scob).

https://doi.org/10.1371/journal.pone.0188888.s003

(TIF)

Acknowledgments

We are grateful to Brian Sidlauskas for counseling and help with the analyses throughout nearly the entire course of this study. We also acknowledge Brian’s willingness to accept our set of dissected opercular bones, including the opercles used in this study, into the Ichthyology Collection at Oregon State University, where they will be archived and made freely available for study. Reida Kimmel contributed key insights about development. We thank Richard Betancur-R for supplying the time tree from [13], Eric Hanneman, Douglas Markle, Kevin Meyer and Kevin Clifford for expert species identifications, and Charline Walker for help with dissections. Patrick Phillips, Bill Cresko, Brian Sidlauskas, Philip Hastings, and anonymous reviewers read and provided comments on earlier drafts of the paper. Bill Cresko contributed to hours of fruitful and wonderful discussion. We especially acknowledge a large number of individuals for supplying without cost to us the dead fish, either entire bodies or just heads, from which we made all of the preparations–Angel Amores, Margie Armstrong, Peter Armstrong, Susie Bassham, Ingo Braasch, Sam Chapman, Bill Cresko, Mark Currey, Thomas Desvignes, John Dowd, Eric Hanneman, Laura Hisrich, Jamie Kidder, Dick Kimmel, Reida Kimmel, Seth Kimmel, David Lains, Juan Manuel, Braedan McCluskey, Bill Mebane, Kevin Meyer, Mike Miller, Jamie Nichols, John Postlethwait, Bonnie Ullmann, Connie Wagner, Dave Wagner, Charline Walker, Sawyer Watson, Monte Westerfield, Harold Zakon, and Steve Zottoli. The authors have no conflict of interest to declare. C.B.K.’s interest in this project was sparked during a sabbatical leave at the University of Cambridge.

Citation: Kimmel CB, Small CM, Knope ML (2017) A rich diversity of opercle bone shape among teleost fishes. PLoS ONE 12(12): e0188888. https://doi.org/10.1371/journal.pone.0188888

References

1. Mitteroecker P, Gunz P. Advances in geometric morphometrics. Evo Biol 2009;36: 235–247.

2. Klingenberg CP. 2010. Evolution and development of shape: integrating quantitative approaches. Nat Rev Genet. 2010;11: 623–635. pmid:20697423

3. Raup DM, Michelson A. Theoretical morphology of the coiled shell. Science 1965;147: 1294–1295. pmid:17790826

4. Raup DM. Geometric analysis of shell coiling: general problems. J Paleont. 1966;40: 1178–1190.

5. Alberch P. Developmental constraints in evolutionary processes. In: Evolution and Development. Dahalem Konferenzen. (ed J. T. Bonner JT). pp. 311–332. Berlin, Springer-Verlag. 1982.

6. Gould SJ. Wonderful life: The Burgess Shale and the nature of history. New York: WW Norton, 1989

7. Foote M. The evolution of morphological disparity. Ann. Rev. Ecol. Syst. 1997;28: 129–152.

8. Erwin DH (2007) Disparity: morphological pattern and developmental context. Palaeontology 2007;50: 57–73.

9. Olson ME (2012) The developmental renaissance in adaptationism. Trends Ecol Evol. 27; 2012: 278–287.

10. Armbruster WS, Pélabonan C, Bolstad GH, Hansen TF (2014) Integrated phenotypes: understanding trait covariation in plants and animals. Phil Trans R Soc B 2014;369: 20130245. pmid:25002693

11. Sidlauskas BL. Continuous and arrested morphological diversification in sister clades of characiform fishes: a phylomorphospace approach. Evolution 2008;62: 3135–3156. pmid:18786183

12. Maddison WP. Squared-change parsimony reconstructions of ancestral states for continuous-valued characters on a phylogenetic tree. Syst Zool. 1991;40: 304–314.

13. Betancur R R, Broughton RE, Wiley EO, Carpenter K, Lopez JA, Li C, et al. The tree of life and a new classification of bony fishes. PLOS Currents Tree of Life. 2013;Edition 1. pmid:23653398

14. Gregory WK. Fish skulls: A study of the evolution of natural mechanisms. Malibar, Florida: American Philosophical Society. 1933.

15. Kimmel CB, DeLaurier A, Ullmann B, Dowd J, McFadden M. (2010) Modes of developmental outgrowth and shaping of a craniofacial bone in zebrafish. PLOS One 2010;5: e9475. pmid:20221441

16. Huycke TR, Eames BF, Kimmel CB. Hedgehog-dependent proliferation drives modular growth during morphogenesis of a dermal bone. Development 2012;139: 2371–2380. pmid:22627283

17. DeLaurier A, Huycke TR, Nichols JT, Swartz ME, Larsen A, Walker C, Dowd J, Pan L, Moens CB, Kimmel CB. Role of mef2ca in developmental buffering of the zebrafish larval hyoid dermal skeleton. Dev Biol. 2014;385: 189–199. pmid:24269905

18. Nichols JT, Blanco-Sánchez B, Brooks EP, Parthasarathy P, Dowd J, Subramanian A, Nachtrab G, Poss KD, Schilling TF, Kimmel CB. Ligament versus bone cell identity in the zebrafish hyoid skeleton is regulated by mef2ca. Development 2016;143: 4430–4440 pmid:27789622

19. Hughes GM. A comparative study of gill ventilation in marine teleosts. J. Exp. Biol.1960;37: 28–45.

20. Anker GC. Morphology and kinetics of the head of the stickleback, Gasterosteus aculeatus. T Zool Soc. Lond. 1974;32: 311–416.

21. Kimmel CB, Cresko WA, Phillips PC, Ullmann B, Currey M, von Hippel F, et al. Independent axes of genetic variation and parallel evolutionary divergence of opercle bone shape in threespine stickleback. Evolution 2012;66: 419–434. pmid:22276538

22. McAllister DE. The evolution of branchiostegals and associated gular, opercular and hyoid bones in teleostome fishes. Bull Nat Hist Mus Canada 1968; 2221: 239 pp.

23. Arratia G, Schultze HF. The urohyal: development and homology within osteichthyes. J Morphol 1974;203; 247–282.

24. Wilson LAB, Colombo M, Sánchez-Villagra MR, Salzburger W. Evolution of opercle shape in cichlid fishes from Lake Tanganyika—adaptive trait interactions in extant and extinct species flocks. Sci Rep. 5, 2015; Article number: 16909. pmid:26584885

25. Stange M, Aguirre-Fernandez G, Cooke RG, Barros T, Salzburger W, Sanchez-Villagra MR. Evolution of opercle bone shape along a macrohabitat gradient: species identification using mtDNA and geometric morphometric analyses in neotropical sea catfishes (Ariidae). Ecol Evo. 2016;6: 1–14. pmid:27547357

26. Friedman M. Explosive morphological diversification of the spiny-finned teleost fishes in the aftermath of the end-Cretaceous extinction. Proc. R. Soc. B. 2010;277: 1675–1683. pmid:20133356

27. Near TJ, Eytan R L, Dornburg A, Eytan RI, Keck BP, Smith WL, et al. Resolution of ray-finned fish phylogeny and timing of diversification. Proc. Natl. Acad. Sci. U.S.A. 2013;110: 12738–12743. pmid:23858462

28. Walker MB, Kimmel CB. A two-color acid-free cartilage and bone stain for zebrafish larvae. Biotech Histochem. 2007;82: 23–28. pmid:17510811

29. Kimmel CB, Aguirre WE, Ullmann B, Currey M Cresko WA. Allometric change accompanies opercular shape evolution in Alaskan threespine sticklebacks. Behavior 2008;14: 669–691.

30. Kimmel CB, Hohenlohe PA, Ullmann B, Currey M, Cresko WA. Developmental dissociation in morphological evolution of the stickleback opercle. Evol Dev. 2012;14: 326–337. pmid:22765204

31. Zelditch ML, Swiderski DL, Sheets HD. Geometric morphometrics for biologists: A primer. 2nd edition. San Diego: Elsevier Academic Press. 2012.

32. Klingenberg CP. MorphoJ: an integrated software package for geometric morphometrics. Mol Ecol Resour 2011;11: 353–357. pmid:21429143

33. Rohlf F. tpsDig2, digitize landmarks and outlines, version 2.12. Department of Ecology and Evolution, State University of New York at Stony Brook. 2008.

34. Maddison WP, Maddison, DR. Mesquite: a modular system for evolutionary analysis. Version 3.04 http://mesquiteproject.org. 2015

35. Adams DC, Otarola-Castillo E. Geomorph: an R package for the collection and analysis of geometric morphometric shape data. Methods Ecol Evol 2013;4: 393–399.

36. Adams DC. Quantifying and comparing phylogenetic evolutionary rates for shape and other high-diemnsional phenotypic data. Syst Biol 2014;63: 166–177. pmid:24335426

37. Denton JSS, Adams DC. A new phylogenetic test for comparing multiple high-dimensional evolutionary rates suggests interplay of evolutionary rates and modularity in lantern fishes (Myctophiformes; Myctophidae. Evolution 2015;69: 2425–2440. pmid:26278586

38. Blomberg SP, Garland T Jr., Ives AR. Testing for phylogenetic signal in comparative data: behavioral traits are more labile. Evolution 2003;57: 717–745. pmid:12778543

39. Adams DC. A generalized K statistic for estimating phylogenetic signal from shape and other high-dimensional multivariate data. Systematic Biology 2014;63: 685–697. pmid:24789073

40. Harmon LJ, Weir JT, Brock CD, Glor RE, Challenger W. GEIGER: investigating evolutionary radiations. Bioinformatics Applications Note Phylogenetics 2008;24: 129–131.

41. Pennell MW, Eastman JM, Slater GJ, Brown JW, Uyeda JC, FitzJohn RG, et al. Geiger v2.0: an expanded suite of methods for fitting macroevolutionary models to phylogenetic trees. Bioinformatics 2014;30: 2216–2218. pmid:24728855

42. Simpson GG. The major features of evolution. New York. Simon and Schustser. 1953.

43. Harmon LJ, Losos JB, Davies TJ, Gillespie RG, Gittleman JL, Jennings WB, et al. Early burst of body size and shape evolution are rare in comparative data. Evolution 2010;64: 2385–2396. pmid:20455932

44. Stayton CT. The definition, recognition, and interpretation of convergent evolution, and two new measures for quantifying and assessing the significance of convergence. Evolution 2015;69: 2140–2153. pmid:26177938

45. Stayton CT. What does convergent evolution mean? The interpretation of convergence and its implications in the search for limits to evolution. Interface Focus 2015;5: 20150039. http://dx.doi.org/10.1098/rsfs.2015.0039. pmid:26640646

46. Revell LJ. Plotting methods for comparative data and phylogenies. http://lukejharmon.github.io/ilhabela/instruction/2015/07/05/plotting-methods/.

47. McGhee GR Jr. The geometry of evolution: Adaptive landscapes and theoretical morphology. Cambridge, UK: Cambridge University Press. 2007.

48. Raff RA. The shape of life. Genes, development, and the evolution of animal form. Chicago: The University of Chicago Press. 1996.

49. Wilkins AS. The evolution of developmental pathways. Sunderland, MA. Sinaur Associates. 2002.

50. Gould SJ, Lewontin R. The spandrels of San Marco and the Panglossian paradigm: a critique of the adaptionist programme. Proc. R. Soc. Lond. B. 1979;205: 581–598 pmid:42062

51. Alberch P. Ontogenesis and morphological diversification. Amer. Zool. 1980;20: 653–667.

52. Wagner GP. Homologues, natural kinds and the evolution of modularity. Am. Zool. 1996;36: 36–43.

53. Wagner G, Altenberg L. Complex adaptations and the evolution of evolvability. Evolution 1996;50: 967–976. pmid:28565291

54. Hallgrímsson B, Jamniczky H, Young NM, Rolian C, Parsons TE, Brougner JC, Marcusio RS. Deciphering the palimpsest: Studying the relationship between morphological integration and phenotypic covariation. Evol. Bio. 2009;36: 355–376.

55. Kimmel CB. Skull developmental modularity: a view from a single bone–or two. J. Appl. Ichthyol. 2014;30: 600–607. pmid:25294950

56. Mitteroecker P, Gunz P, Weber GW, Bookstein FL. Regional dissociated heterochrony in multivariate analysis. Ann. Anat 2004;186: 463–470. pmid:15646279

57. Klingenberg CP. Morphometric integration and modularity in configurations of landmarks: tools for evaluating a priori hypotheses. Evol. Dev. 11: 405–421. pmid:19601974

58. Alberch P. The logic of monsters: evidence for internal constraint in development and evolution. Geobios Memoire Special 1989;12: 21–57.

59. Pigliucci M. Adaptive landscapes, phenotypic space, and the power of metaphors. Biol. Rev. 2008;83: 283–287.

60. Price SA, Schmitz L, Oufiero CE, Eytan RI, Dornburg A, Smith WL et al. Two waves of colonization straddling the K–Pg boundary formed the modern reef fish fauna. Proc. R. Soc. B 2014;281: 20140321. http://dx.doi.org/10.1098/rspb.2014.0321. pmid:24695431