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

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 μm 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).