Introduction

Rates of both phenotypic and species diversification vary widely among animal clades, and a primary goal of comparative biology is to identify and investigate the evolutionary and ecological correlates of these rate differences. Speciation and phenotypic diversification rates may be related to factors such as invasion of novel habitats, habitat complexity, niche specialization, and/or competition. Releases of ecological or competitive constraints can theoretically result in an adaptive radiation, with predicted increases both in species and in phenotypic diversity (Yoder et al. 2010). Among fishes, examples of unusually high species diversification rates have been identified in Lake Baikal sculpins, Rift lake cichlids, Antarctic icefishes, and among larger clades including Ostariophysi and Percomorpha as a whole (Alfaro et al. 2009b; Hunt et al. 1997; Near et al. 2012a, b; Salzburger et al. 2002.; Santini et al. 2009). Similarly, phenotypic diversification has been evaluated in a phylogenetic context in many fish groups (Aguilar-Medrano et al. 2011; Collar et al. 2009; Cooper and Westneat 2009; Price et al. 2011; Sidlauskas 2007, 2008), as well as broadly among acanthomorphs through deeper evolutionary time (Friedman 2010; Sallan and Friedman 2012). Morphological diversification is often assessed by quantifying changes in the shapes of the head and jaws and is thus hypothesized to be linked to feeding adaptations; changes in body shape may also be correlated with novel locomotor patterns. Among the fish and amphibian groups that have been examined, morphological disparity and species diversity have been found to be positively correlated in some cases (Alfaro et al. 2009a; Dornburg et al. 2011; Martin and Wainwright 2011; Price et al. 2010; Rabosky et al. 2013) and negatively correlated in others (Santini et al. 2013; Near et al. 2013); unsurprisingly, evaluation of this relationship depends critically on the parameters used for the evolutionary model (Adams et al. 2009; Rabosky and Adams 2012).

Habitat and ecology are often invoked as causal factors or at least correlates of unusual species diversification in fishes, in particular, switches between marine and freshwater environments. Among pufferfishes, invasion of reef habitats is phylogenetically coincident with higher speciation rates in most clades (order Tetraodontiformes; Alfaro et al. 2007). A finer scale analysis of the family Tetraodontidae indicated an increase in species diversification on coral reefs and a decrease in lineages that returned to freshwater (Santini et al. 2013). Damselfishes have also been the subject of several comparative analyses, with feeding structures shown to have diversified repeatedly, with varying tempos, into a limited set of configurations correlated with feeding strategy (Cooper and Westneat 2009; Frédérich et al. 2013; Litsios et al. 2012). Broader studies of acanthopterygiian diversification have focused on comparing speciation rate with overall measures of shape or body size. They reveal a correlation between size evolution and species diversification among recent fishes (Rabosky et al. 2013), but conversely, a drop in species diversification concordant with elevated shape diversification when patterns are inferred throughout the Cretaceous and Cenozoic (Near et al. 2013). Those studies also separately identified several extraordinarily diverse lineages of fishes: Mormyridae, Ariidae, Mochokidae, Clariidae, Salmonidae, Nomeidae, Cichlidae, Sebastidae + Scorpaenidae, Cottidae, Zoarcidae, Liparidae, and Tetraodontidae in the size-based analyses of Rabosky et al. (2013); Scombriformes, most Blennioidei (Clinidae, Labrisomidae, Dactyloscopidae, Chaenopsidae, and Blenniidae), Afro-American Cichlidae, Liparidae, and gobies [Gobiidae + Gobionellidae sensu Thacker (2009)] in the calibrated phylogeny of Near et al. (2013). No single unique cause can be hypothesized for the elevated diversification in all these disparate groups. As more and finer-scale studies of fish diversification accumulate, a clearer picture emerges of their evolutionary trajectories, featuring intermittent pulses and slowdowns of both speciation and morphological diversification across the phylogeny and in a range of habitats and ecological regimes. Understanding the patterns of morphological change, speciation rates, and ecological shifts through time allows evaluation of broad evolutionary patterns and processes.

One of the most highly diverse groups of recent percomorph fishes, Gobiiformes, particularly Gobioidei, are often cited as an example of exceptional vertebrate species diversity (Birdsong et al. 1988; Nelson 2006; Thacker 2009), and their elevated diversification was confirmed in the broad analysis of teleost diversification of Near et al. (2013). Neither the shape diversity nor the tempo of morphological diversification in gobiiform lineages has been analyzed previously, and it is not known if lineage and morphological diversification are correlated (positively or negatively) or whether or not ecology is correlated with diversification rate shifts. Gobiiformes includes gobies and gudgeons (suborder Gobioidei), cardinalfishes and sweepers (suborder Apogonoidei), and nurseryfishes (suborder Kurtoidei). Ponyfishes (Leiognathidae) may also be close relatives (Thacker 2009). More recent comprehensive phylogenetic studies of fishes have confirmed the close relationships among Gobioidei, Apogonidae, and Kurtidae, but have placed ponyfishes and sweepers separately well apart from Gobiiformes (Near et al. 2012b, 2013; Betancur-R et al. 2013).

