Exploring phylogenetic relationships of Pteraspidiformes heterostracans (stem-gnathostomes) using continuous and discrete characters

Ostracoderms are a paraphyletic group of extinct jawless fishes comprising the gnathostome stem and are fundamental to our understanding of early vertebrate evolution. However, only a handful of these clades have robust phylogenies in place, hindering our interpretation of early vertebrate histories. A new phylogeny is proposed for the Pteraspidiformes – the largest and most-studied clade of heterostracan ostracoderms. Difficulties such as large amounts of missing data and the limited morphological variability within the group have led us to explore different coding strategies such as the inclusion of quantitative data and implied weighting. We present a new comprehensive data set including all described genera of Pteraspidiformes (47 taxa) analysed using discrete characters only, a combination of discrete and continuous characters, and gap-coding strategies (transforming the continuous into discrete characters), along with a Bayesian analysis. Two representatives of the Psammosteidae (Drepanaspis and Psammosteus) are also incorporated within the analysis to elucidate their inclusiveness within the Pteraspidiformes. Well-resolved trees are only achieved under re-weighted (implied weighting) analyses. Here, we show that many ‘classic’ Pteraspidiformes clades hold true under our different coding methods, with the implied weighting of discrete characters and inclusion of continuous characters giving very similar topologies. In all instances, the Psammosteidae are found to belong within the Pteraspidiformes, nested with the Spitsbergen genera Doryaspis, Xylaspis and Woodfjordaspis. Gap coding, however, results in a different tree topology to other analyses, perhaps due to the high sensitivity to missing data. Our results indicate that careful consideration and justifications should be applied to quantitative characters when reconstructing relationships of homoplastic ostracoderms. Superfamily Doryaspidae superfam. nov. is introduced. http://zoobank.org/urn:lsid:zoobank.org:pub:BF5255A1-1417-432B-AC64-BF132A0BCC65


