Shared ecological traits influence shape of the skeleton in flatfishes (Pleuronectiformes)

In the age of phylogenetic comparative methods, evolutionary biologists have been able to explore evolutionary trends in form in unique and extraordinarily diverse groups of animals. Pleuronectiformes, commonly known as flatfishes, is a diverse and specialized order of fishes that have remarkable asymmetry induced by ocular migration and a benthic life style. Although flatfishes are unique from other fishes, species within the group are morphologically diverse. The origin of ocular migration has been a primary focus of research; however, little is known about overall shape diversification among the flatfishes. In this study, we use integrative methods to examine how body shape evolved within the flatfishes. Shape was quantified from X-rays using geometric morphometrics for 389 individuals across 145 species. The most recent and robust phylogeny was overlaid onto the morphospace and phylogenetic signal was calculated to ascertain convergence in the morphospace. In addition, phylogenetic linear models were employed to determine if ecological traits were correlated with shape and if size had an effect on overall body shape. Results revealed that the majority of variation evolved recently, within the past 15–10-million-years in the middle Miocene, and is highly variable within the flatfishes. These changes are best summarized by body depth, jaw length and medial fin length. Dorsal and anal fin length are correlated, which may be due to the unique mode of locomotion used by flatfishes. A phylogenetic linear model and phylomorphospace analysis suggested that several ecological traits are correlated with shape, which indicates an ecological role in the diversification of flatfishes.


INTRODUCTION
Since the publication of Darwin's On the Origin of Species, biologists have sought to identify the evolutionary forces driving form. Form can be influenced by many factors, such as ecological interactions, biomechanical constraints and natural selection (Wake & Larson, 1987;Gould, 2002;Adams & Nistri, 2010). The advent of phylogenetic comparative methods provided a mechanistic means to make hypotheses on how different factors influence the evolution of form (Harmon, 2019). In the past two decades, advancements in phylogenetic comparative methods have enabled evolutionary biologists to explore Nelson, Grande & Wilson, 2016). Likewise, their preferred habitat type is widely variable (Gibson et al., 2007;Eschmeyer & Fricke, 2011). Most species can be found at a depth between 0 and 500 m, but some have been recorded at depths below 1,500 m. Flatfishes prefer a variety of benthic habitats, from the mouths of rivers to beyond the continental shelf. Similarly, preferred sediment type and diets are variable among species, with flatfishes burying themselves in a range of substrates, from mud to rocky sediment and diet preferences ranging from polychaetes to other fishes. The tremendous amount of shape variation and their distinctive asymmetry make flatfishes a unique and interesting group for the analyses of body shape. The processes driving morphological change across the flatfishes are not well understood, and studies examining shape across the phylogeny have yet to be assessed in the group. Previous studies have shown a correlation between shape, phylogeny and ecology. For example, research using cichlid models have examined morphological variation in relation to ecology under an adaptive radiation model. These studies suggested that as ecological niches become available, rapid diversification and convergence in body shape arise (Schluter, 2000;Barluenga et al., 2006;Yoder et al., 2010;Muschick, Indermaur & Salzburger, 2012;Wagner, Harmon & Seehausen, 2012;Elmer et al., 2014;Burress, 2015;Ford et al., 2016;Burress et al., 2018). Phylogenetic hypotheses of flatfishes based on molecular data revealed that flatfishes likely arose through adaptive radiation, yet it is not known if shape is influenced by ecology or other factors (Harrington et al., 2016;Byrne, Chapleau & Aris-Brosou, 2018). By examining shape diversification in the context of genomic hypotheses, we can to observe evolutionary trends in body shape.
The objective of this study was to understand when shape diversification occurred and how ecological traits contributed to shape variation within flatfishes. By incorporating morphological methods in a time-calibrated phylogenetic context, we were able to address how shape diversified across the phylogeny and determine if shape was correlated to diet, the maximum depth zone, climate type, water type and/or sediment type. Using the most robust and recent time-calibrated genomic phylogeny (Byrne, Chapleau & Aris-Brosou, 2018), we generated a phylomorphospace and chronophylomorphospace to visualize shape diversification across the evolutionary history of a representative group of flatfishes (Bookstein, Chernoff & Elder, 1985;Rohlf, 2002;Stone, 2003;Zelditch et al., 2004;Clabaut et al., 2007;Brusatte et al., 2008;Sidlauskas, 2008;Sakamoto & Ruta, 2012). The phylomorphospace allowed us to observe where shape was conserved within a clade and where shape converged. Furthermore, a phylogenetic linear model was employed to determine whether ecological traits correlate to body shape.