Lineages in the order Gobiiformes (as construed here) collectively encompass 2,311 described species distributed worldwide in shallow tropical and temperate waters, including several radiations in freshwater rivers and streams. Among gobiiform families, species diversities vary widely: Gobiidae, the largest family, includes more than a thousand species; Gobionellidae, Apogonidae, and Eleotridae each include one to several hundreds, and the remaining families are much less diverse, with as few as two to four species. These lineages also exhibit great variation in shape, ranging from the rhomboid, highly laterally compressed profiles of Kurtidae to the stout, near-cylindrical aspects of most Gobioidei. Within Gobioidei, shapes vary somewhat among families, but follow a similar general pattern. Apogonidae displays intermediate overall profiles, with laterally compressed, generally rhomboid bodies that vary among lineages primarily in their relative degree of elongation. The genus Pseudamia, one of the more slender forms, has traditionally been classified within Apogonidae, but analysis of molecular phylogeny has indicated that it forms a separate lineage outside the family (Thacker and Roje 2009).

This large radiation encompasses high diversity of species, of form, and of environments inhabited. Gobiiform lineages occupy the full diversity of shelf and continental aquatic habitats, including marine reefs (Gobiidae, Apogonidae), inshore brackish to freshwaters (Gobionellidae, Kurtidae), and exclusively freshwater rivers and streams (Butidae). An earlier study (Thacker 2009) hypothesized a phylogeny for Gobiiformes, noted a correlation between a significantly higher species diversity in Gobiidae as compared to its sister family Gobionellidae, and postulated that this was due to the invasion of fully marine habitats in Gobiidae. Phylogenetic and comparative methods are now available to examine diversification dynamics in detail, including Bayesian time calibration of molecular phylogenies, identification of exceptionally diverse (or depauperate) lineages, and quantification of morphological rates of change. In this study, I utilize these methods to quantify the tempos of lineage and shape diversification and to determine whether or not diversification shifts are correlated with habitat shifts or shape changes. I first construct a time-calibrated phylogeny for Gobiiformes, using fossil calibrations, and test for unusual diversification patterns among lineages. I then analyze shape change among lineages using geometric morphometrics and evaluate the patterns of shape change in context of the phylogeny qualitatively by constructing a phylomorphospace (Sidlauskas 2008; Stone 2003) and quantitatively in two ways. First, using a censored rate test in BROWNIE (O’Meara et al. 2006), I tests for the presence of significant rate shifts between clades specified a priori. Then, I use a Bayesian reversible-jump Markov Chain Monte Carlo (MCMC) approach as implemented in AUTEUR (Eastman et al. 2011) test for rate shifts occurring at any branch across a phylogeny. I then seek out correlations between shifts in species and morphological diversification rates and consider these patterns in the context of broad-scale ecological attributes, with the aim of inferring a comprehensive portrait of the tempo and conditions of gobiiform diversification.

Materials and methods

Phylogenetic reconstruction and calibration

To construct a phylogeny of gobiiform species, I assembled sequence data (mitochondrial cytb, COI, ND1, and ND2) from the studies of Thacker (2009) and Thacker and Roje (2009). I used Mesquite v.2.7.5 (Maddison and Madison 2011) to assemble the DNA sequence matrix for phylogenetic analysis as well as to visualize the phylogenetic hypothesis as a phylomorphospace. I inferred phylogeny using TREEFINDER (Jobb et al. 2004), with maximum likelihood optimization and the dataset partitioned by gene fragment; I also determined appropriate models for each fragment with TREEFINDER, yielding J2+G for COI, TVM+G for ND1, J3+G for ND2, and GTR+G for cytb.

To calibrate the likelihood phylogeny as required for the diversification analyses, I used BEAST v.1.7.5 (Drummond et al. 2012) run with an uncorrelated lognormal-relaxed clock model and a birth/death speciation prior. I assigned divergence time calibrations for Gobioidei, Butidae, Gobionellidae, and Apogonidae based on fossil dates. Two of the fossil calibrations (Gobionellidae and Gobioidei) were based on isolated otolith fossils. For Gobionellidae, the calibration derives from mid-Oligocene otoliths of the genus Pomatoschistus (Reichenbacher and Schwarz 1997; Thacker 2013). The oldest known gobioid fossils are otoliths from an early Eocene deposit in a lignite mine in western India (Bajpai and Kapur 2004). These otoliths were originally described as gobiids, but their morphology more closely matches otoliths of Eleotridae and Odontobutidae (Gierl et al. 2013), so I used them as an age approximation for the base of Gobioidei (Odontobutidae was the most basal taxon considered in this analysis). I based the calibrations for Butidae and Apogonidae (Bannikov 2008; Gierl et al. 2013) on remarkably preserved whole body fossils from the upper Oligocene (Butidae) and middle Eocene (Apogonidae). Calibrations were imposed as exponentially distributed priors, with offset at the calibration age, and conservative means of 10.0. I constrained the analysis with the major nodes present in the likelihood topology (families and lineages) and ran the search for 10.0 × 107 generations, with trees sampled every 1,000 generations. At the end of the analysis, estimated effective sample sizes (ESS) for all parameters exceeded 200; appropriate burn-in fraction was determined using Tracer v.1.5 (Drummond et al. 2012). I constructed a maximum clade credibility consensus of the post burn-in trees using TreeAnnotator v.1.7.5 (Drummond et al. 2012) and combined this time-calibrated phylogeny with taxonomic richness data for each family, derived from FishBase (Froese and Pauly 2014), and the Catalog of Fishes (Eschmeyer 2013) for the species diversification analysis.

Morphological and ecological analyses