Introduction
Phylogenetic relationships of fossil jawless vertebrates (Agnatha) are imperative to our understanding of vertebrate evolution. Bony jawless vertebrates (ostracoderms) constitute the majority of the gnathostome (jawed vertebrate) stem group. Jawed vertebrate novelties such as acellular and cellular bone, a bony dermoskeleton, a differentiated gut, a mineralized neurocranium and paired pectoral fins, arose first within the jawless ostracoderms (i.e. heterostracans, thelodonts, anaspids, galeaspids and osteostracans) (Donoghue & Keating 2014). Identification of ostracoderm evolutionary dynamics along with their ancestral morphotypes will greatly improve our understanding of early vertebrate histories and thus the circumstances behind their diversification and the rise to dominance of gnathostomes.
Heterostraci are taxonomically the largest group of stem-gnathostomes and along with other pteraspidimorphs have been interpreted as the most basal with a bony skeleton (Forey & Janvier 1993;Donoghue et al. 2000;Donoghue & Aldridge 2001;Janvier 2001;Donoghue & Keating 2014). In previous analyses they have been recovered as the sister clade to galaeaspids, osteostracans and jawed vertebrates, or as a paraphyletic group leading to these taxa (Donoghue & Smith 2001;Gess et al. 2006;Sansom et al. 2010;Blom 2012). However, despite being intensively studied over the past 150 years, the evolutionary relationships of the group are poorly understood, impeding our understanding of early vertebrate evolution. A solid phylogenetic framework for the Heterostraci would allow identification of pleisiomorphic states and characters, which could in turn elucidate ancestral conditions on the gnathostome stem. Furthermore, a phylogeny would enable the incorporation of ghost ranges to studies of diversity through time (Sansom et al. 2015) and a test for scenarios of palaeobiogeography (Sansom 2009;Blieck 2011;Zigait_ e & Blieck 2013). Historically, there are five major clades within the Heterostraci; Cyathaspididae, Amphiaspididae, Traquairaspididae, Psammosteidae and Pteraspidiformes, along with many incertae sedis genera of uncertain affinity (Halstead 1973;Janvier 1996). Only a handful of heterostracan phylogenies have been constructed, most of which concentrate on the 'higher heterostracans' À that is, Cyathaspididae and Pteraspidiformes. Containing over 110 species and 47 genera (Dineley & Loeffler 1976;Elliott 1983Elliott , 1984Blieck 1984;Ilyes & Elliott 1994;Elliott & Ilyes 1996;Elliott et al. 2000;Pern egre 2002Pern egre , 2006Novitskaya 2004;Pern egre & Goujet 2007; Pern egre & Elliott 2008;Voichyshyn 2011), Pteraspidiformes are the largest and most-studied group of heterostracans. They are relatively morphologically diverse, with a wide palaeogeographical range (Laurussia, Baltica, Kara and Avalonia) and a stratigraphical range from the Late Silurian to the Middle Devonian (Pr ıdoli-Givetian) (Blieck & Janvier 1999;Elliott et al. 2000;Novitskaya 2007;Sansom 2009;Sansom et al. 2015). The group is characterized as possessing a separate dorsal spine, a dorsal shield composed of several independent plates (dorsal, pineal, rostral and paired orbital, brachial and cornual plates), ornamentation of the dorsal shield in serrated, concentric rings, and supraorbital canals meeting behind the pineal region (Blieck 1984;Janvier 1996;Pern egre & Elliott 2008;Blieck et al. 1991). Within the group is a large variety of forms, some resembling their purported sister group, the cyathaspids, and others taking more unusual forms such as Doryaspis with its ornamented cornual plates and long blade-like pseudorostrum (Pern egre 2002). Anchipteraspididae are the earliest members of Pteraspidiformes, appearing first in the Pr ıdoli of the Canadian Arctic (Elliott 1984). The latest occurring member of the group is Helaspis, found in the Yahatinda Formation of Givetian age in southern Alberta and British Columbia (Elliott et al. 2000). Despite being the most heavily studied and numerous heterostracan group, there is yet to be a phylogeny containing all described genera: only subsets of the taxa have been subjected to phylogenetic analyses (Ilyes & Elliott 1994 & Elliott 2008).
The close affinities of Pteraspidiformes and Psammosteidae have long been recognized, with many phylogenetic reconstructions considering them as sister groups (Halstead 1973;Janvier 1996). Similarities occur in the arrangement of head shield plates, such as separate dorsal, orbital, pineal, rostral, branchial and cornual plates, which are not seen in other heterostracan clades. Also, juvenile psammosteids have been found to possess no fields of tesserae between their plates, thus superficially resembling Pteraspidiformes (Blieck 1984;Janvier 1996). Previous studies have interpreted Psammosteidae as basal or sister taxon to Pteraspidiformes (Blieck 1984;Pern egre & Elliott 2008), while others have recovered them nested well with Pteraspidiformes (Janvier 1996;Pern egre 2002). The uncertainty of their placement means scenarios for heterostracan evolution are somewhat confused. For example, Tarlo (1964) envisaged a transition from tessellated forms such as Tesseraspis and Kallostrakon to Psammosteidae, with their 'fields of tesserae' a vestige of this heritage. If Psammosteidae are found to belong within Pteraspidiformes, then this scenario cannot be corroborated.
Cladistic analysis of Pteraspidiformes is made problematical by similarities in the combination and arrangements of their dermal plates. Many taxonomic distinctions are based on the dimensions of these plates, which are often continuously variable, rather than discrete, states. Anatomical features are often described as short, narrow, elongate, etc., descriptors that without quantification or point of reference are not replicable. One such example of ambiguity is the ratio of dorsal plate width to length, for which 0.7 has been arbitrarily used to distinguish between a broad or narrow dorsal plate (Pern egre & Goujet 2007). Continuous morphological variation such as this is hard to translate into meaningful and discrete cladistic characters and can often lead to uninformative and misleading character coding (Brazeau 2011). Adding to this is the large amount of missing data and poor preservation of some specimens which make the reconstruction of evolutionary histories difficult.
Characters used within phylogenetic analyses are generally discrete; for example, a feature is either present or absent, or is blue or green. There has been considerable debate as to the appropriateness of including continuous characters within cladistic analyses (see Thiele 1993; Rae 1998 for reviews). However, continuous characters can be phylogenetically informative: continuous variation of homologous features meets the requirements of cladistic characters (Thiele 1993;Rae 1998;Garcia-Cruz & Sosa 2006;Goloboff et al. 2006;Sansom 2008;Romano & Nicosia 2015). There are different methodologies for treating continuous data (Thiele 1993;Garcia-Cruz & Sosa 2006). Gap coding involves identifying gaps within a data set that represent boundaries between different character states. One example used by Sansom (2008) was to discretize continuous ratio data in his thyestiid osteostracan phylogeny (another group of jawless vertebrates), identifying gaps in ordered data. These gaps were then translated into changes in character states. However, discretizing continuous data in this manner is problematic as the methods are highly sensitive to missing data and may shoehorn data into evolutionarily uninformative character states. Another method, advocated by Goloboff et al. (2006), is to include the raw continuous data within the phylogenetic analysis (this can be implemented in phylogenetic programs such as TNT; Goloboff et al. 2008).
Here, we undertake a comprehensive analysis of Pteraspidiformes phylogeny that includes all formally described genera. Inclusion of Psammosteidae within Pteraspidiformes is also investigated, with the hope of elucidating the position of this clade. Different coding methods are explored, including analyses with equally weighted characters, implied weighting and inclusion of quantitative data. Quantitative data are further subjected to different treatments including gap coding (transforming continuous data into discrete character states) and using the raw continuous values for taxa. The decision to include quantitative data was made in order to address directly continuous shape variation in Pteraspidiformes dermal plates, which is often used for taxonomic discrimination. With our novel phylogeny, the stratigraphical and palaeobiogeographical occurrences of the Pteraspidiformes are reviewed.

