Early Gnathostome Phylogeny Revisited: Multiple Method Consensus

A series of recent studies recovered consistent phylogenetic scenarios of jawed vertebrates, such as the paraphyly of placoderms with respect to crown gnathostomes, and antiarchs as the sister group of all other jawed vertebrates. However, some of the phylogenetic relationships within the group have remained controversial, such as the positions of Entelognathus, ptyctodontids, and the Guiyu-lineage that comprises Guiyu, Psarolepis and Achoania. The revision of the dataset in a recent study reveals a modified phylogenetic hypothesis, which shows that some of these phylogenetic conflicts were sourced from a few inadvertent miscodings. The interrelationships of early gnathostomes are addressed based on a combined new dataset with 103 taxa and 335 characters, which is the most comprehensive morphological dataset constructed to date. This dataset is investigated in a phylogenetic context using maximum parsimony (MP), Bayesian inference (BI) and maximum likelihood (ML) approaches in an attempt to explore the consensus and incongruence between the hypotheses of early gnathostome interrelationships recovered from different methods. Our findings consistently corroborate the paraphyly of placoderms, all ‘acanthodians’ as a paraphyletic stem group of chondrichthyans, Entelognathus as a stem gnathostome, and the Guiyu-lineage as stem sarcopterygians. The incongruence using different methods is less significant than the consensus, and mainly relates to the positions of the placoderm Wuttagoonaspis, the stem chondrichthyan Ramirosuarezia, and the stem osteichthyan Lophosteus—the taxa that are either poorly known or highly specialized in character complement. Given that the different performances of each phylogenetic approach, our study provides an empirical case that the multiple phylogenetic analyses of morphological data are mutually complementary rather than redundant.