To assess the overall shape diversity among gobiiform lineages, I gathered morphometric landmark data from radiographs or cleared and stained preparations of 400 adult specimens representing 86 species of 75 genera in the families Apogonidae (including Pseudamia), Butidae, Eleotridae, Gobiidae, Gobionellidae, Kurtidae, Leiognathidae, Odontobutidae, and Pempheridae (Rhyacichthyidae and Milyeringidae were not used due to the difficulty of obtaining undistorted radiographs in lateral view). All specimens examined were preserved in 10 % buffered formalin and stored in ethanol. For each species, I sampled between one and ten individuals, with variable sampling due to rarity of specimens, difficulty in obtaining undistorted radiographs, or small size. For the smallest bodied species, for which clear radiographs could not be obtained, I cleared and stained specimens following a modification of the method of Pothoff (1984) and photographed them using a Samsung NX10 digital camera mounted on a Zeiss SteREO Discovery v8 dissecting microscope.

I used a suite of 21 landmarks chosen to represent the ranges in overall gobiiform body shape; these are shown on a representative individual in Fig. 1. Landmarks were chosen in order to capture the overall dimensions of shape change, as well as being comparable and consistent across species. All landmarks are skeletally based, such that they may be unambiguously identified on both radiographs and cleared and stained specimens. I digitized landmarks from the right side of each specimen unless that side was damaged or impossible to image, in which case, I imaged the left side and inverted the image. In a few additional cases, I used previously existing radiographs that also required inversion. I used adult specimens only to avoid introducing possible confounding effects due to allometry. I only used individuals that were not distorted or bent, and wherever possible, I imaged multiple specimens. Species examined for morphometric analysis are shown in Table 1, including whether a radiograph or photo of a cleared and stained specimen was used and listed with catalog numbers in “Material Examined.” For the most diverse families, particularly Apogonidae, Eleotridae, Gobiidae, and Gobionellidae, I sampled a relatively small fraction of overall species diversity. However, despite the shape variation that is present among gobiiform families, shape variation within each family is much more restricted, and the overall diversity is well-represented among the sampled species. One notable exception is the lack of extremely elongate anguilliform morphologies which are found in some Gobiidae (Microdesmus) and Gobionellidae (Gobioides). These species could not be used primarily because it was impossible to generate non-distorted radiographs due to deformations of specimens from fixation and storage. I also selected taxa based on whether or not DNA sequence data were available in order to provide an independent data source for inference of phylogeny, restricting the species diversity available for morphometric examination. In most cases, I imaged and morphometrically quantified the same species as was present in the molecular phylogeny, but if that species was unavailable, I substituted a similar congener (Table 1).

Fig. 1
figure 1

Locations of morphometric landmarks. Locations of landmarks used in this study shown on right lateral view of Cheilodipterus isostigmus (USNM 171260). Landmarks digitized are: (1) tip of ascending process of premaxilla, (2) anterior tip of premaxilla, (3) dorsal tip of maxilla, (4) anterior tip of dentary, (5) posteroventral tip of maxilla, (6) posteroventral tip of dentary, (7) dorsal tip of supraoccipital crest, (8) articulation point of first vertebra with basioccipital, (9) dorsal extent of dorsalmost pectoral fin radial, (10) anterior tip of pelvis, (11) anteroventral tip of cleithrum, (12) posterior extent of branchiostegal rays, (13) posterior tip of pelvis and articulation of pelvic fin rays, (14) pterygiophore of first dorsal spine, (15) pterygiophore of last dorsal spine, (16) pterygiophore of first dorsal ray, (17) pterygiophore of last dorsal ray, (18) pterygiophore of first anal spine, (19) pterygiophore of last anal ray, (20) articularion point of last vertebra with hypurals, (21) anterodorsal extent of hypurals

Full size image
Table 1 Species imaged for morphometrics and sequenced for molecular phylogeny
Full size table

For morphometric analysis, I loaded digitized radiographs or digital photos of cleared and stained specimens for each species into tpsUtil, exported them as .tps files, and then imported them into tpsDig2 (Rohlf 2010). I assigned landmarks to each image and exported coordinates into tpsRelw for construction of a generalized least-square Procrustes consensus configuration for each species in order to correct for individual variation in size, rotation, and translation. I combined species files in tpsUtil and forwarded them to MorphoJ 1.02d (Klingenberg 2011) for calculation of covariance matrices and principal components (PC) analysis of Procrustes coordinates (equivalent to relative warp analysis). Axes that explained more than 5 % of the total variance were used in the morphological diversification rate analysis. I combined the calibrated phylogeny with the PC scores in Mesquite v.2.7.5 and plotted the phylomorphospace using the Rhetenor module, which superimposes phylogeny onto a graph of PC1 vs. PC2.

To examine the distribution of environmental preference, I then used Mesquite v.2.7.5 to assemble a matrix of habitat condition (freshwater, brackish, or marine) for each taxon. I coded characters as discrete (0, 1, or 2), with multistate coding used if the species inhabits more than one of the environments, as indicated in FishBase (Froese and Pauly 2014). Many gobiiform species are found in all three coded habitat types; to accommodate polymorphic character reconstruction, I used parsimony to estimate ancestral conditions of lineages and to visualize the character states on the phylogeny.

Comparative analyses of diversification dynamics