Previous phylogenetic analyses
One of the most comprehensive and thorough studies of pteraspidiform phylogeny and evolutionary history is the monograph of Blieck (1984) in which he systematically described all currently known Pteraspidiformes and considered their stratigraphical and palaeobiogeographical ranges. Two major groupings of Pteraspidiformes were recognized based on characters such as the position and coverage of sensory canals, shape and contact of the pineal plate, rostral size and shape and margins of the plates: Pteraspididae and Protopteraspididae. Blieck (1984) also suggested that Psammosteidae is the sister lineage to Pteraspidiformes (Fig. 1A). Blieck et al. (1991) revised the relationships of the major groupings and placed Pteraspidiformes (pteraspids and anchipteraspids) with psammosteids in a group labelled 'APP'. APP was defined as having two separate dorsal and ventral plates, pineal and rostral plates and paired orbital, branchial and cornual plates, along with a complete pineal canal, which loops between the two branches of the supra-orbital canal, and three pairs of transverse commissures. However, the pineal canal of psammosteids is unknown, perhaps due to preservation, disarticulated nature of psammosteid fossil material or, as Blieck et al. (1991) suggested, migration of the canal from on the plates to the margins between the pineal and rostral plates. The Pteraspidiformes phylogeny of Janvier (1996) was an amalgamation of evolutionary theories for the group at the time (Fig. 1B). Contrasting with Blieck's (1984) theory of Psammosteidae-Pteraspidiformes sister group relationships, Janvier (1996) proposed that Psammosteidae is a derived clade within Pteraspidiformes, with protaspids as a sister group based on the shared possession and arrangement of plates and general similarities to protaspid Pteraspidiformes.
The first cladistic study to utilize global parsimony in a computerized analysis was undertaken by Ilyes & Elliott (1994). They focused on forms from the western USA ( Fig. 1C), citing concerns of convergent evolution amongst Pteraspidiformes relating to ecology rather than phylogeny. A second phylogenetic analysis was undertaken by Pern egre (2002) (Fig. 1D) following a redescription of Doryaspis. Drepanaspis (a psammosteid heterostracan) was included within this analysis and placed as a sister taxon to Doryaspis nested well within the Pteraspidiformes. The analysis included a small sample of Pteraspidiformes but supported Janvier's (1996) view that the Psammosteidae were derived Pteraspidiformes. Pern egre & Goujet's (2007) analysis following re-evaluation of Gigantaspis included many more characters exploring the variation and differences of Pteraspidiformes morphology (Fig. 1E). In their analysis, Gigantaspis was recovered as the sister group to Europrotaspis and protaspid taxa. The most recent study (Pern egre & Elliott 2008) focused on a subset of Pteraspidiformes taxa that are morphologically well known (i.e. are known from many specimens, are articulated, or a reliable reconstruction can be made). They recognized five major clades within the Pteraspidiformes: Anchipteraspididae, a paraphyletic Protopteraspididae, Pteraspididae, Gigantaspididae and Protaspididae (Fig. 1F). Drepanaspis was also included in a second analysis (Pern egre & Elliott 2008, fig. 5) and recovered fairly close to the Pteraspidiformes root, which contrasts with the more derived position recovered by Pern egre (2002). However, when we used the methods and published data matrix of Pern egre & Elliott (2008), Drepanaspis was placed outside the Pteraspidiformes clade, as its sister group. Despite these advances in the systematics of Pteraspidiformes, there has yet to be a phylogeny including all taxonomically described genera, and the position of the Psammosteidae is still ambiguous.

Ingroup taxa
All 47 genera described as belonging to the Pteraspidiformes are included in this analysis (see Supplemental  Table 1). This includes all genera outlined in Blieck 's (1984) extensive Pteraspidiformes review, along with newly described, or redescribed, taxa (Elliott 1983(Elliott , 1984Ilyes & Elliott 1994;Elliott & Ilyes 1996;Pern egre 2002;Pern egre & Goujet 2007;Voichyshyn 2011). Where possible, the type species of each genus has been used, unless that species is inadequately preserved, the material is missing, or another species from within the genus is substantially better known (see Supplemental Table 1). Character coding is based on the selected species within a genus rather than generalizations for a genus or ambiguity due to intraspecific variation. Two representatives of both Protopteraspis and Zascinaspis have been included within the study. Protopteraspis contains many species, some of which exhibit notable variation (e.g. presence of internal organ impressions). The two species of Zascinaspis (Z. heintzi and Z. carmani) are also included within the analysis to test the monophyly of the genus. Zascinaspis carmani was first described extensively by Denison (1960) and originally placed within Pteraspis but considered an intermediate between Pteraspis and Protaspis. A re-evaluation by Blieck (1984)  similarities between the overall shape, pineal-orbital belt and oral region of Z. carmani and Z. heintzi, which he considered enough to include these species within the same genus, Zascinaspis. On examination, similarities such as the scale-like dorsal spine, ornamentation pattern and proportions of the dorsal plate suggest that Z. carmani perhaps does not belong within the genus Zascinaspis; therefore, Z.carmani has been included within the analysis to elucidate its classification. Drepanaspis gemuendenensis and Psammosteus megalopteryx have been chosen as two representatives of the Psammosteidae. These taxa are known from articulated or bountiful material, enabling a clear comparison of the morphology of these forms and the more traditional Pteraspidiformes.

