An integrative phylogenetic approach for inferring relationships of fossil gobioids (Teleostei: Gobiiformes)

The suborder Gobioidei is among the most diverse groups of vertebrates, comprising about 2310 species. In the fossil record gobioids date back to the early Eocene (c. 50 m.y. ago), and a considerable increase in numbers of described species is evident since the middle Miocene (c. 16 m.y. ago). About 40 skeleton-based gobioid species and > 100 otolith-based species have been described until to date. However, assignment of a fossil gobioid species to specific families has often remained tentative, even if well preserved complete specimens are available. The reasons are that synapomorphies that can be recognized in a fossil skeleton are rare (or absent) and that no phylogenetic framework applicable to gobioid fossils exists. Here we aim to overcome this problem by developing a phylogenetic total evidence framework that is suitable to place a fossil skeleton-based gobioid at family level. Using both literature and newly collected data we assembled a morphological character matrix (48 characters) for 29 extant species, representing all extant gobioid families, and ten fossil gobioid species, and we compiled a multi-gene concatenated alignment (supermatrix; 6271 bp) of published molecular sequence data for the extant species. Bayesian and Maximum Parsimony analyses revealed that our selection of extant species was sufficient to achieve a molecular ‘backbone’ that fully conforms to previous molecular work. Our data revealed that inclusion of all fossil species simultaneously produced very poorly resolved trees, even for some extant taxa. In contrast, addition of a single fossil species to the total evidence data set of the extant species provided new insight in its possible placement at family level, especially in a Bayesian framework. Five out of the ten fossil species were recovered in the same family as had been suggested in previous works based on comparative morphology. The remaining five fossil species had hitherto been left as family incertae sedis. Now, based on our phylogenetic framework, new and mostly well supported hypotheses to which clades they could belong can be presented. We conclude that the total evidence framework presented here will be beneficial for all future work dealing with the phylogenetic placement of a fossil skeleton-based gobioid and thus will help to improve our understanding of the evolutionary history of these fascinating fishes. Moreover, our data highlight that increased sampling of fossil taxa in a total-evidence context is not universally beneficial, as might be expected, but strongly depends on the study group and peculiarities of the morphological data.