I used phylogenetic comparative methods to evaluate both the rate of morphological diversification and the rate of species diversification among lineages. I compared species diversification rates across the phylogeny using MEDUSA (implemented in R version 3.0.0; R Development Core Team 2011), a utility that estimates the likelihood of given species diversities first under a single rate of diversification, then progressively with rate shifts postulated at various nodes. For each iteration, the likelihood scores are compared using a stepwise Akaike Information Criterion (AIC) procedure continued until scores do not significantly improve with inference of an additional rate shift (Alfaro et al. 2009b; Santini et al. 2009). For the MEDUSA analysis, I trimmed the time-calibrated phylogeny down to a diversity tree, with single representatives of the families Leiognathidae, Kurtidae, Pempheridae, Apogonidae, and the genus Pseudamia. Gobioidei, Odontobutidae, Eleotridae, and Butidae, were each represented by single terminals, and within Gobiidae and Gobionellidae, I selected terminals to represent each of the lineages identified by Thacker and Roje (2011) and Thacker (2013). Thus, terminals for Gobionellidae represented the Mugilogobius, Stenogobius, Periophthalmus, and North Pacific lineages. For Gobiidae, 9 of the 13 lineages of Thacker and Roje (2011) were present: American seven-spined gobies, lagoon gobies, reef shrimp gobies, burrowing paired gobies, coral gobies, tiny banded gobies, crested gobies, wormfishes and dartfishes, and inshore gobies. Representatives of the Mediterranean/Ponto-Caspian/Eastern Atlantic gobies, sanddivers, silt shrimp gobies, and flapheaded gobies were not included, so for the purpose of the diversity analysis, I combined the species richness counts for those lineages with those of their closest relative in the hypothesis of Thacker and Roje (2011). Thus, counts for the Mediterranean/Ponto-Caspian/Eastern Atlantic gobies were combined with the American seven-spined gobies, those for the sanddivers were included with the inshore gobies, those for the silt shrimp gobies were included with the lagoon gobies, and those for the flapheaded gobies were included with the crested gobies. Species richness counts for each lineage were compiled as for the families described above using the generic lists of Thacker and Roje (2011) and Thacker (2013), with species counts for each genus derived from FishBase (Froese and Pauly 2014) and the Catalog of Fishes (Eschmeyer 2013). The lineages used in the MEDUSA analysis, along with their species richnesses, are given in Table 2.

Table 2 Species diversity of lineages used in the diversification rate analysis
Full size table

To test whether rates of morphological change (as assessed with the geometric morphometric PC axes) were uniform across the phylogeny, I performed a censored rate test with BROWNIE version 2.1.2 (O’Meara et al. 2006). This analysis compares the fit, under a Brownian motion model, of a single rate of character evolution across a phylogenetic tree compared with a two-rate model in which a shift occurs on one node. I postulated rate shifts at nodes for the clades indicated by MEDUSA to have experienced a diversification rate increase. Likelihoods of the single rate and two-rate models were compared using a modified AIC corrected for small sample sizes (AICc). Significance for the differences between single rate and two-rate AICc’s were assessed using both the Χ 2 p value and the p value estimated from 1,000 replications of a parametric bootstrap. Additionally, in order to investigate the possibility that rate shifts obtained a more complex pattern on the phylogeny, I used AUTEUR (Eastman et al. 2011 implemented in R version 3.0.0 R Development Core Team 2011) to evaluate and identify rate shifts across the entire phylogeny without specifying candidate rate-shifted clades. AUTEUR uses a Bayesian reversible-jump MCMC approach to optimize rate shifts on the phylogeny; for these analysis, I performed 100,000 randomizations of each permutation test, separately for scores on each of the four significant morphometric axes. Finally, I compared the evolutionary timing of shifts in speciation and morphological diversification to the history of habitat shifts across the phylogeny to detect any possible ecological correlates of evolutionary diversification dynamics.

Results

The phylogeny yielded by the likelihood analysis, partitioned by gene fragment, agreed in most respects with those depicted in previous separate studies (Thacker 2009; Thacker and Roje 2009). The one notable difference between this hypothesis and more comprehensive ones previously published was that Odontobutis obscura was placed within Eleotridae rather than as sister to the remainder of Gobioidei. This positioning is likely due to incomplete data for that taxon, so the position was altered using Mesquite v.2.7.5 to accord with that shown by Thacker (2009). Within Gobionellidae, the relationships among basal taxa were weakly supported in the same pattern as in Thacker (2009). These relationships should most conservatively be depicted as polytomies; however, the diversification rate analyses used (BROWNIE, AUTEUR, and MEDUSA) require a fully resolved input phylogeny, so the positions of Mugilogobius chulae and Acanthogobius flavimanus were adjusted to accord with the groupings identified using both morphological and molecular characters in Thacker (2013). This phylogeny was used in all of the time calibration and comparative analyses and is given in Fig. 2. On this calibrated phylogeny, the age of Gobiiformes is estimated at 77.0 MYA, with a 95 % confidence interval of 65.4–91.2 MYA.

Fig. 2
figure 2

Time-calibrated phylogeny of Gobiiformes. This phylogeny is based on four mitochondrial genes, with major clades and lineages labeled. Age calibrations were assigned for Gobiiformes, Eleotridae, Gobionellidae, and Gobiidae indicated on the phylogeny with asterisks. Bars indicate 95 % highest posterior density (HPD) distributions for the time estimates

Full size image