Outgroup taxa
The Pteraspidiformes are often grouped with Cyathaspididae as the 'higher heterostracans' (Janvier 1996), the general consensus being that cyathaspids are the sister group of pteraspids (Elliott 1984;Janvier 1996); they are therefore the logical choice as the outgroup for this analysis. Following Lundgren & Blom (2013), we use Athenaegis chattertoni from the Wenlock of the Mackenzie Mountains, Northwest Territories, Canada, and following Pern egre & Elliott (2008) we use Nahanniaspis mackensiei from the Lochkovian of the Mackenzie Mountains, Northwest Territories, Canada and Anglaspis maccoulloughi from the Pr ıdoli and Lochkovian of the Welsh Borders, UK and the Artois-Ardenne regions of France and Belgium. Athenaegis is the oldest heterostracan known from articulated material, and was placed by Janvier (1996) as sister taxon to his 'higher heterostracan' clade.
Athenaegis was originally placed within Cyathaspididae by Soehn & Wilson (1990), but is believed by Janvier (1996) to possess more plesiomorphic characters. Anglaspis and Nahannisapis are both known from relatively complete material and share many characteristics and states with some Pteraspidiformes, and are thus very useful in polarizing character states within the analysis.

Characters and coding methods
Characters used in the phylogenetic analysis were obtained from direct observations of specimens, and species descriptions from published literature, and were adapted from the character lists from previous phylogenies (Ilyes & Elliott 1994 Goloboff et al. 2006). Twenty-two continuous characters were used in this analysis in the form of ratio data to standardize for size variation (definitions are given in the character list; characters 66À87). Palaeontological data sets are notoriously difficult for morphological phylogenetic analysis with poor preservation and missing data.
Here, different strategies are used to explore the Pteraspidiformes data set including equally weighted characters versus implied weighting (Goloboff 2014), and gap coding to discretize quantitative characters versus inclusion of raw continuous characters (Rae 1998;Wiens 2001;Goloboff et al. 2006;Sansom 2008). Implied weighting is a tool for greater resolution in data sets that contain a large amount of homoplasy and is commonly used in morphological phylogenetics (Dias-da-Silva & Marsicano 2011; Legg et al. 2013;Xu & Pol 2014). This method reweights characters based on their homoplasy index (Goloboff 2014), such that highly homoplastic characters are given less weight. However, it has recently received criticism from studies of simulated data, with concerns raised about the increasing amount of error associated with reweighted character analyses compared With respect to continuous characters, some morphological measurements are used in multiple characters (often as a reference measurement for a ratio) and this can lead to logical dependence or correlation between characters (Goloboff et al. 2006). Therefore, a Pearson's correlation coefficient test was used to identify correlated data (RStudio Team 2015). Significantly correlated characters were removed sequentially, starting with the character with the highest number of correlations. After a character was removed the number of correlations per character was then recalculated, and the highest again removed, and so on. When only single pairs of correlations remained, the character with the least amount of data was then removed, leaving the character with the least amount of missing data. This was repeated until all the remaining characters were uncorrelated and operationally independent.
An alternative strategy to address shape and ratio characters is to discretize them into ordinal character states. This can be done arbitrarily or more formally, for example using gap coding. Traditional gap-coding methods used on neontological data sets (see Garcia-Cruz & Sosa 2006) rely on complete data sets with many replicates of the continuous traits in order to obtain a consensus of the range of variability within a trait, which is not applicable to most palaeontological data sets. The gap coding method used here uses the sequential difference between ordered data points and identifies gaps when the difference is greater than two standard deviations of the gap difference data. This sometimes identified uninformative characters, which were not included within the analysis.
Dermal plates of psammosteids and Pteraspidiformes are treated as homologous following general consensus (Obruchev 1967;Blieck 1984;Blieck et al. 1991). It is not clear, however, whether these plates are functionally analogous, or if the quantitative measurements are evolutionary informative (due to the fields of tesserae between the plates of psammosteids). To account for these two possibilities, analyses were conducted with the continuous characters treated as either homologous or inapplicable for psammosteids (see Supplemental Material for analysis of homologous psammosteid data).

Data analysis
Parsimony tree searches were conducted in TNT (Goloboff et al. 2008) using a heuristic search method (multiple tree bisection rearrangement) with 1000 replications and space for holding 1000 trees. TNT was the favoured program as it can process discrete and continuous data both separately and together (Goloboff et al. 2006). Implied weighting is often used to mitigate the effect of homoplasy (Goloboff 2014). To investigate this possibility we ran searches with and without implied weighting. There is no general consensus as to which concavity value should be used in an analysis; therefore, the default concavity value (k D 3) was implemented (other values for k are presented in the Supplemental Material). Bayesian analysis was performed in MrBayes 3.2 (Ronquist & Huelsenbeck 2003). Priors were kept to their default for standard (morphological) data except the coding function, which was set to 'informative' as the data set used only contains parsimony informative characters. The analysis was run for 2 £ 10 7 generations with samples taken every 1000 generations; the burn-in was set to 0.25.