Chondrichthyans lack dermal bones and are characterised by an endoskeleton of eventually calcified cartilage. They are represented by the modern sharks and rays (Elasmobranchii) and chimaerids (Holocephali). Fossil remains of chondrichthyans, apart from isolated teeth and fin spines, are rare and the anatomy of early chondrichthyans (Early Devonian or older) is known only from one isolated braincase [5] and one nearly complete specimen, Doliodus [6]. Osteichthyans are characterized by endochondral bone within their endoskeleton, and include about 98% of extant vertebrate species. Osteichthyans diverged along two major lineages, namely actinopterygians (bichirs, sturgeons, gars, bowfins and teleosts) and sarcopterygians (coelacanths, lungfishes and tetrapods). ' Acanthodians' are a group of jawed vertebrates with small, square-crowned scales, spines before the dorsal, anal and paired fins, and a heterocercal caudal fin. This exclusively Paleozoic group exhibits a mosaic of shark-and bony fish-like characters that has long given them prominence in discussions of early gnathostome evolution [2,4,[7][8][9][10][11]. Their relationships with modern gnathostomes have remained mysterious, partly because the detailed endoskeletal structure is only known by the latest, highly specialized Acanthodes bronni [7,8,10,[12][13][14][15]. Placoderms, with their characteristic armor of bony plates, were the most successful and diverse group of jawed fishes during the Late Silurian and Devonian. They have an excellent fossil record because their dermal bones were generally robust and easily preserved. Placoderms are of great significance as a model for the ancestral gnathostome condition.
The phylogeny of early gnathostomes is one of the puzzling issues in the study of early vertebrates. Our understanding of early gnathostomes has improved greatly in recent years as a result of new discoveries [16][17][18][19][20] and the re-examinations of available fossils [10,11,21]. Although some areas have remained controversial, such as the interrelationships of placoderms, recent studies recovered consistent phylogenetic scenarios of early gnathostomes, such as the paraphyly of placoderms, and antiarchs as the sister to all other jawed vertebrates [4,9,18,[20][21][22].

Conflicting Phylogenies of Early Gnathostomes
Recently, Long et al. [20] applied maximum parsimony (MP) analysis to the dataset of Dupret et al. [21] with the addition of four extra characters and 14 additional placoderm taxa, but recovered different results from other phylogenies [18,22] in the positions of Entelognathus, ptyctodontids and the Guiyu-lineage (Guiyu, Psarolepis and Achoania) (Fig 1).
Entelognathus was discovered from the Late Silurian in China. It combines typical placoderm characters of dermal skeleton and braincase with osteichthyan-like marginal jaw bones, and has been considered in a polytomy with arthrodires, ptyctodontids and crown gnathostomes or as the sister group of crown gnathostomes [18,21]. Long et al. [20] recovered Entelognathus as a stem osteichthyan. However, they have admitted that "this position on the tree could be an artefact" caused by the absence of dermal bone jaw characters for chondrichthyans and acanthodians. Furthermore, they stated "the braincase and palatoquadrate of this taxon clearly distinguish it from the Osteichthyes" ( [20]: supplementary information). This new position impacts our understanding of various character acquisitions and character polarities related to the origins of gnathostomes, such as the innovation of marginal jaw bones (premaxilla, maxilla and dentary) and operculogular series. For example, under the recently resolved framework in which all acanthodians are placed in the stem group of chondrichthyans [18,20,22], given Entelognathus as the sister group of crown gnathostomes, the marginal jaw bones and operculogular series are possibly present in the common ancestor of chondrichthyans and osteichthyans, and thus might have been secondarily lost in acanthodians and chondrichthyans. However, if we accept Entelognathus as a stem osteichthyan as in Long et al. [20], the marginal jaw bones and operculogular series are regarded as synapomorphies of osteichthyans.
Three members of the Guiyu lineage are all found from the Late Silurian-Early Devonian of South China and northern Vietnam [16,[23][24][25]. Guiyu and Psarolepis manifest a combination of features found in both sarcopterygians and actinopterygians (e.g. pectoral girdle structures, the cheek and operculo-gular bone pattern, and scale articulation) [16,24]. They also reveal a combination of osteichthyan and non-osteichthyan features, including spine-bearing pectoral girdles and spine-bearing median dorsal plates found in non-osteichthyan gnathostomes, as well as cranial morphology and derived macromeric squamation found in crown osteichthyans [16]. They were referred to stem sarcopterygians in most earlier studies [16,26,27]. The phylogenetic analysis in Zhu et al. [24] assigned two possible positions for Psarolepis, either as a stem sarcopterygian or as a stem osteichthyan. Qu et al. [28] reported that the absence of tooth enamel in Psarolepis, contrasting with its presence in both actinoptergians and sarcopterygians. This character distribution corroborates Psarolepis as a stem osteichthyan. The phylogenetic hypothesis proposed by Long et al. [20] focused on placoderms and did not add many actinopterygian and sarcopterygian taxa in their analysis, thus they placed the Guiyu-lineage among stem actinopterygians. This unexpected placement, which was followed by Schultze [29], poses a challenge to attempts at polarizing the plesiomorphic osteichthyan and sarcopterygian characters and reconstructing osteichthyan morphotype.
The position of ptyctodontids is another major discrepancy between the phylogenies of Zhu et al. [18] and Long et al. [20]. Miles and Young [30] proposed ptyctodontids as the sister group of all other placoderms. Goujet [31], and Goujet and Young [32], later set ptyctodontids as the sister group of petalichthyids, and Long et al. [20] placed ptyctodontids within a paraphyletic Petalichthyida. By comparison, Zhu et al. [18] assigned ptyctodontids crownward to the petalichthyid Macropetalichthys, although the relationships of ptyctodontids, arthrodires, Entelognathus and crown gnathostomes were unresolved as a polytomy. To find the sources of these conflicting phylogenies, we re-scrutinized the datasets of Dupret et al. [21] and Long et al. [20], which were in turn modified from Zhu et al. [18], and found that some states in the former were erroneously coded from Zhu et al. [18] by an inadvertent copy (see the Materials and Methods). These confused codings have a significant effect on the phylogenetic results, especially for the positions of Entelognathus and the Guiyu lineage. The clarification of this problem will help to evaluate the phylogeny of Long et al. [20]. Based on the revised dataset of Long et al. [20], we revisited the phylogeny of early gnathostomes and further expanded the dataset from other recent works by Giles et al. [22], Brazeau and de Winter [11] and Lu et al. [33].
For a long time, the MP method had been employed in most phylogenetic studies of paleontological morphological data. Along with their extensive adoption in molecular phylogenetic studies, probability-based methods have been used more often in palaeontological phylogenetic analyses recently [34]. The probability-based methods differ from the MP method mainly in each branch of the recovered tree having a given 'length' that is equally applied to all characters in order to determine the probability of change along that branch [34][35][36]. Due to these different principles, the positions of some taxa vary across the trees recovered from different phylogenetic methods. To evaluate the robustness and confidence of the analysed data, we conducted two types of probabilistic analysis: maximum likelihood (ML) and Bayesian inference (BI) analyses to the expanded dataset, in addition to the MP method.

Materials and Methods
The character data entry and formatting were performed in Mesquite (version 2.5) [37]. All characters were unordered, with the exception of two multistate characters which had been ordered as they were first defined (Characters 294 and 302).

Maximum Parsimony Analysis
The dataset was subjected to the MP analysis in TNT software package [38]. The analyses were conducted using a traditional search strategy, with default settings apart from the following: 10,000 maximum trees in memory and 1,000 replications. Bremer support values were generated in TNT by applying the 'New Traditional Search' using TBR and collecting suboptimal topologies with 1,000 replicates. Bootstrap values were generated in TNT using 1,000 replicates. We used the "unambiguous changes only" option in WinClada [39] to optimize character state changes on the cladograms.

Bayesian Inference Analysis
The BI analyses were conducted using MrBayes 3.1.2 [40,41]. Osteostracans were set as the outgroup, and the coding showing polymorphisms were changed to '?' . The analyses used the Lewis stochastic model [42] for standard datatype with coding set to 'variable' and a gamma parameter to account for rate variation across traits. Four simultaneous analyses were run, each including four independent Markov chains: one cold chain and three incrementally heated chains. The analysis was run for 1×10 7 generations. Samples were taken every 1×10 3 generations, resulting in a total of 1×10 4 samples for each of the parallel analyses. Convergence of all four runs was confirmed by superposition of likelihood and posterior probability traces, as viewed in Tracer [43], and effective sample size >1,000 for every parameter. Standard deviation of split frequencies was approximately 0.005 across the four runs. Convergence of the MCMC was graphically checked by plotting frequencies of splits found in Bayesian run 1 against frequencies found in the other three runs, using the online tool AWTY [44]. The first 1.0×10 3 samples for each run, representing the 'burn-in' period, were discarded. The 50% majority-rule consensus tree was computed for the sampled generations of one analysis.

Maximum Likelihood Analysis
The ML analyses were conducted using RAxML v. 8.2.4 [45], with 1,000 bootstraps followed by a maximum likelihood search. Commands used were '-f a -m MULTIGAMMA -K MK -#1,000' . The 50% majority-rule consensus tree was computed for the 1,000 bootstrap trees. ML analysis used the Lewis stochastic model [42] with a gamma parameter to account for rate variation across traits.

Revision of Long et al.'s [20] Dataset
The following are the revisions of seven characters, which were erroneously coded in Dupret et al. 's [21] dataset by an inadvertent copy from Zhu et al. [18] and then adopted by Long et al. [20]. Other revisions see the supplementary materials online.
Character 155-Shape of parasphenoid denticulated field (Character 240 in [18]): Zhu et al. [18] identified three states for this character. By oversight, Dupret et al. [21] mixed states 1 and 2 of this character as a single state, and the third state (slender, splint-shaped parasphenoid, state 2 in [18]) was not coded. The formulation and state codings of this character were restored from Zhu et al. [18].
Character 157-Resorption and redeposition of odontodes (Character 139 in [18]): Dupret et al. [21] missed the codings of the second state (developed resorption and redeposition of odontodes) in their dataset, thus rendered this character uninformative. The codings of the second state were restored from Zhu et al. [18].
Characters 188-Course of ethmoid commissure (Character 183 in [18]) and Character 211-Parasymphysial plate (Character 209 in [18]): Zhu et al. [18] identified three states for each character. By oversight, Dupret et al. [21] mixed states 0 and 1 of each character as a single state, and rendered their codings confused. The formulation and state codings of these two characters were restored from Zhu et al. [18].

Results and Discussion
Comparison between the Phylogeny based on the Revised Dataset and the Phylogeny of Long et al. [20] The MP analysis used the revised dataset with the exclusion of six placoderms (Sinolepis, Gavinaspis, Sigaspis Parayunnanolepis, Quasipetalichthys, and Diandongpetalichthys) with >75% missing data as in the earlier version of this analysis [20]. If only the seven erroneously-coded characters by inadvertent copy are revised, Entelognathus falls into a polytomy with osteichthyans and chondrichthyans in the strict consensus tree (SCT) and is recovered as the sister taxon of all crown gnathostomes in the 50% majority-rule consensus tree (MCT), consistent with the position as it was first described [18] (S1b Fig). Therefore, the six erroneously-coded characters affect the phylogenetic position of Entelognathus, but the Guiyu-lineage is still close to actinopterygians.
The SCT makes a few departures from the original result of Long et al. [20] (Fig 2b). Brindabellaspis and antiarchs are placed in a polytomy with other gnathostomes. Wuttagoonaspis is resolved as the sister taxon to all other members in a major clade grouping the remaining gnathostomes except Brindabellaspis and antiarchs. By comparison, Wuttagoonaspis is placed in the arthrodire clade in the original SCT (Fig 2b). The grouping of petalichthyids and ptyctodontids is not supported as in the original analysis. Lunaspis and Macropetalichthys fall into a polytomy with the remaining gnathostomes except Brindabellaspis, antiarchs and Wuttagoonaspis.
Entelognathus is recovered as the sister taxon of all crown gnathostomes, consistent with the position as it was first described [18]. The Guiyu lineage is placed as stem sarcopterygians rather than stem actinopterygians as in the phylogeny of Long et al. [20]. All sampled acanthodians are placed as stem chondrichthyans (contra [9] and [10]), but as a paraphyletic assemblage [18,20]. Dupret et al. [21] placed acanthodians as the monophyletic sister group to chondrichthyans. However, that analysis was incorrect as pointed out by Long et al. [20]. The real results also support the acanthodians as a paraphyletic assemblage ( [20]: extended data figure 7).
The original and revised matrices with all 91 taxa (i.e. without deleting those six taxa with >75% missing data) were also analyzed using the MP method. The SCT ( The SCT places Brindabellaspis and antiarchs in a polytomy with other gnathostomes. Petalichthyids and ptyctodontids form a clade which also includes Wuttagoonaspis as the sister taxon to all other members in this clade. Arthrodires, which were placed as the closest relatives of crown gnathostomes in previous analysis, are sister to the clade uniting the paired taxa of Jagorina and Gemuendina and crown ganthostomes, consistent with the phylogeny of Giles et al. [22]. Entelognathus is recovered as the sister taxon of all crown gnathostomes, as it appeared on the phylogeny recovered from the revised dataset and as it was first described [18]. The relationships of the crown gnathostomes are not well resolved in the SCT, although the osteichthyans clade remain intact. Janusiscus is in a polytomy with Lophosteus, osteichthyans and totalgroup chondrichthyans, different from the result from Giles et al. [22]. The Guiyu lineage is placed in a polytomy with the clade uniting other sarcopterygians and the actinopterygian clade. Most acanthodians collapse into polytomies, while conventionally-defined chondrichthyans remain monophyletic. On the 50% MCT, the Guiyu lineage is resolved as stem sarcopterygians [18] rather than stem actinopterygians [20,29]. Ramirosuarezia is placed in a monophyletic group including the acanthodian taxa Tetanopayrus, Diplacanthus, Gladiobranchus and Rhadinacanthus in the 50% MCT. All acanthodians are placed as stem chondrichthyans (contra [10]), but as a paraphyletic assemblage [18,20].

BI Analysis of the Expanded Dataset
The BI analysis (Fig 4a) places antiarchs as the sister group to all other gnathostomes. Brindabellaspis and Romundina, together with the sister pair of petalichthyids and ptyctodontids, form a series of plesions that are successive sister groups of crown gnathostome node.
Wuttagoonaspis is assigned a position sister to a clade comprising arthrodires Jagorina and Gemuendina, and crown gnathostome. This clade was unresolved as a polytomy.
Entelognathus is close to the crown gnathostome node. Janusiscus is situated in a polytomy with Lophosteus, osteichthyans and total-group chondrichthyans. The Guiyu lineage is resolved as stem sarcopterygians.  Ramirosuarezia is placed as the sister taxon to all other members along the acanthodianchondrichthyan branch. Ramirosuarezia has previously been compared with rhenanid placoderms and holocephalan chondrichthyans [38] and was recovered in a polytomy with Janusiscus and the gnathostome crown [22]. All acanthodians fall into the stem segments of the chondrichthyan total group.

ML Analysis of the Expanded Dataset
The ML analysis (Fig 5) places Brindabellaspis and antiarchs in a polytomy with the main lineage including other gnathostomes. Romundina is the sister taxon to all other members of this main lineage. The sister pair of petalichthyids and ptyctodontids is retained. Wuttagoonaspis is assigned to the arthrodire cluster, as the sister taxon of Groenlandaspis and closely related to the brachythoracids.
Entelognathus is close to the crown gnathostome branch, as in the BI analysis tree. The clade comprising Janusiscus and Ramirosuarezia is sister to the crown gnathostomes. The Guiyu-lineage is resolved as stem sarcopterygians.
Lophosteus is placed as the sister taxon of all other members of total-group osteichthyans. All acanthodians collapse and fall into the stem segments of the chondrichthyan total group.

Consensus of the Three Different Methods (MP, BI and ML)
Based on the expanded dataset, the three phylogenetic methods come to a consensus on the following solutions.
Entelognathus is placed as the sister group of crown gnathostomes. The clade uniting Entelognathus and crown gnathostomes is supported by the Bremer index of 3; this clade is supported by a posterior probability of 1 in the Bayesian tree and bootstrap value of 0.88 in the RAxML analysis. The crown gnathostome clade, excluding Entelognathus, is supported by a Bremer index of 2, a posterior probability of 0.98 in the Bayesian tree and bootstrap value of 0.75 in the RAxML analysis. The clade uniting Entelognathus and crown gnathostomes in one of the MPTs is supported by four homoplastic characters and one unique character (Character 244: metapterygoid portion of palatoquadrate with developed medial ventral protrusion). The reinstatement of Entelognathus as a stem gnathostome supports the hypothesis that marginal jaw bones (premaxilla, maxilla and dentary) and operculogular series are possibly present in the common ancestor of chondrichthyans and osteichthyans, and have been secondarily lost in chondrichthyans.
The Guiyu lineage is placed in the sarcopterygians rather than the actinopterygians. All the MPTs support the placement of Guiyu lineage as stem sarcopterygians. The clade comprising the Guiyu lineage and other sarcopterygians is supported by five homoplastic characters and one unique character (Character 221: unconstricted cranial notochord). The placement Guiyu lineage on the sarcopterygian stem is strongly supported by a posterior probability of 0.99 in the BI analysis and a bootstrap value of 0.60 in the ML analysis.
All acanthodians are a paraphyletic assemblage of stem chondrichthyans in every analysis [18,20] (contra [10]). Dupret et al. [21] placed acanthodians as the monophyletic sister group of chondrichthyans. However, as discussed by Long et al. ([20]: extended data figure 7), the reported tree in Dupret et al. [21] was not the shortest possible for their data dataset. A stem chondrichthyan status of acanthodians would imply that the common ancestor of chondrichthyans and osteichthyans carried a macromeric dermal skeleton, that the macromeric skeleton is homologous in placoderms and osteichthyans, and that it has been replaced by a micromeric skeleton in chondrichthyans [18].
The relationships of petalichthyids and ptyctodontids remain uncertain. This group is supported by two unambiguous characters herein: sensory canals enclosed within dermal bones (Character 15) and no pineal perforation in skull roof (Character 24). The BI analysis results in relatively low support for this group (0.6 posterior probability) and the ML analysis returns a bootstrap value of 0.57. This pairing was found by Long et al. (2015) with 85 taxa included and the 50% MCT resulted from the revised total dataset (91 taxa included) (S5b Fig). Zhu et al. [18] assigned ptyctodontids crownward to the petalichthyid Macropetalichthys. Our analyses suggest that the close relationship of petalichthyids and ptyctodontids remains plausible, although the sensitivity of this pairing to minor changes in the dataset or analysis suggests that the evidence is weak at best.

Incongruence among Different Methods
The main controversy of early gnathostome phylogeny using different methods mainly relates to the positions of the placoderm Wuttagoonaspis, the stem chondrichthyan Ramirosuarezia, and the stem osteichthyan Lophosteus-taxa that are either poorly known or highly mosaic in character complement. The nodes around these taxa are either poorly resolved or weakly supported within each method.
Wuttagoonaspis is resolved as the sister taxon of the clade uniting petalichthyids and ptyctodontids in the MP method. The clade is supported by two unambiguous homoplastic characters: central dermal skull bone (nuchal) with converging but not meeting posterior pit-line canals and supraorbital canals (Character 248); cheek plate divided (Character 278). It is placed as the sister taxon to the clade uniting arthrodires, Jagorina, Gemuendina and other crown gnathostomes in the BI tree. Wuttagoonaspis is assigned in the arthrodire cluster, as the sister taxon of Groenlandaspis and closely related to the brachythoracids in the ML tree. Wuttagoonaspis is an unusual form from Australia [74,75,122], which has been compared to petalichthyids [74]. However, these similarities have been disputed [30] and it is generally agreed that Wuttagoonaspis shares specialized features with arthrodires [74,[122][123][124] but could not be accomodated within a certain arthrodire classification. It is placed as the sister group of phyllolepids ( [123]: three of six the MPTs, figure 1; [125]), or grouped with arctolepids ( [123]: three of six the MPTs, figure 1). Long et al. [20] placed Wuttagoonaspis within the brachythoracid clade. By comparison, Wuttagoonaspis is resolved as the sister taxon of all other members in a major clade grouping most gnathostomes excluding Brinabellaspis and antiarchs in the revised results (Fig 2a). The unstable placement of Wuttagoonaspis is most likely caused by conflicting character sets suggesting different relationships, caused by convergent evolution with other specialised groups. Deletion of this taxon has no impact on the tree topology under the MP method.
Ramirosuarezia is placed in a monophyletic group including acanthodian taxa Tetanopsyrus, Diplacanthus, Gladiobranchus and Rhadinacanthus in 92% of the MPTs (S6 Fig). The other 8% MPTs recover Ramirosuarezia as a stem gnathostome either as the successive plesion to or the sister group of Janusiscus (Fig 6), the latter result also found in Giles et al. [22]. The BI analysis (Fig 4) places Ramirosuarezia as the sister taxon to all other membersin the Early Gnathostome Phylogeny Revisited acanthodian-chondrichthyan branch. The ML analysis (Fig 5) recovers Ramirosuarezia and Janusiscus in a clade, as the sister of crown gnathostomes. Ramirosuarezia is an enigmatic gnathostome with characteristics in common with elasmobranch and holocephalan chondrichthyans, some placoderms and osteichthyans [54]. It was stated that no character clearly supports affinities of Ramirosuarezia to any of the currently known major gnathostome groups [54] and therefore its phylogenetic instability is hardly a surprise. There cannot be said to be much evidence supporting a relationship of Ramirosuarezia to particular acanthodians, as found in the MP analysis, as the character sets that can be scored for these two groups are largely non-overlapping. Placement of Ramirosuarezia within acanthodians is largely due to a single character (presence of pronounced dorsal process on Meckelian bone or cartilage), shared with Diplacanthus, Tetanopsyrus and Gladiobranchus. The deletion of Ramirosuarezia, which has more than 85% missing data, does not affect the tree topology using the MP method. Thus the unstable placement of Ramirosuarezia might be caused by its incompleteness and the difficulty of interpreting its characters against any model.
Lophosteus is situated in a polytomy with Janusiscus, osteichthyans and total-group chondrichthyans in the MP and BI trees. The ML analysis recovers Lophosteus as a member of totalgroup osteichthyans (Fig 5), with a bootstrap value of 68. Lophosteus is known from very limited materials including a fragment of jaw bone, some isolated scales, teeth and minute bone fragments [86][87][88][89][90][91]. Because of the limited knowledge about this genus and its mosaic character combination, the affinity of Lophosteus has long been debated. It was referred to osteichthyans on the basis of its diamond-shaped scales when it was first discovered [126,127]. On the basis of histological information, Burrow [89] later reported its possible affinities with placoderms. Its large fin spines also suggested a relationship to acanthodians [87,88,90]. When the denticle-bearing jaw bone was described, it was noted that Lophosteus may be more closely related to crown-group osteichthyans [91]. The unresolved position of Lophosteus in the MP and BI trees might be due to its large proportion of missing data (93.7%). Although its position on the strict consensus tree is unresolved, all the 2496 MPTs confirm its affinities to stem osteichthyans (Fig 6). If Lophosteus is removed from the expanded dataset, the resolution among the stem osteichthyans becomes better-Dialipina, and Ligulalepis form two plesions that are successive sister groups of crown osteichthyans.

Multiple Phylogenetic Methods for Paleontological Data
The MP method is generally used in phylogenetic analyses of paleontological morphological data. While the BI and ML methods are widely used to analyze molecular data, they are comparatively less explored in exclusively morphological data, paleontological data in particular. This is partly due to the small size of most palaeontological dasets which reduces the robustness of parameter estimates. Recently, the sizes of some paleontological datasets are growing (e.g. [128]), due to the accumulation of fossil materials and technique advances (e.g. high resolution computed tomography). On the other hand, some researchers have followed the 'total evidence' approach [129], and treated all available characters simultaneously in a single dataset or combined morphology-only datasets with DNA datasets (e.g. [129][130][131]). Recently, there has been an increase in the application of these probability-based approaches to paleontological morphological data.
The BI and ML methods adopt evolutionary models to evaluate the probability of character change along branches in order to determine the tree topology [34-36]. However, while the ML method determines the single tree that maximizes the probability of the observed data, and bootstrap values are applied to evaluate the robustness of the recovered clades, the BI method relies on the posterior distribution-the probability of the hypothesis given the data and a specified prior probability, posterior probability being the measure of credibility of the recovered clades. Although MP, BI and ML methods are based on different principles, the overall tree topologies recovered by three methods are often similar [132]. It is also noted that the results of an MP analysis can be reproduced by a ML approach in which each character has its own set of branch lengths [133]. Some systematists therefore argue that the usage of multiple analytical methods is largely redundant [132]. However, this is not found in the study presented here.
Firstly, the consistent tree topologies recovered by different methods suggest that the BI and ML approaches do not perform worse than the MP method. At present, the most widely used model in both BI and MP for estimating phylogeny from discrete dataset is the Mkv model (where k is the number discrete character states, and v refers to likelihood calculations being conditioned on only variable characters being present in the data) proposed by Lewis [42]. This updated model is modified from the Jukes-Cantor (JC69) model of nucleotide sequence evolution and was specifically designed to imitate the step-counting method of MP under a likelihood framework. This model assumes a Markov process for character change, allowing for multiple character-state changes along a single branch. Historically, this model has been argued that even in ideal cases in which the evolutionary model perfectly fits the assumed model of the probabilistic method, there are particular cases or biases that affect the efficiency of both ML [134] and BI [135]. Recently, the BI and ML approaches are proved to outperform parsimony for estimation of phylogeny from discrete morphological data [136], whereas, O'Reilly et al. [137] stated that Bayesian methods outperform parsimony but at the expense of topology resolution. Here, based on the tree topologies recovered by BI, ML and MP, the tree resolution of BI and ML are better than MP method. Therefore, we conduct that the BI and ML approaches could be useful in morphological phylogenetics.
Secondly, the BI and ML methods are statistically well understood and could provide biologists with not only better resolution of trees with which to formulate evolutionary hypotheses [138], but also new ways to study evolution [139]. This includes the estimation of divergence times to produce time-calibrated phylogenies [140], which can be produced using only morphological data [141], and optimization of continuous characters [142]. BI also allows the inclusion of phylogenetic uncertainty in comparative analysis [143]. Most phylogenetic comparative methods do not account for uncertainties in phylogenies and intraspecific variations. Not accounting for these sources of uncertainty leads to false perceptions of precision [144].
Thirdly, there are always conflicts among the results of the MP, ML and BI methods [34,128]. The discrepancies might be invoked as ancillary stability measures. As to our dataset, the discrepancies around Wuttagoonaspis, Ramirosuarezia and Lophosteus require more thorough morphological investigation of these relevant taxa. The BI and ML methods have been developed as alternatives to parsimony reconstruction, and could reveal a surprising amount of uncertainty in ancestral reconstructions [145]. Furthermore, although these methods usually achieve concordant topological results, they may generate discordant inferences of character evolution from the same dataset [146]. This indicates that method selection has profound impacts in evolutionary scenarios and taxonomic arrangements. Given that different performances of each phylogenetic approach, our study provides an empirical case that the multiple phylogenetic analyses of morphological data are mutually complementary rather than redundant.
After all, an accurate data dataset is the foundation of any phylogenetic analysis. It is important to first make a more thorough investigation on the morphology of relevant taxa and their associated character states before any phylogenetic analysis. A long-standing idea holds that autapomorphies should be added to the data set. Müller and Reisz [147] noted that Mk model only worked well when morphological data sets incorporated autapomorphies. However, this can be problematic because it is untestable if all potential autapomorphies for a given taxon have been documented. This is particularly troublesome with fossil data, especially when dealing with taxa that are known from fragmentary material or only known from one individual. Alternatively, Müller and Reisz [147] noted that if a gamma parameter was included to account for rate variation, then Mk worked well with or without autapomorphies.

Conclusions
The revised data from a recent phylogenetic study [20] recover Entelognathus as a stem gnathostome rather than a stem osteichthyan, and the Guiyu lineage as stem sarcopterygians rather than stem actinopterygians.
An expanded dataset with 103 taxa and 335 characters, which is the most comprehensive morphological dataset for early gnathostomes constructed to date, is presented.
Both the parsimony and probability-based studies based on the expanded dataset of early gnathostomes confirm the paraphyly of placoderms with respect to crown gnathostomes, and place all acanthodians as a paraphyletic stem group of the conventionally defined chondrichthyans. However, some discrepancies among different methods, including the placements of Ramirosuarezia, Lophosteus, and Wuttagoonaspis, indicate focus areas for further research in the study of early vertebrates. Many areas of the phylogeny are highly unstable. For example, among placoderms, a number of core groups are clearly monophyletic, incuding brachythoracid arthrodires, antiarchs and ptyctodontids. However, their relationships to each other are highly unstable, analogous to the poor ability of morphological data to resolve the relationships of the orders of mammals [148] or birds [34,128]. Accurate resolution of the true tree of relationships of placoderms may not be possible on current evidence, but Bayesian phylogenetic methods at least allow this uncertainty to be accounted for when the tree is used to study evolutionary process.
A better knowledge of the influence of different phylogenetic analyses is surely needed. The present study shows that the combination of parsimony analysis, Bayesian and maximum likelihood inferences are mutually complementary rather than redundant.