Specimens coded for morphometric landmark data are listed in Table 1. The majority of morphometric variation (85.3 %) was explained by the first four morphospace axes (PC’s: 50.7 % by PC 1, 21.7 % by PC 2, 7.3 % by PC 3, and 5.6 % by PC4). Each of the remaining PC axes accounted for less than 5 % of the variance. Principal component one describes differences in the overall depth of the body, as well as variation in the length of the head from premaxilla to opercle. This axis separates the gobioid families (Gobiidae, Gobionellidae, Butidae, Eleotridae, and Odontobutidae) from the remainder of Gobiiformes. Changes described by PC two include a localized compression of the anterior portion of the head, as well as shifts in the lengths and relative positions of the median fins, particularly the anal fin; variation on this axis primarily distinguishes Kurtidae, Leiognathidae, and Pempheridae from Apogonidae and Pseudamia. The third PC describes shifts in the dorsal fin points of origin and termination, as well as the length of the caudal peduncle and further distinguishes the compact Kurtidae and Pempheridae from Apogonidae, Leiognathidae, and Gobioidei. Principal component four summarizes small changes in the overall shape of the head, as well as depth of the body and inflection of the caudal region. Plots of principal components (PC1 vs. PC2, PC3, and PC4) and their wireframes are shown in Fig. 3. The gobioid families are separated from the remainder on all plots except for a small overlap between Apogonidae and Gobioidei on PC3. Within Gobioidei, the families overlap each other completely; in particular, Gobiidae and Gobionellidae do not appear to be morphometrically distinguishable at all, occupying overlapping locations in the morphospace. Variation in PC1 accounts for half of the overall shape variation and describes change in the depth of the body. PC2, associated with variation in the extent and position of the median fins and the width of the head, separates Apogonidae and Gobioidei from a cluster that contains Leiognathidae, Kurtidae, and Pempheridae (Fig. 3a). The plot of PC1 vs. PC3 (Fig. 3b) is similar to the PC2 results in that PC3 does not distinguish Apogonidae from Gobioidei. In this comparison, there is also a very slight overlap between the apogonid and gobioid regions of morphospace. The PC1 vs. PC4 results (Fig. 3c) also show overlap among groups on PC4. In all of these plots, the genus Pseudamia falls among the Apogonidae, although at the margin of that family’s morphospace. Pseudamia has traditionally been grouped among the cardinalfishes, but a recent molecular phylogenetic analysis (Thacker and Roje 2009) indicates that it is a distinct lineage outside the family. Here, Pseudamia is shown to be included in the morphospace of Apogonidae. The phylomorphospace plot (Fig. 4) most clearly reveals the evolutionary pattern in shape space. The inferred ancestral shape is central in the figure, with Apogonidae, Gobioidei, and a cluster consisting of Kurtidae, Leiognathidae, and Pempheridae each invading distinct portions of overall morphospace.

Fig. 3
figure 3

Graphs of principal component axes for gobiiform families. Points are color-coded by taxon, as depicted in legend at right. APO Apogonidae, BUT Butidae, ELE Eleotridae, GBN Gobionellidae, GOB Gobiidae, KUR Kurtidae, LEI Leiognathidae, ODO Odontobutidae, PEM Pempheridae, PSE Pseudamia. Wireframes describe the shape changes for each principal component, PC1 on the X-axis, with PC2–PC4 on the Y-axis. Light blue wireframe indicates the average specimen shape, and dark blue wireframe shows deformation of 0.1 units on that axis. a PC1 vs. PC2, b PC1 vs. PC3, c PC1 vs. PC

Full size image
Fig. 4
figure 4

Phylomorphospace for Gobiiformes. The phylogeny is superimposed on a plot of PC1 vs. PC2, with points color coded as indicated on the figure. Pseudamia is labeled to distinguish it from the other Apogonidae because the phylogeny indicates that it forms a lineage outside the other apogonids (Thacker and Roje 2009). Arrow indicates the root of the tree as depicted in Fig. 2. Line drawings of fishes are modified from Nelson (2006)

Full size image

The diversity tree for Gobiiformes is shown in Fig. 5. It shows the relationships among major gobiiform lineages, with net diversification rate (r = λ − μ) and ∆AIC listed for each rate-shifted node, as well as the background rate. MEDUSA inferred rate shifts at three nodes: a significant slowdown at the node subtending Kurtidae and accelerations at nodes subtending both Apogonidae and Gobiidae + Gobionellidae. The background net diversification rate of 0.063 lineages/Myr was well exceeded by both Gobiidae + Gobionellidae (r = 0.118) and Apogonidae (r = 0.104). Diversities among the constituent lineages in both Gobiidae and Gobionellidae do not significantly differ; the rate shift was not localized in any particular lineage or group of lineages within the larger clade. In contrast, the estimated net diversification for Kurtidae (r = 0.009) falls well below the background rate, consistent with its markedly low diversity (two extant species).

Fig. 5
figure 5

Time-calibrated diversity tree for Gobiiformes. This tree shows the phylogenetic backbone and collapsed lineages used in the MEDUSA diversification rate analysis. Species diversity for each lineage used is listed in Table 2. Lineages with unusual rates of diversification are indicated in color: blue for Kurtidae, green for Apogonidae, and red for Gobiidae + Gobionellidae. Estimates of net diversification rate and ∆AIC values for each rate-shifted clade are listed in the table

Full size image