Results
Searches including all Pteraspidiformes genera recover very little signal, principally due to two taxa (Mylopteraspis and Palanasaspis) that are very incomplete (95 and 96% missing, respectively), acting as wild cards. These two taxa were excluded from all subsequent analyses.
Heuristic searches using discrete characters only resulted in 275 most parsimonious trees (MPTs; tree length 276), giving a partially resolved strict consensus tree (Fig. 4A). Psammosteidae are placed well within the Pteraspidiformes, in a clade with Woodfjordaspis and Doryaspis. The anchipteraspidid taxa Rachiaspis, Anchipteraspis and Ulutitaspis are recovered in a clade sister to all other Pteraspidiformes (including Psammosteidae). A large number of taxa are in a polytomy within the consensus tree with only one other clade apparent containing Larnovaspis and Unarkaspis. Of the suspect genera for which two species were included in the analysis (Zascinaspis and Protopteraspis), neither was found to be monophyletic. Bayesian analysis of the same data set (see Supplemental Fig. 1A for majority rule tree) has better resolution towards the root of the tree, with the Protopteraspis taxa and Stegobranchiaspis placed as a paraphyletic successive sister assemblage leading towards a large polytomous clade containing the remaining taxa. The Psammosteidae are recovered within Pteraspidiformes but within the unresolved polytomy.
Application of implied weighting (k D 3) results in four MPTs with a score of 23.11 (Fig. 4B is the strict concensus). Relationships at the base of the tree are much better resolved than those of the equally weighted analysis. The two Protopteraspis taxa along with Stegobranchiaspis, Loricopteraspis, a clade containing Gigantaspis and Zascinaspis heintzi, Escharaspis and a clade containing Candapteraspis and Miltaspis all form a paraphyletic assemblage leading to a clade containing the majority of the Pteraspidiformes. Seen within this tree is a dichotomous split between the remaining taxa, with the majority of taxa from the western USA (Psephaspis, Cosmaspis, Tuberculaspis, Cyrtaspidichthys, Lampraspis, Zascinaspis carmani, Oreaspis, Lamiaspis, Pirumaspis, Panamintaspis and Blieckaspis) in one clade along with Helaspis and some Podolian forms such as Dnestraspis, Alaeckaspis, Mylopteraspidella, Djurinaspis and Semipodolaspis. The rest of the taxa are contained within the other large clade, which includes the psammosteid taxa and classic Welsh-Border-Artois-Ardenne taxa such as Pteraspis, Errivaspis and Rhinopteraspis. Consistent with the equally weighted analysis is the position of taxa in the clade containing the psammosteids with the addition of Xylaspis, Europrotaspis and Eucyclaspis.
The results for the two different methods for including quantitative characters are shown in Figures 5 and 6. Inclusion of the raw continuous characters recovers a single MPT (Fig. 5A). However, this result is dubious as it has a topology very different from all other analyses. In this tree, taxa that are usually positioned towards the root of the tree (as seen in the Bayesian analysis, Supplemental Figure 1 and the reweighted analyses (Figs 4B, 5B, 6B)), such as Protopteraspis and Stegobracnhiaspis, are replaced by western USA taxa such as Cosmaspis, Lampraspis, Oreaspis and Eucyclaspis, which are usually found in a fairly derived position. Application of implied weighting (k D 3) also results in a single MPT (tree length 26.53) (Fig. 5B) that compares more favourably with the tree resulting from the implied weighted analysis of discrete-only characters. The majority of the unresolved taxa from the discrete-only analysis with implied weighting (Fig. 4B) are resolved following addition of continuous characters. The composition of the clade containing the psammosteid taxa is unchanged (compared to Fig. 4B) but with much better resolution (Fig. 5B); Doryaspis is now placed as the sister group to the psammosteid taxa, with Woodfjordaspis and Xylaspis forming a clade. Of differences relating to other taxa, Blieckaspis is now placed in a clade with Panamintaspis rather than with Pirumaspis, Djurinaspis, Lamiaspis, Mylopteraspidella, Oreaspis and Semipodolaspis.
Discretizing the continuous characters (Fig. 6A) results in a loss of resolution, with many taxa belonging to an unresolved polytomy in the equally weighted analysis, as is also true for the Bayesian analysis of discretized characters (Supplemental Fig. 1B). Application of implied weighting (k D 3) results in much better resolution, with one most parsimonious solution (tree length of 27.86) (Fig. 6B). However, the resulting topology is quite different from that seen in the previous two analyses of weighted characters (Figs 4B, 5B). Taxa towards the root of the tree are unchanged in their positions when compared to the other weighted analyses. However, Miltaspis and Canadapteraspis have migrated towards the root, and occupy a position between Protopteraspis and Loricopteraspis which in previous analyses was occupied by Stegobranchiaspis. Similarly to previous analyses, the psammosteid taxa occur in a clade with Doryaspis, Xylaspis and Woodfjordaspis. However, also occurring in this clade are Pteraspis and Escharaspis, both of which occupied quite different positions in the discrete and discretewith-continuous analyses. The majority of the derived pteraspidiform taxa occur within two clades. One contains taxa mainly from the western USA and Europrotaspis, whereas the other contains all other taxa. The major difference between the discretized analyses (Fig. 6B) and the discrete-only and discrete-with-continuous analyses (Figs 4B, 5B) is the movement of Podolian taxa along with Oreaspis, Pirumaspis and Lamiaspis from the clade containing the majority western USA taxa into the other major clade with the psammosteid clade, Rhinopteraspis clade and other Podolian, Welsh Borders-Artois-Ardenne forms.