Morphometric analyses
A total of 389 individuals representing 145 species within the Pleuronectiformes were radiographed from fish collections at the University of Kansas Natural History Museum and Smithsonian National Museum of Natural History. Twelve families were represented by the following number of species: Achiridae n = 10, Achiropsettidae n = 1, Bothidae n = 23, Citharidae n = 4, Cynoglossidae n = 7, Paralichthyidae n = 23, Pleuronectidae n = 32, Poecilopsettidae n = 3, Psettodidae n = 2, Samaridae n = 2, Scophthalmidae n = 3 and Soleidae n = 24 (Table S1). When possible, representatives with minimal visible damage and the most recent collection dates were chosen to reduce chances of bone degradation. Between one and six adult individuals for each species were radiographed. The sex of specimens is unknown; however, sexual dimorphism is generally correlated to larger body sizes in females and rarely in variation of shape in flatfishes (Gibson et al., 2007). Additionally, shape variation associated with sex was overpowered by the high diversity in shape across a wide range of species. Radiographs were taken using a Thermo Scientific Kevex PXS5-927EA Microfocus X-ray source, with a focal spot of 4 mm at 2 W, on a Varian PaxScan 4030 E with Kodak Lanex Fine Screen scintillator digital panel with a 40 × 28 cm dimension for Smithsonian specimens. Radiographs of specimens at the University of Kansas were taken using a GE Picker X-ray head in a Technology for Industry controller on Kodax Mammography X-ray film. Images were captured using VIVA K.03 Image Acquisition/Control Software. To reduce distortion of the body caused during the preservation process, specimens were flattened using a sheet of acrylic glass and fabric hook-and-loop fastener straps. Images were manipulated in Photoshop to improve clarity of radiographs by adjusting brightness and color levels.
The diversity in cranial morphology within flatfishes made it difficult to determine homologous landmarks across the species included in this study. However, we were able to identify ten landmarks defining the body outline (Fig. S1). These landmarks were chosen based on standard landmarks used in geometric morphometric studies of fishes and the ability to capture the overall outline of the flatfish (Zelditch et al., 2004). Flatfishes also have a curvature of the spine which primarily involves the abdominal vertebrae. As the number of vertebrae change across species, the curvature was captured using a series of landmarks and semi-landmarks. The semi-landmarks were evenly spaced between the vertebral landmarks along the length of the spinal column. Semi-landmarks are not individually homologous, instead they sample points along the homologous curve of the spinal column (Zelditch et al., 2004). Landmarks and the curve were digitized using the software TPSdig 2.16 (Rohlf, 2010) and semi-landmarks were appended to landmark files using tpsUtil (Rohlf, 2010).
Specimens were superimposed using a generalized least squares Procrustes superimposition to remove non-shape related information (translation, orientation and size) using geomorph 3.0.5 in R (Adams & Otárola-Castillo, 2013). The superimposed landmarks were then averaged in the base package in R (R Development Core Team, 2016) for each species.