BROWNIE was used to assess morphological rate changes in the exceptionally speciose clades Gobiidae + Gobionellidae and Apogonidae. It could not be used to compare the low-diversity lineage Kurtidae because morphometric data was only available for one species, and so a morphological range could not be estimated. BROWNIE indicated significant rate shifts for both Apogonidae and Gobiidae + Gobionellidae on different shape space axes (Table 3). Apogonidae experienced a rate decrease on the PC2 shape axis, which describes the shape of the anterior portion of the head and the positions and lengths of the median fins (∆AIC = 4.508, ∆AICC = 4.303, Χ 2 p value = 0.011, parametric p value = 0.008). In the case of Gobiidae + Gobionellidae, rate of shape change was diminished on two of the four shape space axes. A shift was indicated on PC3, the axis that describes changes in both dorsal fin and caudal peduncle length (∆AIC = 7.081, ∆AICC = 6.877, Χ 2 p value = 0.003, parametric p value = 0.004), as well as PC4, a complex measure encompassing changes in the overall shape of the head, body depth, and caudal inflection (∆AIC = 3.246, ∆AICC = 3.042, Χ 2 p value = 0.022, parametric p value = 0.024). In all three of these cases, the rate shifts are slowdowns relative to the overall rate of diversification, and all were of roughly the same magnitude, with the new rate ranging from 37 to 48 % of the background rate.

Table 3 Rates of morphological diversification in Apogonidae and Gobiidae + Gobionellidae
Full size table

Analysis of phylogenetic patterns of shape change with AUTEUR showed a complex mosaic of shape rate shifts among gobiiforms, with no notable single changes inferred except for a shift in change in PC3 for Pempheridae. The depauperate Kurtidae shows small rate shifts on axes PC2 and PC3, and for both Apogonidae and Gobiidae + Gobionellidae, a variety of shifts are postulated on various branches among the sampled taxa (Fig. 6). Additional rate shifts have also occurred outside Gobiidae + Gobionellidae for all shape axes, such that no single shift in shape change may be localized to correlate with the diversification in speciation. Environmental preferences are optimized on the phylogeny in Fig. 7. Ecological changes from brackish/fresh environments to fully marine are postulated for both Gobiidae and Apogonidae. There is not an unambiguous shift hypothesized at the node subtending Gobiidae + Gobionellidae: that clade as well as its sister Butidae is optimized with an ancestral condition of inhabiting both brackish and freshwater. However, it is notable that Butidae includes predominantly freshwater taxa that may occasionally range into brackish habitats, whereas within Gobionellidae, the habitat preference for various species is generally much broader and more estuarine.

Fig. 6
figure 6

Results of evolutionary rate-shift analyses in AUTEUR. Inferred shifts in rates of shape change are shown optimized across the gobiiform phylogeny. The four trees represent rate-shift analyses for each major principal component axis of shape change: a PC1, b PC2, c PC3, d PC4. In all cases, no single strong rate shift is hypothesized, rather, a mosaic of rate shifts has manifested among all the gobiiform lineages

Full size image
Fig. 7
figure 7

Environmental preference optimized on the phylogeny. The phylogeny is presented with taxa color-coded as they are in Fig. 3, except taxa in Gobiidae, which are rendered in gray. Branches are colored in accordance with optimized environmental preference: blue freshwater, gray brackish/estuarine water, black saltwater

Full size image

Discussion

My results provide the first densely sampled time-calibrated phylogeny of Gobiiformes. My estimation of the age of Gobiiformes at 77MYA (95 % confidence interval range 65.4–91.2 MYA) is comparable to those indicated for Gobiiformes in most previous larger analyses (Betancur-R et al. 2013; Cowman and Bellwood 2011; Near et al. 2012b, 2013; Santini et al. 2009), although the estimates of Chakrabarty et al. (2012) are much older. Those analyses used mitochondrial DNA, nuclear DNA, or both for their inferences of phylogeny, and in most cases did not employ node calibrations within Gobiiformes. The fossil calibration priors used in this study are newly described (Gierl et al. 2013) or newly interpreted in the context of recent phylogenetic work on Gobioidei (Reichenbacher and Schwarz 1997; Bajpai and Kapur 2004; Thacker 2013). Thus, the overall concordance among the majority of these analyses is notable and supports the general validity of the estimates. Gobiiformes are among the older groups of acanthomorphs and relatively early diverging outside the radiation of most acanthomorph fishes (Near et al. 2012b, 2013).

Using a combination of the calibrated phylogenetic analysis and comparative analyses of both species diversity and shape, I show that accelerated speciation in both gobies and cardinalfishes is evolutionarily coincident with a decrease in the rate of shape change manifested differently in the two clades. The comparative analyses of diversity in both species and form indicate that species diversification was amplified relative to background rates in Apogonidae and Gobiidae + Gobionellidae and slowed down in Kurtidae, confirming the results of Near et al. (2013), that show elevated diversification among gobies as compared broadly to other acanthomorphs. Congruent with these diversification rate changes, morphological evolution in both Apogonidae and Gobiidae + Gobionellidae has also shifted, involving decreases in both groups but on different PC axes. Lower rates of diversification are supported on morphometric PC2 for Apogonidae, an axis that describes changes in head width and median fin position. The clade of Gobiidae + Gobionellidae shows decreases on both PC3 and PC4, axes that account for smaller changes in dorsal fin position, head shape, and caudal peduncle shape and inflection. This inverse correlation between species and morphological diversification complements both the results of Santini et al. (2013), in their study of pufferfish lineages, and those of Near et al. (2013) for Acanthomorpha (spiny rayed fishes). Both those studies identified the reverse pattern: drops in speciation rate that were correlated with increases in overall shape or size diversity. An inverse relationship between speciation and morphological diversity is the opposite of what is predicted by both punctuated equilibrium and adaptive radiation theory (Yoder et al. 2010; Rabosky et al. 2013). The pattern may be due to a relatively constant overall rate of morphological change coinciding with, but not linked to, a period of rapid species diversification. In that case, over a given time period, the morphological change that accumulated would be partitioned more finely among more species in a rapidly radiating clade, as compared to one that was diversifying more slowly. It would also suggest that newly diverged sister species would be expected to be morphologically similar. However, the morphological change being investigated at this level, changes in the overall shape or size of the body, are likely not the changes most important in speciation. A pattern of rapid diversification resulting in many species with small distinctions among them would also be consistent with a common phylogenetic pattern, that of poor resolution and weak support at the backbone of a rapidly diversifying clade. This pattern has been observed in Gobiidae, within which distinct lineages are resolved, but the relationships among them are uncertain (Thacker and Roje 2011; Agoretta et al. 2013; Tornabene et al. 2013).