Discussion
A major issue of contention with respect to Pteraspidiformes phylogeny is the inclusion of the Psammosteidae. In all instances, Drepanaspis and Psammosteus fall well within the Pteraspidiformes clade; thus, our results support the hypotheses of Blieck et al. (1991), Janvier (1996, Pern egre (2002) and Pern egre & Elliott (2008) that Psammosteidae do indeed belong within the Pteraspidiformes. The psammosteids are consistently placed in a clade with the Spitsbergen genera Doryaspis, Woodfjordaspis and Xylaspis. In the discrete-only and discrete-with-continuous analyses (Figs 4B, 5B), the psammosteids are also united with Eucyclaspis and Europrotaspis. However, in the discretized analysis these two genera move into the clade containing the western USA taxa and are replaced by Escharaspis and Pteraspis. Psephaspis was once believed to belong within the Psammosteidae (Orvig 1961;Tarlo 1964), but was recognized by Denison (1968) and Blieck (1984) as belonging to the Pteraspidiformes based on morphological features such as the possession of a free dorsal spine. Psephaspis is recovered here as belonging to a clade containing psammosteids but only in the equally weighted discrete-withcontinuous analyses; in the remaining analyses it is placed with the western USA protaspids. The position of the psammosteids presented here is most consistent with the tree presented by Pern egre & Elliott (2008). However, it must be noted that under different weighting strategies (i.e. k D 1), searches using discrete or discrete plus continuous characters, psammosteids move into a highly derived position within the tree, within the traditional Protaspididae clade, which agrees with the interpretations of Janvier (1996) and Pern egre (2002). This could indicate that the characters supporting the position of Psammosteidae are highly homoplastic. Nevertheless, a pteraspidiform affinity for Psammosteidae is strongly supported in all analyses. The Anchipteraspididae (Anchipteraspis, Rachiaspis and Ulutitaspis) taxa are consistently positioned at the base of the tree as sister to all other Pteraspidiformes; this supports the hypothesis of Elliott (1984) that they are transitional forms between Cyathaspididae and Pteraspidiformes. However, this could potentially change with the inclusion of more taxa from other heterostracan clades (i.e. Traquairaspididae) which possess primitive 'oak leaf' ornamentation but are similar to the Pteraspidiformes in that their headshield is composed of separate dermal plates (Dineley & Loeffler 1976;Tarrant 1991).
Recovered towards the root in the majority of analyses are paraphyletic Protopteraspis, Stegobranchiaspis, Loricopteraspis, Zascinaspis heintzi, Gigantaspis, Escharaspis, Canadapteraspis and Miltaspis. Exceptions occur in the equally weighted analysis including continuous characters in which these taxa are placed in a clade in a more derived position, and in the implied weighted discretized analysis where the latter three taxa are placed with the western USA taxa. These taxa are in a position analogous to the Protopteraspididae of Blieck (1984) and Pern egre & Elliott (2008). However, Blieck's (1984) Protopteraspididae (Fig. 1A) also contained Unarkaspis, which is recovered in a more derived position in all three parsimony analyses (Figs 4B, 5B, 6B). Also, in Pern egre & Elliott's (2008) phylogeny, Doryaspis, Unarkaspis and Panamintaspis are within the Protopteraspididae, which is not recovered here. Gigantaspis in the majority of other analyses (Fig. 1A, E, F), together with Europrotaspis, is recovered as the sister clade to Protaspididae (Pern egre & Elliott 2008), but not in any of our results. Gigantaspis, like the majority of the taxa of Protaspididae, is much larger and wider than traditional taxa of Pteraspidiformes, and Pern egre & Goujet (2007) suggested that Gigantaspis is perhaps the ancestral morphotype of the Protaspididae.
Zascinaspis carmani in all instances is placed with Lampraspis in a fairly derived position (with the protaspid taxa), which is also supported in analyses using different implied weightings (see Supplemental Figs 1 and 2). The results from these analyses suggest that Z. carmani needs taxonomic re-evaluation.
Phylogenetic analyses of equally weighted characters are poorly resolved. Application of implied weighting (default k D 3) improves the resolution of all the different treatments (Figs 4B, 5B, 6B). Exploring different levels of implied weighting produced different tree topologies for each analysis À that is, discrete or discrete and continuous characters (see Supplemental Figs 2 and 3). Trees arising from analyses with the same implied weighting (e.g. k D 2) were more similar to each other (compare Supplemental Figs 2A and 3A) than trees arising from the same character treatments with different weightings (compare Supplemental Fig. 2A and B). This indicates that the data sets have a fairly high level of homoplasy. Most notable is the movement of Drepanaspis and Psammosteus (Psammosteidae) away from the root under higher or more severe weighting (k < 3).
Inclusion of the quantitative characters within the analyses consistently improved the resolution. Inclusion of raw continuous characters (implied weighted) resulted in one MPT (Fig. 5B). Poor resolution in some areas of the tree, such as the taxa preceding the protaspid clade, in the discrete-only analysis, is most likely due the homogenous nature of the Pteraspidiformes bauplan. This was also found to be the case within Sansom's (2009) phylogeny of Osteostraci, in which he stated "the group is diverse in terms of number of taxa yet very few anatomical features can be used to distinguish between them in any characterbased analysis" (Sansom 2009, p. 112). This appears to be a common problem within ostracoderms and perhaps indicates the need to explore other methods for distinguishing their evolutionary relationships, such as geometric morphometrics. The resolution of the analysis containing raw continuous characters suggests that this information is phylogenetically informative and suitable for analysis of ostracoderm relationships. Interestingly, the results when including continuous characters (k D 3, Fig. 5B) compare closely with the discrete characters only analysis but with improved resolution (k D 3, Fig. 4B). Uncertainty in the placing of Oreaspis, Semipodolaspis and Mylopteraspidella, along with Larnovaspis, Podolaspis, Woodfjordaspis, Xylaspis and Doryaspis, is resolved with the inclusion of the continuous data. It appears that the continuous characters improve the resolution that describes morphological variations used to distinguish different taxonomic genera.
The alternative method of including quantitative characters (i.e. discretizing) was found to perform relatively poorly in the equally weighted analysis (Fig. 6A). Application of implied weighting (Fig. 6B) greatly improved the resolution with one most parsimonious solution. Similar results are seen in the analysis when measurements of Psammosteidae are coded as homologous (Supplemental Fig. 5). Similar to the results of the other analyses (discrete and discrete with continuous characters) is the recovery of a clade containing Drepanaspis, Psammosteus, Xylaspis, Woodfjordaspis and Doryaspis. The position and relationships of taxa recovered here (i.e. a major split of Pteraspidiformes taxa into two clades) is most similar to that of Blieck (1984), Pern egre & Goujet (2007) and Pern egre & Elliott (2008) (Fig. 1A, E, F). A caveat of this method is that it is highly sensitive to missing data, which is prevalent in palaeontological data sets. One of the biggest problems is that it is not always clear if a gap is a genuine evolutionary phenomenon (phylogenetically informative), or the result of missing data or poor taxon sampling. One of the few areas of accordance between this method and the others presented here is the relationships of taxa towards the Pteraspidiformes root. Consideration of Drepanaspis and Psammosteus continuous measurements in the gap-coding process appears to produce a less-resolved tree (see Supplemental Fig. 5A). Incorporation of their measurements changed the coding states in some of the discretized characters. Inclusion of just two of the 18 described genera of Psammosteidae is not very informative and, as stated previously, the configuration of the psammosteid headshield with its fields of tesserae or platelets means that the continuous character measurements are probably not analogous.
Inclusion of the psammosteid measurements in the discrete and continuous analysis (equally weighted, see Supplemental Fig. 4A) resulted in three MPTs; however, the tree topology and resulting relationships differed from those of all other analyses. This is evident by the psammosteid taxa along with Doryaspis, Woodfjordaspis and Xylaspis being placed in a more derived position within a clade containing western USA taxa, which is only seen in implied weighted analyses with a concavity constant of k D 1. Application of implied weighting (k D 3) (Supplemental Fig. 4B) resulted in a topology similar to that where the psammosteid measurements were coded as missing (Fig. 5B). This indicates that the psammosteid measurements (i.e. plate dimensions) may not be homologous to those of other Pteraspidiformes.

Stratigraphical and biogeographical congruence
The trees shown in Figure 5B (implied weighted analysis including continuous characters) and Figure 6B (implied weighted including discretized characters) were calibrated against the stratigraphical and palaeobiogeographical ranges of taxa to explore scenarios of Pteraspidiformes evolution and assess congruence (Fig. 7). Based on the relationships recovered here, along with stratigraphical and biogeographic occurrences, it is apparent that Pteraspidiformes did not show strong endemism, unlike other jawless vertebrate clades (i.e. Osteostraci and Galeaspida) (Sansom 2009). It appears that Pteraspidiformes originated in Laurentia (Canadian Arctic, Mackenzie Mountains, Spitsbergen and western USA localities) and had dispersed to the other palaeocontinents (Avalonia (Welsh Border-Artois-Ardenne), Baltica (Podolia) and Kara (Severnaya-Zemlya)) by the Early Devonian, a very similar dispersal pattern to that seen in osteostracans (Sansom 2009). However, these similar patterns could be due to inherent biases associated with facies differences and an incomplete fossil record (Sansom et al. 2015).
Based on the stratigraphical occurrences and palaeobiogeography of the genera, a more reliable phylogeny À between the discrete-and-continuous phylogeny and the phylogeny resulting from the discretized data À cannot be discerned. There appears to be just as much conflict (based on ghost ranges) and congruence in the two phylogenies (Fig. 7A, B). Palaeobiogeographically, the result from the discretized data appears to be more parsimonious as the majority of taxa from the western USA are placed in one clade, whereas those from Podolia and the Welsh Borders-Artois-Ardenne are placed in another. The discretized phylogeny also has slightly fewer conflicts in terms of ghost ranges and time bins crossed than that from the raw quantitative data. Many forms are common to localities of the Welsh Borders-Artois-Ardenne basin and Podolia, indicating there was faunal exchange between these two palaeobasins, as suggested by Blieck (1984) and Voichyshyn (2011).
The most root-ward and earliest Pteraspidiformes (Anchipteraspididae) occur in the Pr ıdoli of the Canadian Arctic (Elliott 1984). This first occurrence does not coincide with the first occurrences of other clades of Heterostraci, such as Cyathaspididae, Tesseraspis and Athenaegis, which are found in rocks of considerable older age (Wenlock and Ludlovian) (Denison 1963;Dineley & Loeffler 1976;Blieck et al. 2002;Soehn & Wilson 1990). This suggests either Pteraspidiformes evolved later than other clades, or their early history is not recorded perhaps due to facies preservational biases (Sansom et al. 2015).
The majority of derived taxa in the tree are from western USA localities, most of which are resolved as belonging within the Protaspididae of Pern egre & Elliott (2008). These happen to be some of the youngest members of the Pteraspidiformes, spanning from the Pragian to the Eifelian, indicating some congruence between this phylogeny and the stratigraphical record of the Pteraspidiformes. It has long been suggested that the western USA fauna was endemic to this area (Ilyes & Elliott 1994), but our phylogenies indicate that there was some faunal exchange between different palaeobasins (Fig. 7A, B). In the discrete-and-continuous character analysis (Fig. 7A), Eucyclaspis is placed in a clade with Spitzbergen forms, and other taxa are interspersed with taxa from Podolia. In the discretized analysis the western USA taxa appear to be less interspersed with other localities and taxa, apart from Pirumaspis, Lamiaspis and Oreaspis, which are nested within a clade containing predominantly Podolian and Welsh Borders-Artois-Ardenne taxa. Alternatively, as Ilyes & Elliott (1994) D 3). B, discretized analysis with implied weighting (k D 3). Colours relate to palaeobiogeographical provinces. morphology due to similar environmental pressures could have resulted in their positions within the tree. The inclusion of Psammosteidae within the Pteraspidiformes extends the range of the latter by 15 million years; the Psammosteidae are the youngest members of the Heterostraci. Drepanaspis is the first Psammosteidae to appear within the fossil record (Lochkovian; Tarlo 1965) with the last occurring just before the end-Devonian mass extinction in the Frasnian (Tarlo 1965;Moloshnikov 2009;Blieck et al. 2002;Novitskaya 2004).
Some discrepancies between stratigraphical ranges and phylogenetic positions are apparent, which necessitates long ghost ranges. Helaspis from Alberta and British Columbia is the youngest genus of Pteraspidiformes sensu stricto, and the placement of this taxon with Pesphaspis, in a clade with protaspids, results in a fairly long ghost range. These are, however, some of the youngest members of the Pteraspidiformes and this ghost range is most likely due to incompleteness of the rock record (Sansom et al. 2015). The relatively early appearance of the psammosteid Drepanapis in the Lochkovian causes many of the ghost ranges seen within the clade containing the psammosteid taxa in both phylogenies. However, the long ghost range of Psammosteus is most probably due to poor taxon sampling; inclusion of more representatives of psammosteids in the phylogeny would break up the long branch between Drepanaspis and Psammosteus. In the discrete and discretized phylogeny (Fig. 7B), Europrotaspis potentially causes ghost ranges in the wide-shielded protaspid clade.