The suborder Gobioidei is among the most diverse groups of vertebrates, comprising about 2310 species. In the fossil record gobioids date back to the early Eocene (c. 50 m.y. ago), and a considerable increase in numbers of described species is evident since the middle Miocene (c. 16 m.y. ago). About 40 skeleton-based gobioid species and > 100 otolith-based species have been described until to date. However, assignment of a fossil gobioid species to specific families has often remained tentative, even if well preserved complete specimens are available. The reasons are that synapomorphies that can be recognized in a fossil skeleton are rare (or absent) and that no phylogenetic framework applicable to gobioid fossils exists. Here we aim to overcome this problem by developing a phylogenetic total evidence framework that is suitable to place a fossil skeleton-based gobioid at family level. Using both literature and newly collected data we assembled a morphological character matrix (48 characters) for 29 extant species, representing all extant gobioid families, and ten fossil gobioid species, and we compiled a multi-gene concatenated alignment (supermatrix; 6271 bp) of published molecular sequence data for the extant species. Bayesian and Maximum Parsimony analyses revealed that our selection of extant species was sufficient to achieve a molecular 'backbone' that fully conforms to previous molecular work. Our data revealed that inclusion of all fossil species simultaneously produced very poorly resolved trees, even for some extant taxa. In contrast, addition of a single fossil species to the total evidence data set of the extant species provided new insight in its possible placement at family level, especially in a Bayesian framework. Five out of the ten fossil species were recovered in the same family as had been suggested in previous works based on comparative morphology. The remaining five fossil species had hitherto been left as family incertae sedis. Now, based on our phylogenetic framework, new and mostly well supported hypotheses to which clades they could belong can be presented. We conclude that the total evidence framework a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 taxa of the 5brG. Pezold [27] recognized a distinct head pore configuration as an autapomorphy of the subfamily Gobiinae (now in Gobiidae). Hoese and Gill [28] recognized the family Odontobutidae and used 16 characters to propose interrelationships among Rhyacichthyidae, Odontobutidae, and Gobiidae; their Gobiidae included the Butinae, Eleotridinae and Gobiinae, each of which is now assigned family rank [12,29]. Other authors used morphological data for phylogenetic studies of subgroups of various families, as Larson [30] has done for the Mugilogobius group of the Gobiidae, and Murdy [31] for the Oxudercinae. Most of the more recent works that include morphological characters have focused on smaller groups within Gobioidei (e.g. [29,[32][33][34][35][36]).
The known fossil record of Gobioidei is relatively modest in comparison to their recent diversity and comprises about 40 species that are known from skeletal material. In older studies, these species have been assigned to genera or families by comparing meristic characters, such as counts of vertebrae and fin rays (e.g. [37][38][39][40]). More recent authors have also considered the systematic value of certain osteological features (e.g. shape of the palatine, first dorsal fin-pterygiophore formula, number of branchiostegal rays, configuration of the palatineectopterygoid complex, presence of the entopterygoid), which enabled a tentative systematic placement of some fossil gobioids in relation to extant taxa [41][42][43][44][45][46][47]. However, determining the systematic context of fossil gobioids remains a challenge. Firstly, fossil gobioids apparently include several extinct lineages (e.g. [42,43,45,[48][49][50][51]). Secondly, a phylogenetic analysis of extant gobioids based on morphological characters is hampered by the rarity of synapomorphic characters [28], and becomes even more difficult in fossil gobioids, which usually preserve only skeletal traits and otoliths [11]. Thirdly, all comparative approaches suffer from a limited knowledge of the range of skeletal characters of extant gobioids, which is also due to the sheer number of species (currently 2310, see above) and the inaccessibility of some rare taxa [52,53]. This explains why, with the exception of †Pirskeniidae [11], no phylogenetic analyses have yet been conducted to position fossil gobioids within the phylogeny of extant taxa.
The objective of this study was to place ten selected fossil gobioid species within an integrative ("total evidence") phylogenetic framework including published molecular sequence data from extant species in combination with morphological data for both fossil and extant species. The choice of fossils was largely based on our previous works (see section Material and Methods). The overall idea was that placement of fossil gobioids in the context of a rigorous phylogenetic analysis would greatly help to improve our understanding of the evolutionary history of these fascinating fishes. Bannikov and Carnevale 2016 [64] n/a n/a n/a n/a n/a Fossil †Eleogobius brevis (Agassiz,  1839) NHMUK PV OR 42779-80; Gierl and Reichenbacher 2015 [42] n/a n/a n/a n/a n/a Fossil †Eleogobius gaudanti Gierl and Reichenbacher, 2015 Gierl and Reichenbacher 2015 [65] n/a n/a n/a n/a n/a Fossil †"Gobius" francofurtanus Koken, 1891 Weiler 1963; Gierl 2012 [62,66] n/a n/a n/a n/a n/a Fossil †Gobius jarosi Prikryl and Reichenbacher, 2018 Reichenbacher et al. 2018 [46] n/a n/a n/a n/a n/a Fossil †Lepidocottus aries (Agassiz, 1839) Gierl et al. 2013 [44] n/a n/a n/a n/a n/a Fossil †Paralates bleicheri Sauvage, 1883 Gierl and Reichenbacher 2017 [43] n/a n/a n/a n/a n/a Fossil †Paralates chapelcorneri Gierl and Reichenbacher, 2017 Gierl and Reichenbacher 2017 [43] n/a n/a n/a n/a n/a Fossil †Pirskenius diatomaceus Obrhelová, 1961 Obrhelová 1961; Přikryl 2014; Reichenbacher et al. 2020 [11,49,67] n/a n/a n/a n/a n/a Fossil †Pirskenius radoni Přikryl 2014 Přikryl 2014; Reichenbacher et al. 2020 [11,67] n/a n/a n/a n/a n/a Gobiidae, was sufficient to achieve a molecular 'backbone' as its phylogenetic analysis produced a well-resolved tree that fully agrees with published hypotheses (see Results). Ten fossil gobioid species were added to the extant taxon set (Table 1). We selected those that we had examined in previous works [11, 42-44, 46, 62], and added also the oldest known putative gobioid so far, †Carlomonnius quasigobius (Table 1). A short overview of all fossil species is provided in S1 Appendix.

Study and compilation of morphological characters
Literature data, when available, were used to compile phylogenetically informative morphological characters for the extant species (see Table 1 for details). Additionally, X-ray images were produced for the extant species with a Faxitron Ultrafocus facility at the ZSM and examined to determine numbers of vertebrae, fin elements and pterygiophores, and configuration of the caudal skeleton. After radiography, otoliths were extracted from the same specimens and prepared for scanning electron microscopy (SEM) imaging (using a HITACHI SU 5000 Schottky FE-SEM at the Department of Earth-and Environmental Sciences, LMU Munich). SEM images served as the basis for the identification of the otolith characters (Fig 2), which are used here for the first time within a phylogenetic matrix. They include (i) the overall otolith shape (six character states: trapezoid/triangular; long rectangular; high rectangular; quadratic; rounded; longish-ovate), (ii) the posterodorsal projection (two states: present; absent), (iii) the sulcus shape (three states: perciform-like, shoe sole, shoe sole/specialized), and (iv) the sulcus shape and position (three states: shoe sole + centred, shoe sole + shifted anteriorly, not shoe sole). Plesiomorphic otolith character states were defined according to the condition seen in the otolith of the outgroup (see S1

PLOS ONE
Inferring relationships of fossil gobioids is depicted in Fig 2. For the fossil species, skeletal and otolith characters were largely compiled from previous works, but some additional characters could also be added (see Results). A total of 48 morphological characters were assembled. Thirty-eight characters concern bony structures of the skeleton, four characters relate to cartilage, membrane, or tendon  configurations, four refer to otolith morphology, two concern the presence and type of ctenii on the scales, and one is a morphometric character (see S1 Table Part B for all characters). Character states were determined according to literature data and our morphological investigations (based on X-rays, SEM images of otoliths) (see S1 Table Part A for details).
All taxa and characters were assembled in Mesquite 3.61 [87]. We used presence/absence coding (1/0) or up to ten states, depending on the character (see S1 Table).

Preparation of the molecular data matrix
Molecular sequence data for extant species were assembled based on previously published sequences of five markers: rDNA (12S rRNA, tRNA-Val, 16S rRNA), cytb, rag1, zic1, and sreb2. The sequences were mainly from the study of Agorreta et al. [12], supplemented by some data from Near et al. [63] and Thacker et al. [10]. Molecular data were downloaded from GenBank, aligned in AliView 1.26 [88] with MUSCLE [89], followed by manual adjustment and exclusion of ambiguous regions where necessary. Individual gene alignments were then concatenated into a supermatrix (6271 bp) with SeaView 5.0.4 [90]. For details and GenBank accession numbers see Table 1. All data matrices are available on figshare (https://figshare. com/s/ed10b9a5ac382f856a20).

Phylogenetic analyses
We conducted phylogenetic analyses for the extant species based on the morphological character matrix, the molecular supermatrix, and based on the combined molecular and morphological (total evidence) datasets. Adding all ten fossils to the total evidence character set resulted in a collapse of the molecular backbone (see Discussion for possible reasons). Therefore, we used a step-wise approach: (i) a single fossil species was added, (ii) two fossil species of the same genus were added, (iii) four fossil species were added. We used the Bayesian Markov Chain Monte Carlo (MCMC) approach implemented in MrBayes versions 3.2.6 and 3.2.7a [91], as well as implied-weights maximum parsimony (IW-MP, [92]) implemented in TNT 1.5 [93] to infer phylogenies.
In MrBayes, we assigned the Mkv+G model [94,95] to the morphological data, and separate GTR+G models [95,96] to each molecular partition. For each analysis, we ran 2 x 4 MCMC chains in parallel for 5 x 10 6 generations. We used the "sump" command in MrBayes as well as Tracer 1.7.1 [97] to check for convergence and discarded the first 10% of samples of each analysis as burn-in before summarizing the remaining samples in 50% majority-rule consensus (MRC) trees including posterior probability (PP) values for all clades.

Skeleton data
The character state of each character of each species (extant and fossil) is provided in the S1 Table; unknown character states were coded with a question mark.
Two characters could be coded for the first time for †Eleogobis brevis based on the specimens NHMUK PV OR 42779 and 42780, deposited in the Natural History Museum in London: the presence of a single anal fin pterygiophore inserting before the haemal spine of the first caudal vertebra (AP = 1), and the position of the penultimate branchiostegal on the ceratohyal. Re-inspection of the type specimens of †Gobius jarosi revealed that also in this species the penultimate branchiostegal is located on the ceratohyal. Finally, a count of three to four anal fin pterygiophores (AP = 3-4) could be determined for †Lepidocottus aries based on reinspection of the specimens used in Gierl et al. [44]. Aphia minuta (Risso, 1810).-Rounded, margins smooth, short projection slightly below level of ostium tip. Sulcus shifted anteriorly, rather shoe-sole shape, cauda strongly reduced, ostium with dorsal lobe and pointed tip.
Cryptocentrus cinctus (Herre, 1936).-High-rectangular, margins smooth, with a projection in the middle of the dorsal margin and a median constriction on anterior and posterior rims. Sulcus centered, shoe-sole shape with small, rounded cauda and a broad ostium with ostial lobe.

Phylogenetic analyses
Phylogenetic relationships of the extant species. Molecular data.-The Bayesian phylogeny based on molecular data (Fig 4A) recovers all families as monophyletic. Notably, although we used a restricted number of species, the tree is completely congruent with the trees published by Agorreta et al. [12] and Thacker et al. [10]: Rhyacichthyidae and Odontobutidae are sister groups, and together they are sister to the rest of the gobioid families. Milyeringidae, Eleotridae, Butidae, and Thalasseleotrididae are successive sister groups to the 5brG clade, which is composed of well-supported Oxudercidae and Gobiidae.
Maximum Parsimony analysis produced a single most parsimonious tree (S1 Fig in S1 File), which is similar to the Bayesian tree, but some nodes have poor bootstrap support (BS). The tree is topologically identical to the Bayesian tree on family-level, only within Gobiidae there are some minor differences concerning the positions of Amblygobius and Asterropteryx; BS for several nodes is rather low compared to the support in the Bayesian tree.
Morphological data.-In the Bayesian phylogeny restricted to the extant species, the 5brG clade (Gobiidae + Oxudercidae) is recovered with maximum support, but the internal structure of the clade is completely unresolved (Fig 4B). The thalasseleotridid species G. radiatus and Th. iota group together and are sister to 5brG. The Thalasseleotrididae + 5brG clade is maximally supported, consistent with molecular phylogenetic results. Also, similar to the established molecular phylogeny, Eleotridae, Butidae, and Milyeringidae are closely related to that clade (1.00). However, only Eleotridae and Milyeringidae are recovered as monophyletic (0.91, 0.82), whereas Butidae are paraphyletic (Fig 4B). Monophyly of Odontobutidae is also not resolved, and in contrast to the molecular phylogeny the odontobutid species are closer to the above assemblage (1.00) instead of being sister to Rhyacichthyidae, which are recovered paraphyletic (Fig 4B).
The Maximum Parsimony analysis recovered one most parsimonious tree (S2 Fig in S1  File), which is overall similar to the Bayesian tree, but many nodes have very poor bootstrap support (BS). Within the 5brG clade (BS = 95%) the reciprocal monophyly of Gobiidae and Oxudercidae is not recovered. Thalasseleotrididae is monophyletic and highly supported as sister to 5brG (99%). As in the Bayesian tree, among the remaining families Eleotridae (71%) and Milyeringidae (52%) are supported as monophyletic.
Total-evidence approach.-The Bayesian phylogeny inferred from the total-evidence dataset (= combined molecular and morphological data) including only extant species is topologically identical to the molecular phylogeny but shows maximum support for Rhyacichthyidae and Rhyacichthyidae + Odontobutidae

PLOS ONE
Inferring relationships of fossil gobioids some clades within 5brG that are weakly supported. Concerning the topology, the only difference is that the positions of the oxudercid genera Eucyclogobius, Pomatoschistus and Chlamydogobius are resolved while they form a polytomy in the Bayesian tree (see S3 Fig vs. S4 Fig in  S1 File).
Phylogenies including extant and all ten fossil species. Morphological data.-In the Bayesian phylogeny, the 5brG clade (Gobiidae + Oxudercidae)-including the two fossil Gobius spp.-is recovered only with moderate support (0.70), and internal relationships are again unresolved (Fig 5A). The rest of the tree is topologically similar to the one based on the morphology of the extant species only (Fig 4B), Thalasseleotrididae are now recovered with slightly higher support (0.80 vs. 0.77). The two †Eleogobius spp. are placed between Thalasseleotrididae and 5brG (0.60), but monophyly of the genus is not resolved. The †Pirskenius and †Paralates spp. form a weakly supported clade (0.64) that is sister to the above assemblage (0.73), but only †Pirskenius (i.e., †Pirskeniidae) is resolved as monophyletic (0.87). The remaining two fossil species, †Lepidocottus aries and †Carlomonnius quasigobius, are clearly more related to Thalasseleotrididae + 5brG than to Rhyacichthyidae and Odontobutidae (0.84), but their exact placements are not resolved.
The Maximum Parsimony analysis using the same data set recovered six most parsimonious trees. The resulting 50% majority rule consensus tree (S5 Fig in S1 File) is overall consistent with the one based on extant species only (S2 Fig in S1 File) but is very poorly resolved on a deeper level; most nodes have BS < 50%. The two fossil Gobius spp. are included in the 5brG in a clade with Gobius niger, Discordipinna and Lesueurigobius. The two †Eleogobius spp. are sister to 5brG but monophyly of the genus remains unresolved. †Carlomonnius quasigobius is placed in a clade containing Eleotridae and Kribia nana (Butidae). †Pirskenius is monophyletic (64%) but its position, as well as that of the remaining three fossil species, is not further resolved.
Total-evidence approach.-The Bayesian phylogenetic analysis of the total-evidence dataset including all 29 extant and the ten fossil species produced a consensus tree with a largely collapsed backbone (Fig 5B) compared to the tree based on only the morphological data set of the same taxa (Fig 5A). It contains numerous polytomies and shows poor support for many deeper nodes. The 5brG clade (Gobiidae + Oxudercidae) forms a polytomy with the Thalasseleotrididae, the Butidae, the clade of †Paralates + †Pirskenius and three further fossil species ( †Eleogobius brevis, †E. gaudanti, †"Gobius" francofurtanus). Nevertheless, several groups can be recovered within this polytomy. Oxudercidae, which were not resolved in the phylogeny based only on morphology (Fig 5A), are now recovered with high support (0.93), as well as some clades within Gobiidae that correspond to the molecular phylogeny, i.e. the Aphia (0.89) and Glossogobius (0.72) clades. †Gobius jarosi is weakly (0.55) placed as sister to the extant G. niger. The clade of †Pirskenius + †Paralates (see above) is recovered with slightly weaker support (0.59 vs. 0.64), while support for monophyly of †Pirskeniidae is slightly increased (0.93 vs. 0.87). †Lepidocottus aries and †Carlomonnius quasigobius are resolved as members of Butidae, albeit with low support (0.61).
The consensus tree of the two most parsimonious trees of the same data set (S6 Fig in S1 File) is even less well resolved, with very weak BS for many nodes. The 5brG clade is recovered, with reciprocally monophyletic Oxudercidae and Gobiidae. In contrast to the Bayesian tree, †"Gobius" francofurtanus is sister to Gobius niger (< 50%) (vs. not resolved in the Bayesian tree), whereas †Gobius jarosi is sister to Lesueurigobius (also < 50%) (vs. sister to G. niger).
†Eleogobius is monophyletic and sister to the 5brG clade, but again BS for this is negligible. †Pirskenius is monophyletic (65%) but its position, as well as that of the remaining four fossil species, is not further resolved.
Total-evidence phylogenies including one to four fossil species. In the following analyses, a single fossil was added to the data set of the extant species and the trees were inferred based on the morphological and molecular data (total-evidence approach).
†Carlomonnius quasigobius.-The Bayesian phylogeny is topologically identical to the molecular and the total evidence phylogenies based on extant species only; †Carlomonnius

PLOS ONE
Inferring relationships of fossil gobioids quasigobius is placed within Butidae (Fig 6A). Support values are high (>0.90) for most clades, but not for Butidae (0.74 vs. 1.00 in the molecular tree). The single Maximum Parsimony tree retained recovered all recent families, like the Bayesian tree ( S7 Fig in S1 File). However, in this tree †C. quasigobius is placed as sister to Thalasseleotrididae + 5brG, albeit with very low support (< 50%).
†Lepidocottus aries.-The Bayesian phylogeny (Fig 6B) and the single Maximum Parsimony tree retained (S8 Fig in S1 File) reveal the same topology and almost the same support values as described above for the trees including †C. quasigobius. †Lepidocottus aries is placed within Butidae in both the Bayesian and the Maximum Parsimony tree, with high support values in the former (0.95), but very low support in the latter (< 50%).
†Gobius jarosi and †"Gobius" francofurtanus.-In the Bayesian tree (Fig 7A), topologies and support values are similar as described above. In comparison to the Bayesian total-evidence phylogeny based on extant taxa only (S3 Fig in S1 File), somewhat decreased support values concern two internal gobiid clades: the Aphia-lineage (0.95 vs. 1.00) and the clade containing Glossogobius and two members of the Cryptocentrus-lineage (0.90 vs. 1.00). The two fossil species are recovered as successive sister groups to G. niger with moderate support (0.77) ( Fig  5D). In the consensus tree of two Maximum Parsimony trees (S9 Fig in S1 File), †"Gobius" francofurtanus is sister to Gobius niger while the position of †G. jarosi is unresolved in a clade with Lesueurigobius sanzi, Tigrigobius multifasciatus, Asterropteryx semipunctata, Amblygobius phalaena, and Aphia minuta (S9 Fig in S1 File). When only one of the two fossil species is included, a sister-relation with G. niger is apparent in each case, in both the Bayesian and the Maximum Parsimony tree (S10-S13 Figs in S1 File).
†Eleogobius brevis and †E. gaudanti.-The Bayesian phylogeny inferred from the total-evidence dataset including all extant species and either †E. brevis or †E. gaudanti (S14, S15 Figs in S1 File) is topologically identical to the molecular phylogeny, and almost no decrease of support values is seen. †Eleogobius brevis is recovered as sister to Gobius niger with good support (0.87), while †E. gaudanti shows a well-supported (0.87) sister relation to the Thalasseleotrididae. When both species of †Eleogobius are added, they are recovered in a polytomy with the Thalasseleotrididae (0.88) (Fig 7B), and there is slightly decreased support for the 5brG clade (0.90 vs. 1.00), and also for some of the gobiid clades (e.g. Tigrigobius + G. niger: 0.91 vs. 1.00) compared to the total evidence tree using only the extant species (S3 Fig in S1 File).
In the Maximum Parsimony trees, when only one of the †Eleogobius spp. is included, its position matches that recovered in the Bayesian analyses, albeit with poor support (S16, S17 Figs in S1 File). When both species are included, they form a clade and are sister to the Thalasseleotrididae (< 50%, S18 Fig in S1 File).
†Pirskenius diatomaceus and †P. radoni.-In the Bayesian tree, topologies are the same as described above; the two species of †Pirskenius are recovered as sister to Thalasseleotrididae (Fig 8A). However, in comparison to the Bayesian total-evidence phylogeny based on extant taxa only, decreased support values occur for the 5brG clade (0.89 vs. 1.00), the Thalasseleotrididae (0.65 vs. 1.00), and for the gobiid clade of Tigrigobius + Gobius niger (0.90 vs. 1.00). Support for the clade Thalasseleotrididae plus one of the †Pirskenius species is similar when only †P. radoni is added (0.69) (S19 Fig in S1 File), but higher when only †P. diatomaceus is involved (0.91) (S20 Fig in S1 File).
The single Maximum Parsimony tree including both species of †Pirskenius recovers all recent families as monophyletic (S21 Fig in S1 File) and the genus is recovered as sister to Thalasseleotrididae + 5brG; its monophyly is supported with 74% BS. Adding solely †P. radoni or †P. diatomaceus results in the same topology and similar support values (S22, S23 Figs in S1 File) as seen in the tree including both species.

PLOS ONE
Inferring relationships of fossil gobioids  https://doi.org/10.1371/journal.pone.0271121.g008 †Paralates bleicheri and †Pa. chapelcorneri.-Inclusion of both species of †Paralates does not recover a relationship of these two fossil taxa with any of the extant clades in the Bayesian phylogeny, rather they form a polytomy with the Odontobutidae + Rhyacichthyidae clade and the clade containing all other families (Fig 8B). The Bayesian phylogeny that includes only †Pa. bleicheri resolves this species as sister to Odontobutidae (S24 Fig in S1 File), albeit with relatively low support (0.60). In comparison with the total-evidence phylogeny including only extant taxa, support for Odontobutidae decreased (0.61 vs. 1.00), while the high support for all other deeper nodes was not affected. The Bayesian phylogeny that contains only †Pa. chapelcorneri recovers this species as sister to the clade containing the 5brG, Thalasseleotrididae, Butidae and Eleotridae, but with very low support (0.52) (S25 Fig in S1 File; note also decreased support for some backbone nodes in this tree).
In the Maximum Parsimony analyses, inclusion of both species could not resolve their phylogenetic position, as in the Bayesian tree, but here the resolution of the backbone of the tree is even more severely reduced (S26 Fig in S1 File). The Maximum Parsimony results for †Pa. bleicheri (S27 Fig in S1 File) match those of the Bayesian analyses, whereas when only †Pa. chapelcorneri is included, it is placed within Thalasseleotrididae (< 50%, S28 Fig in S1 File).
†Pirskenius spp. and †Paralates spp.-When all four †Paralates and †Pirskenius species were included in the Bayesian analysis, they were recovered in a †Paralates + †Pirskenius clade (0.74), which was resolved as sister to Thalasseleotrididae with moderate support (0.79) (S29 Fig in S1 File). In the Maximum Parsimony tree †Pirskenius is recovered monophyletic (64%) but its position and the positions of the two †Paralates species within a clade together with Butidae, Thalasseleotrididae and 5brG are not resolved; overall, this tree shows very poor resolution at its backbone ( S30 Fig in S1 File).

Discussion
In this study, we have assembled for the first time a dataset comprising molecular and morphological data for Gobioidei that encompasses both extant and fossil species. This approach was necessary as a phylogenetic analysis of the extant species based solely on their morphological characters could only resolve those clades for which morphological apomorphies are known, i.e. 5brG, Thalasseleotrididae, Thalasseleotrididae+5brG, and Eleotridae [see 15,28], while Butidae, Odontobutidae and Rhyacichthyidae each were recovered as paraphyletic. The overall objective was to investigate whether a fossil gobioid species can be confidently placed at family level in the tree of the extant Gobioidei using a Bayesian or Maximum Parsimony totalevidence phylogenetic approach. The results reveal mostly well supported placement at family level when a single fossil species is added to the total evidence data set of the extant species, especially in the Bayesian setting.
Five of the fossil species used here had previously been assigned at family level based on a comparative approach: †Lepidocottus aries had been placed within Butidae [44], †"Gobius" francofurtanus and †Gobius jarosi had been proposed as members of Gobiidae [46,100], †Pirskenius spp. had been placed in its own family †Pirskeniidae [49] and a sister group relation of †Pirskeniidae to Thalasseleotrididae + 5brG has been suggested [11]. Each of those fossil taxa have been recovered in corresponding positions in the present study (Figs 6B, 7A and 8A). This implies that using comparative morphology has been a very appropriate method to classify those fossils. The family assignment of the remaining five fossil species analyzed here ( †Carlomonnius quasigobius, †Eleogobius spp., †Paralates spp.) had been left as incertae sedis in previous work because they possess a mosaic set of characters that is not known among extant gobioids [42,43,64]. †Carlomonnius quasigobius originates from the Eocene of Monte Bolca in northern Italy [64], from the lower to middle Eocene (Ypresian to Lutetian, c. 50-40 Ma, see [101]). It is placed in Butidae in our study (Fig 6A), and it shares with some Butidae (especially with Kribia) a very small size, but any comparative approach would not have assigned †C. quasigobius to this family because it has only five branchiostegal rays (vs. six in Butidae) and a continuous dorsal fin (vs. divided). The feature that seems to have placed †C. quasigobius within the Butidae and close to Kribia is the number of 11 branched and segmented caudal fin rays, which is uncommon among other Gobioidei. Nevertheless, given that †C. quasigobius is the oldest gobioid species currently known [64], an assignment to the Gobiidae (with which it shares the number of five branchiostegal rays) seems unlikely and its classification within the Butidae appears to be more plausible. It would expand the fossil record of Butidae from the early Oligocene (c. 30 Ma, [102]) to the early-middle Eocene (c. 50-40 Ma). However, an additional possibility is that †C. quasigobius is a member of an extinct gobioid family or a "stem gobioid"showing a mixture of derived characters (five branchiostegal rays, dorsal postcleithrum absent, 11 (7+6) segmented and branched caudal-fin rays, four pelvic-fin rays) and plesiomorphic ones (dorsal fin continuous, 24 (10+14) vertebrae, autogenous haemal spine of the second preural centrum, first two abdominal centra shortened, first dorsal-fin pterygiophore inserting in the second interneural space) [64].
In case of †Eleogobius and †Paralates, the resulting phylogenies indicate that these genera are not monophyletic and their species may not even belong to the same family (see Figs 7B and 8B, S14-S18 Figs and S24-S28 Figs in S1 File). The two species of †Eleogobius have been reported from the lower and middle Miocene (c. 17-14 Ma) of Central Europe, specifically southern Germany [42], Austria ( [103], as Gobius), Switzerland [104,105], and Croatia [106]. They have been interpreted as belonging to the same genus because they share a T-shaped palatine, absence of an endopterygoid and presence of six branchiostegal rays, and their otoliths are superficially similar, but show differences to recognize the two species [42]. Of those characters, the T-shaped palatine and absence of an endopterygoid can be considered as apomorphic [23,28], and would support assignment to Gobiidae, which is proposed, based on our phylogenetic results, for †E. brevis (Fig 7A). In contrast, the phylogenetic position of †E. gaudanti near Thalasseleotrididae (Fig 7B) is difficult to understand as no potential synapomorphies are known that can be recognized in a fossil. However, it is more or less consistent with the hypothesis of Gierl and Reichenbacher [42] that †Eleogobius is somewhat "in-between" the 5brG clade and the 6brG gobioids. Furthermore, Bradić-Milinović et al. [41] have recognized a difference in the arrangement of the branchiostegals in the two Eleogobius species, which appears to support the possibility that †Eleogobius is not monophyletic.
In the case of †Paralates, a family assignment had not previously been proposed. †Paralates bleicheri has only been found in lower Oligocene deposits of the southern Upper Rhinegraben [40]. An assignment of †P. bleicheri to the Odontobutidae, as indicated in our analysis, receives little support based on its skeletal traits as all shared characters with the Odontobutidae represent plesiomorphic character states (e.g. seven spines in the first dorsal fin, 8-9 rays in the second dorsal fin, presence of postmaxillary process) (see S1 Table). Finds of fossil skeletons of †P. bleicheri with otoliths preserved in situ would be necessary to reinforce this hypothesis.
†Paralates chapelcorneri originates from the upper Eocene "Chapelcorner Fish bed" of southern England (Isle of Wight) [38,43]. No otoliths have yet been reported from the "Chapelcorner Fish bed". The family assignment of this species remains a topic of future research based on new material of this species. Moreover, a possible relationship between †Paralates spp. and †Pirskenius spp. was indicated in the Bayesian and Maximum Parsimony trees (S25, S26 Figs in S1 File), which is possibly due to their specific combination of plesiomorphic (e.g. seven spines in the first dorsal fin, nine anal-fin rays) and apomorphic traits (e.g. 12 abdominal vertebrae, presence of interneural gap).
Our study also yields some new insights from a methodological point of view. Adding morphological data from only extant species to the molecular dataset had practically no influence on the tree topology and support values (S3 Fig in S1 File). Likewise, the inclusion of a single fossil or of two congeneric species did not change the tree topology, only sometimes some support values decreased (Figs 6-8). However, when all fossils were included in the total evidence phylogenetic framework, the resolution of relationships between families and most fossil taxa dramatically collapsed (Fig 5B). Notably, the morphological phylogeny including the extant and all fossil species was less severely collapsed in the backbone of the tree as the 5brG clade and the Thalasseleotrididae were resolved (Fig 5A). It seems that in the case of the total evidence phylogeny the fossil taxa added a high level of conflicting phylogenetic signals, which could not be overcome by the molecular data despite the latter having orders of magnitude more characters and harbouring strong signal for resolving gobioid phylogeny. A possible explanation is that the fossil taxa do not only add morphological information, but also a lot of question marks to the matrix, because, depending on their preservation, some morphological traits cannot be determined. An additional (or alternative) explanation is that many extinct gobioid clades and families, each with a unique character combination, existed in the past [11,41,45,49]. These cannot be 'forced' in the tree of extant species and eventually may also be responsible for the collapse of the molecular backbone of the extant families. This highlights that increased sampling of fossil taxa in a total-evidence context is not universally beneficial, as might be expected, but strongly depends on the study group and peculiarities of the morphological data.

Conclusion
We have presented a total-evidence dataset comprising molecular and morphological data of 29 extant gobioid species representing all families. Bayesian and Maximum Parsimony analyses revealed that this dataset is sufficient to achieve a molecular 'backbone' that fully conforms to previous molecular work. The new dataset can be used to analyze the family assignment of fossil skeletal-based gobioid species using Bayesian and Maximum Parsimony total-evidence phylogenetic approaches, which has not been possible before.
Our phylogenetic analyses confirmed the family assignment of those fossil gobioid species for which such a placement had been proposed in previous works. It is thus evident that comparative morphology remains an appropriate method to classify some gobioid fossils. However, our phylogenetic analyses could also suggest relationships of fossil gobioid species for cases where the comparative approach did not yield conclusive results. An example is †Carlomonnius quasigobius, which is the oldest gobioid fossil to date and our phylogeny suggests that it could be a possible member of the Butidae, which would expand the known age of fossil butids from the early Oligocene (30 Ma) to the early-middle Eocene (40-50 Ma). Although such positioning of †C. quasigobius remains somewhat speculative for now, it can give hints to look at certain fossil species from a different and new perspective.
We think that the total evidence framework presented here will be beneficial for all future work dealing with the phylogenetic placement of fossil gobioids and thus will help to improve our understanding of the evolutionary history of these fascinating fishes.
Supporting information S1 Appendix. Short description of the 10 fossil specimens included in this study. (DOCX) S1