A final consideration is whether or not there is any environmental characteristic that may be correlated with the shifts in speciation and morphological evolution in these groups. Ecologically, cardinalfishes and gobies not only share a preference for nearshore marine habitats, particularly coral reefs, but also mudflats and estuaries. In a phylogenetic context, there is a greater ecological contrast between the primarily marine and estuarine Gobiidae + Gobionellidae vs. their sister taxon, Butidae, a family restricted to freshwaters of Oceania and eastern Asia. The pertinent ecological distinction between Butidae and Gobionellidae may simply be one of degree: although a few Butidae inhabit brackish or marine water, most are exclusively freshwater. In contrast, the vast majority of Gobionellidae are partially or completely salt tolerant and most spend at least part of their life cycle in the ocean. Even the waterfall-climbing stream gobies known from oceanic islands are amphidromous; they spawn in freshwaters but larvae are transported to the ocean, where they mature until returning to freshwater to spawn (Tallebois et al. 2013). The association of higher diversification rate with invasion of marine environments was suggested among gobies by Thacker (2009), but more narrowly, with reference to Gobiidae as compared to Gobionellidae, not their combined clade in the context of Gobiiformes more broadly. Here, I am able to more precisely identify the locations and timing of evolutionary shifts among Gobiiformes, as well as consider how they are distributed relative to morphological and ecological changes, showing that the gradual evolutionary transition to saltwater habitats is linked to elevation of speciation and reduction in morphological diversification. In terms of species diversity, the primarily freshwater Butidae (33 species) is less diverse than the euryhaline, but generally estuarine Gobionellidae (606 species), which is less diverse in turn than the cosmopolitan, almost exclusively marine Gobiidae (1,090 species). This staged evolutionary transition from fresh to salt water is correlated with a staged increase in diversification across these three gobioid clades. The diversification shifts are also congruent with environmental changes and development of the coral reef framework and fauna. The ages inferred for both Apogonidae (51.5MYA) and Gobiidae + Gobionellidae (55MYA), in the early Eocene, correspond roughly to the emergence of both the modern reef fish taxonomic assemblage, as well as the first modern coral reefs (Wood 1999; Bellwood and Wainwright 2002). At that time, an abundance of shallow marine habitats, both reef and non-reef, is likely to have promoted speciation at the West Tethyan hotspot, where reef biodiversity was then concentrated (Renema et al. 2008). Gobies undergoing a transition from exclusively freshwaters into estuarine and marine habitats would have encountered a wealth of new habitats, potentially promoting extensive speciation and giving rise to the remarkably speciose clades found circumglobally today.

Material examined