Systematic palaeontology
Following the phylogenetic analyses, a new classification is proposed. It must be noted that there is ambiguity in the classification of some Pteraspidiformes, which reflects instability in our phylogenetic results.
Revised diagnosis. A Protaspididae with wide shield and branchial opening in advance of the posterior end of dorsal plate; ornament of undulating tuberculated ridges; preoral field which lacks a pre-oral surface; rostral plate with rounded orbital notch; rostral plate well rounded anteriorly and separated from the dorsal plate by an orbitopineal belt; scale-like dorsal spine and scale-like cornual plates.
Remarks. Lampraspis tuberculata has an unornamented pre-oral field rather than the pre-oral surface described by Denison (1970).
Remarks. Phylogenetic analyses placed Zascinaspis heintzi and Z. carmani at different places in the tree (Figs 4À6), indicating they do not belong within the same genus. Blieck (1984) assigned 'P.' carmani to Zascinaspis due to similarities in the proportions of the rostral and dorsal plates and the orbito-pineal belt, as well as shared oral region arrangement. In our analysis, 'Z.' carmani falls within the Protaspididae as sister taxon to Lampraspis (both recovered from western USA localities). Lampraspis carmani also shares with protaspids an ornamentation of undulating ridges and a branchial opening placed towards the posterior end of the dorsal plate, with reduced cornual plates.