Morphospace methods
A principal component analysis (PCA) was performed for averaged data in the R package geomorph and principal component backtransformations were generated (Olsen, 2017) to view theoretical shape of the morphospace for 145 species. To explore evolutionary trends of body shape within the flatfishes, a phylomorphospace was generated in the R package geomorph (Adams & Otárola-Castillo, 2013). This method projects the phylogeny onto the multivariate morphospace so the magnitude and direction of shape change can be interpreted in a phylogenetic context (Sidlauskas, 2008). The phylogenomic tree ( Fig. S2) (Byrne, Chapleau & Aris-Brosou, 2018) was downloaded from GitHub (https://github.com/sarisbro) and input into the R environment. Non-corresponding specimens were pruned from the tree in the R package ape (Paradis, Claude & Strimmer, 2004), leaving 98 corresponding species in 15 identifiable clades. Where necessary, clades that separated known families were labeled as 1 and 2 (Paralichthyidae 1 and Paralichthyidae 2). The phylomorphospace was generated by overlaying the pruned tree onto the PCA. To visually understand temporal changes in shape across the phylogeny, time was added to the z-axis to generate a chronophylomorphospace. This method plots reconstructed ancestral shapes in the morphospace and across time based on known relationships (Sakamoto & Ruta, 2012).

Phylogenetic signal, allometry and phylogenetic linear model
Phylogenetic signal was calculated using the K mult method in geomorph (Adams & Otárola-Castillo, 2013). The K mult method is a mathematical generalization of the Kappa statistic (Blomberg, Garland & Ives, 2003), and uses a Brownian motion model to evaluate the degree of phylogenetic signal in a dataset (Adams & Otárola-Castillo, 2013). This is the most appropriate method to use for multivariate data (Adams, 2014;Adams & Collyer, 2018a).
An allometric regression was performed to estimate the effect of centroid size and shape using a simple allometric linear model (shape coordinates~log(size)) and a unique family allometric model (shape coordinates~log(size) × family) with procD.lm in geomorph (Adams & Collyer, 2018a). To test the amount of shape variation affected by size, we calculated the morphological disparity for shape with and without size correction (shape coordinates~log(size) × family vs. shape coordinates~family) in geomorph and preformed a simple linear regression to determine if size significantly affected shape. Additionally, we obtained fin lengths using interlmkdist in geomorph and performed a simple regression to determine if dorsal and anal fin lengths correlated to one another.
Primary ecological data (diet, maximum depth zone, climate type, water type and preferred sediment type) were compiled from FishBase (Froese & Pauly, 2019) and a phylogenetic linear model was performed in geomorph using the procD.pgls and pairwise functions (Adams & Collyer, 2018b). When FishBase suggested multiple ecological variables, the primary variable(s) was used. In the case that FishBase suggested two or more variables were equally primary, the ecological variables were considered independent from other variables to account for the unique ecological rank. For example, if diet type was 40% fish, 40% crustaceans and 20% polychaetes, the assigned diet type would be "fish and crustaceans". Where possible, ecological traits in question were cross referenced to the literature, and species without reliable traits were pruned from the dataset.
Phylogenetic linear models calculate the amount of shape variation and the estimated probability of variation attributed to ecological factors in a linear model to detect relationships between shape and ecological traits. A linear model comparing shape to ecological traits was used (shape coordinates~depth + climate + water type + diet + sediment) and type III (marginal) sums of squares (SS) was calculated. Type III SS was computed as the effect of each variable was evaluated after other factors, which means the order of factors in the linear model does not affect the outcome. However, this method is not appropriate for missing data, so we removed taxa which were missing ecological data.
A pairwise test was implemented to identify which ecological types were different from one another. For all tests, a randomized residual permutation procedure with 1,000 permutations was used. Ecological traits that significantly correlated to shape were plotted to the phylomorphospace to visualize trends in shape and ecology.

Morphospace
The first four principal components (Table S2) account for 85.3% of shape variation and were associated with body depth at the location of the first caudal vertebrae, length of the jaw and origin and insertion points of the caudal and anal fins ( Fig. 2; Fig. S3). Principal component one describes 52.0 ± 6.4% of the variation, and principal component two explains 22.2 ± 4.2% of shape variation. Species that fell towards the negative ends of PC1 and PC2 are very round in body shape, whereas species at the positive ends of PC1 and PC2 have oblong body shapes (Fig. 2). Similarly, species that are on the negative end of PC1 and positive end of PC2 have long jaws and short dorsal and anal fins, whereas species on the positive end of PC1 and negative end of PC2 have short jaws and long dorsal and anal fins (Fig. 2). Principal component three describes 6.4 ± 2.2% of shape variation and principal component four describes 4.7 ± 1.9% of shape variation. Shape change across PC3 and PC4 is similar to PC1 and PC2, yet less extreme. Species that fell to the negative end of PC3 have long jaws and short dorsal and anal fins whereas the positive end shows species with short jaws and short dorsal and anal fins. PC4 shows species which fell to the negative end have long jaws with oblong bodies, and species that fell to the positive end have short jaws with deep bodies (Fig. S3).

Phylomorphospace
Several families within the pleuronectiforms were clustered together in the phylomorphospace; however, there is overlap within space for most families. Psettodidae (maroon), Citharidae (pink) and Poecilopsettidae (light blue), Samaridae (gold), Cynoglossidae (orange) and Soleidae (medium blue) show clustering within the family and separation from other families across the PCA ( Fig. 3; Fig. S4). Conversely, Achiridae (black), Bothidae (medium green), Paralichthyidae 1 (dark purple), Paralichthyidae 2 (dark blue), Pleuronectidae (light purple), Rhombosoleidae (brown) and Schophthalmidae 1 (lime green) are widespread across the PCA and display body shapes and jaw lengths of all types ( Fig. 3; Fig. S4). Furthermore, Bothidae and Rhombosoleidae are widespread across PC3 and PC4, whereas other families cluster closely (Fig. S5). Families with shapes on the extreme ends include Cynoglossidae and Soleidae which cluster to the far positive end of PC1, sharing oblong bodies with small jaws and long fins. Psettodidae has an oblong body with a long jaw on the far positive end of PC2, whereas Achiridae has a round body with a short jaw toward the negative end of PC2 ( Fig. 3; Fig. S4). Clustering of families is more clearly shown in a three-dimensional PCA (Fig. S6).

Chronophylomorphospace
Ancestral states were inferred at the nodes of the genomic tree (Fig. S2) (Byrne, Chapleau & Aris-Brosou, 2018) and time was plotted as axis z ( Fig. 4; Fig. S7). Early divergence led to changes in fin length and jaw shape between approximately 40 and 30 MYA with the lineage leading to Psettodiae becoming slightly more oblong with shorter fins and longer jaws, and the lineages leading to Cynoglossidae and Soleidae becoming oblong with longer fins and shorter jaws. At approximately 30 MYA body shape changed across the PC2 axis with changes in jaw length and body depth. The majority of shape diversification occurred 15-10 MYA, where species experienced changes across PC1 in addition to PC2 to become widespread across the phylomorphospace (Fig. 4; Fig. S7).

Phylogenetic signal allometry and phylogenetic linear model
The observed phylogenetic signal was lower than the expected signal (1.0) under a Brownian motion model at 0.6161 and was significant with a p-value of 0.001. A comparison of allometric models showed that unique family allometries (p = 0.001) are appropriate. Larger species in the families Achiridae, Citharidae, Paralichthyidae, Pleuronectidae, Psettodiae, Rhombosoleidae, Samaridae and Scophthalmidae have deeper bodies, longer jaws and short dorsal and anal fins, whereas larger species in the families Bothidae, Cynoglossidae, Poecilopsettidae and Soleidae have oblong bodies with shorter jaws and long dorsal and anal fins (Fig. 5). A linear regression of morphological disparity with and without size correction showed that size did not significantly affect the Procrustes variance (p = 0.08) which suggests shape should not be corrected for allometry. Additionally, a simple regression showed that the dorsal and anal fin lengths were significantly correlated (p < 2E−16); as the dorsal fin increases in length the anal fin increases (Fig. S8). A phylogenetic linear model was used to determine the relationship between ecology and shape while accounting for phylogeny. Ecological factors that were significant included water type (p = 0.003), diet (p = 0.003), and sediment type (p = 0.023) (Table S3). Pairwise tests determined the body shape of flatfishes which are found on mud sediment types were significantly different (p = 0.003) from flatfishes which are generalists (found on mud, sand, or rock) and the body shapes of flatfishes which preferred sand sediments were also significantly different (p = 0.005) from generalists which can be found on mud, sand and rock sediment types (Table S4). Although water type and diet were significant in the phylogenetic linear model, pairwise tests showed no significant differences between groups (Table S4). Significant ecological traits did not visually cluster across the phylomorphospace (Fig. S9).

DISCUSSION
This study is the first to use a well-supported and robust genomic phylogeny to analyze skeletal body shape in a phylogenetic context of the Pleuronectiforms. Geometric morphometric analyses captured body shape diversification across the order and incorporated phylogenetic hypotheses to reveal that evolutionary history and ecological traits are important for body variation. Shape is highly variable within the flatfishes, and is best summarized by changes in body depth, jaw length and medial fin length. Dorsal and anal fin length are correlated, either both elongated or shortened across the morphospace, which may be associated to the unique mode of locomotion in flatfishes. A phylogenetic linear model showed that shape is correlated to ecological traits which may suggest ecology is driving shape. Finally, the majority of variation evolved recently within the past 15-10-million-years in the middle Miocene. Flatfishes are morphologically diverse and range in shape from oblong to rounded with short to long jaws and short to long dorsal and anal fins. Although flatfishes can have any combination of body depth and jaw length, fin length is restricted; the dorsal fin is always longer than the anal fin, and as the dorsal fin increases in length, the anal fin always increases in length (Fig. S8). A linear model shows that dorsal and anal fin length are correlated which may be related to their unique mode of locomotion. Flatfishes use a Tetraodontiform mode of locomotion, the simultaneous use of dorsal and anal fins, in swimming, walking and burial behaviors (Sfakiotakis, Lane & Davies, 1999;Munroe, 2015;Fox et al., 2018). Fin length may be restricted to optimize movement across the sediment which may not be as effective if the dorsal fin was shorter than the anal fin. However, this remains to be tested and awaits further research. Flatfishes also range in length from 4.5 cm (Tarphops oligolepis) to 2.5 m (Hippoglossus hippoglossus) (Chapleau & Amaoka, 1998). An allometric regression showed that size affects shape although the amount of variation in shape that can be attributed to size is very small. Furthermore, our dataset is biased toward smaller specimens (less than 45 cm in length) as we were restricted the size of the X-ray machine. Interestingly, larger species in the families Achiridae, Citharidae, Paralichthyidae, Pleuronectidae, Psettodiae, Rhombosoleidae, Samaridae and Scophthalmidae have deeper bodies, longer jaws and short dorsal and anal fins, whereas larger species in the families Bothidae, Cynoglossidae, Poecilopsettidae and Soleidae have oblong bodies with shorter jaws and long dorsal and anal fins (Fig. 5). Although more data is needed to clarify patterns within families, we were able to demonstrate that size/shape patterns vary across families within the flatfishes and that the effect of size on shape was negligible.
Furthermore, shape and ecological traits are correlated, suggesting that ecology may drive shape in the flatfishes. This is supported by a small phylogenetic signal and significant relationships for water type, diet, and sediment type indicated by the phylogenetic linear model. Although the phylogenetic linear model was significant for several ecological traits, pairwise tests showed no significate differences between groups with the exception for sediment type. The shape of flatfishes which are found on mud sediment and the shape of flatfishes which are found on sand sediment were significantly different from generalists (prefer all sediments types). This suggests that flatfishes which are sediment specialists are different in shape from flatfishes which are generalists. We are unable to identify specific trends in shape and ecology as there was a lack of clustering by ecological trait in the phylomorphospace (Fig. S9C) and a lack of support for pairwise distances in most ecological types. More robust analyses that focus on ecological traits in flatfishes are needed to address how shape is correlated to ecological traits among species.
Chronophylomorphospace results suggest that the majority of shape variation evolved within the past 15-10-million-years during the middle Miocene ( Fig. 4; Fig. S7). Early divergence led to changes in fin length and jaw shape between approximately 40 and 30 MYA for lineages leading to the Psettodiae and the clade containing Cynoglossidae and Soleidae. The Psettodiae became slightly more oblong with shorter fins and longer jaws, and the ancestor of Cynoglossidae and Soleidae became oblong with longer fins and shorter jaws. The middle Miocene marks a time of decreasing temperatures and is often referred to as the middle Miocene disruption. During this time there was a wave of aquatic extinctions which may have led to a speciation and diversification event in the flatfishes.

CONCLUSIONS
The Pleuronectiformes is a highly diverse order, with variation in shape best summarized by changes in body depth, jaw length and medial fin length and is likely influence by ecological traits. Further research into how dorsal and anal fin length influence locomotion is needed to determine if and why dorsal and anal fin lengths are correlated. Overall, the Pleuronectiformes are incredibly diverse in both shape and ecological traits resulting in a vast geographical and ecological distribution.