Numbers in parentheses indicate numbers of specimens examined from that lot, not necessarily total number of specimens in the lot. Species are sorted alphabetically within each family. Institutional abbreviations follow (Leviton et al. 1985).

  • APOGONIDAE

    • Apogon aurolineatus: USNM 179245 (1); USNM 360593 (3)

    • Apogon maculatus: LACM 6837–4 (7)

    • Archamia biguttata: USNM 261086 (3); USNM356414 (4)

    • Archamia fucata: FMNH 4789 (4)

    • Astrapogon puncticulatus: USNM 38471 (1); USNM 195691 (2); USNM 206226 (1)

    • Cercamia eremia: USNM 332334 (1); USNM 334712 (2)

    • Cheilodipterus artus: USNM 319155 (2); USNM 334717 (2); USNM 369953 (3)

    • Cheilodipterus isostigmus: USNM 171260 (3); USNM 276592 (1); USNM 373698 (1)

    • Cheilodipterus macrodon: MNHN 1981–1236 (1); MNHN 1992–0513 (3); LACM 4024 (2); LACM 4249–66 (1); LACM 37388–2 (1); LACM 42463–32 (1)

    • Fowleria aurita: FMNH 75608 (4); FMNH 76694 (4)

    • Fowleria marmorata: USNM 324802 (3); USNM 347168 (2); USNM 347169 (3)

    • Fowleria variegata: FMNH 110391 (2); MNHN 1980–0267 (3)

    • Glossamia aprion: LACM 44628–3 (5); LACM 44633–1 (1)

    • Gymnapogon urospilotus: FMNH uncat. (Westneat collection; 1)

    • Ostorhinchus angustatus: LACM 33723–10 (1); LACM 39986–2 (2)

    • Ostorhinchus cooki: LACM 30859–52 (7)

    • Ostorhinchus holotaenia: USNM 212301 (3); USNM 213305 (1)

    • Phaeoptyx conklini: MNHN 2002–0675 (1); LACM 5380 (7)

    • Pristiapogon exostigma: LACM 42480–16 (2); LACM 51859–9 (2); LACM 42480–16 (3)

    • Pristiapogon kallopterus: LACM 35977–7 (1); LACM 42489–41 (2); LACM 44789–1 (2); LACM 56969–1 (1)

    • Pseudamia gelatinosa: FMNH 75685 (3)

    • Rhabdamia cypselura: FMNH 22278 (1); FMNH 22281 (1); FMNH 22287 (1); FMNH 22288 (1); FMNH 22289 (1); FMNH 22290 (1); FMNH 22291 (1)

    • Siphamia cuneiceps: FMNH 44940 (3)

    • Sphaeramia orbicularis: MNHN 1963–0546 (5)

    • Zoramia fragilis: USNM 205774 (3); USNM 213073 (3); USNM 225698 (4)

  • BUTIDAE

    • Bostrichthys sinensis: FMNH 44908 (2); FMNH 44909 (2)

    • Butis amboinensis: LACM 51860–12 (5)

    • Ophiocara porocephala: LACM 37429–1 (1); LACM 51860–22 (1)

    • Oxyeleotris urophthalma: FMNH 44925 (1); FMNH 63030 (1)

  • ELEOTRIDAE

    • Calumia godeffroyi: LACM 42491–65 (2)

    • Dormitator latifrons: LACM 56197–4 (6)

    • Dormitator maculatus: LACM 2811 (3); LACM 7746 (3)

    • Eleotris fusca: LACM W65-30 (5)

    • Eleotris picta: LACM 2887 (2); LACM 56198–1 (1)

    • Eleotris pisonis: LACM 32461–9 (5)

    • Erotelis armiger: LACM 44426–3 (3)

    • Erotelis smaragdus: LACM31007-21 (4)

    • Gobiomorphus australis: LACM 42621–2 (10)

    • Gobiomorus dormitor: LACM 56192–5 (6); LACM 56195–5 (2)

    • Guavina guavina: LACM 6869–5 (1)

    • Hemieleotris latifasciatus: LACM 4782 (7)

    • Hypseleotris cyprinoides: LACM 51857–3 (10)

    • Leptophilypnus fluviatilis: LACM 56192–7 (1)

    • Mogurnda adspersa: LACM 44630–1 (3); LACM 444637–5 (3)

    • Ophieleotris aporos: LACM 44639–1 (1); LACM 51857–2 (1); LACM 56010–1 (2)

    • Philypnodon grandiceps: LACM 42621–1 (10)

    • Xenisthmus africanus: UMMZ 186061 (2)

  • GOBIIDAE

    • Amblyeleotris wheeleri: LACM 57162–1 (1)

    • Amblygobius phalaena: LACM 56969–2 (2); LACM 541241–1 (1)

    • Asterropteryx semipunctata: LACM 31005–27 (6)

    • Barbulifer ceuthoecus: LACM 6024 (2)

    • Bathygobius lineatus: LACM 21962 (1); LACM 43690–1 (6); LACM 43690–27 (3)

    • Cabillus tongarevae: LACM 54123–002 (1)

    • Coryphopterus dicrus: LACM 2549 (4)

    • Ctenogobiops feroculus: LACM 57164–1 (2)

    • Eviota prasina: UMMZ 186039 (2)

    • Gobiodon citrinus: LACM 42491–69 (1)

    • Lophogobius cyprinoides: LACM 7847 (2)

    • Nemateleotris magnifica: ANSP 130708 (1)

    • Priolepis cinctus: LACM 26566 (3)

    • Ptereleotris zebra: LACM 26559 (9)

    • Risor ruber: LACM 8939–39 (1)

    • Trimma macrophthalma: LACM 33723–54 (1)

    • Valenciennea muralis: LACM 34711–8 (1); LACM 35969–6 (3); LACM 42484–35 (1)

  • GOBIONELLIDAE

    • Acanthogobius flavimanus: LACM 42710–1 (4); LACM 56287–1 (4)

    • Awaous tajasica: LACM 6096 (1); LACM 9196–2 (2); LACM 32358–3 (1)

    • Chaenogobius gulosus: LACM 44905–2 (6)

    • Ctenogobius criniger: LACM 35978–2 (1); LACM 37409–1 (1)

    • Eucyclogobius newberryi: LACM 36657–3 (3)

    • Evorthodus lyricus: LACM 43423–2 (3)

    • Gillichthys mirabilis: LACM 50502–1 (10)

    • Gnatholepis anjerensis: LACM 55560–4 (3)

    • Gobiopterus brachypterus: FMNH 51646 (1); FMNH 51647 (7)

    • Mugilogobius chulae: FMNH 44934 (3)

    • Periophthalmus barbarus: LACM 54540–1 (1)

    • Pseudapocryptes elongatus: FMNH 90909 (6)

    • Scartelaos histiophorus: LACM 38012–2 (2); LACM 38137–1 (2)

    • Stenogobius singeri: FMNH 68459 (10)

    • Stiphodon elegans: LACM 35511–1 (2)

  • KURTIDAE

    • Kurtus gulliveri: TMB01-10-17 (1); TMB01-11-6 (1); TMB01-15-4 (1); TMB01-15-5 (1); TMB0120-10 (1); TMB01-20-11 (1); TMB01-25-3 (1); TMB01-25-3A (1)

  • LEIOGNATHIDAE

    • Gazza squamiventralis: USNM 345525 (1); USNM 345526 (2)

    • Leiognathus equulus: LACM 369 (1); UMMZ 40502 (6)

    • Secutor megalolepis: USNM 346687 (2)

  • ODONTOBUTIDAE

    • Odontobutis obscura: FMNH 55376 (1); FMNH 57304 (2)

  • PEMPHERIDAE

    • Pempheris vanicolensis: FMNH 22609 (1)