Conclusions
The first phylogenetic analysis to include all genera belonging to the Pteraspidiformes is presented here, using different coding strategies to explore morphological shape variations and relationships of the taxa. In all phylogenies resulting from different parsimony coding methods, application of implied weighting gave the most consistent results. The inclusion of quantitative data (i.e. raw continuous characters or discretized characters) led to better resolution relative to searches using discrete qualitative characters only. The implied weighted phylogeny including continuous characters was most comparable to that of the discrete only (implied weighted) phylogeny, whereas discretizing the continuous characters (using objective gap coding) resulted in a different tree topology. Bayesian analyses resulted in similar topologies and resolution to the equally weighted discrete and discretized analyses. In all iterations the Psammosteidae fall within the Pteraspidiformes clade. Congruence between the phylogeny, stratigraphical ranges and palaeobiogeography supports relationships proposed and suggests a Laurentian origin for the Pteraspidiformes. With regards to the different coding approaches, it is not immediately apparent as to whether raw continuous characters are more appropriate than discretized quantitative characters for reconstructing ostracoderm relationships. However, we found discretized character analyses to be highly sensitive to missing data, incongruent with the trees using qualitative data only, but marginally more congruent with stratigraphical data than results using continuous characters. In all cases, homoplasy was found to be widespread, mirroring to some degree other clades of ostracoderms. As such, we recommend implied weighting as an objective tool to mitigate this. Furthermore, given the limited number of variable discrete characters within ostracoderm clades and the high degree of shape variation, quantitative characters need careful consideration and explicit justification.