Functional aspects of the headshield processes in ostracoderms

Abstract

Ostracoderms, Paleozoic jawless stem-gnathostomes, are characterized by distinctive bony shields covering the front of their bodies. These headshields exhibit significant variations in morphology across species, boasting frontal, lateral, and dorsal processes. Ostracoderms represent pivotal intermediaries between modern jawless and jawed vertebrates, so understanding their biology and ecology is crucial for unraveling the selective pressures that shaped the early evolution and diversification of jawed vertebrates, which now dominate vertebrate diversity. This study employs virtual paleontology techniques and phylogenetic comparative methods to explore the hydrodynamic and ecological implications of these processes, focusing on pteraspidomorphs, the most diverse ostracoderm group. The analysis reveals widespread convergence in the arrangement and development of headshield processes. Lateral processes enhance hydrodynamic efficiency and generate lift, while combined lateral and dorsal processes provide stability in rolling, yawing, and pitching. Frontal processes reduce drag in many cases. These findings illuminate the enigmatic roles of ostracoderm headshields, showing how the dimensions and arrangement of their processes are biomechanically linked to a range of functions and ecological roles. Collectively, this highlights the intricate evolutionary pathways of lifestyles and ecologies within stem-gnathostomes, challenging the idea of a unidirectional trend toward more active lifestyles in vertebrate evolution and suggesting diverse ecological roles for ostracoderms.

Non-technical Summary

Ostracoderms, ancient jawless vertebrates from the Paleozoic era, are known for the unique bony shields covering their bodies. These shields vary greatly in shape and structure among different species, featuring frontal, lateral, and dorsal processes. Because ostracoderms represent transitional forms between jawless and jawed vertebrates, understanding their biology is crucial for deciphering the early evolution of vertebrates. This study utilizes virtual paleontology and phylogenetic methods to investigate the hydrodynamic and ecological implications of these shield processes, particularly focusing on pteraspidomorphs, a diverse ostracoderm group. The analysis uncovers widespread convergence in the arrangement and function of these processes. Lateral processes enhance hydrodynamic efficiency and lift, while combined lateral and dorsal processes aid in stability during movement. Frontal processes reduce drag in many cases. These findings shed light on the complex roles of ostracoderm shields, revealing how their dimensions and arrangement are linked to various functions and ecological niches. This challenges the notion of a linear trend toward more active lifestyles in vertebrate evolution and suggests a diverse array of ecological roles for ostracoderms. Overall, this study provides insight into the evolutionary pathways and ecological diversity of early vertebrates.

Introduction

The origin and evolution of the vertebrate body plan have, for centuries, been central topics of debate in evolutionary biology (Janvier 1996; Ahlberg 2001; Donoghue and Keating 2014 and references therein). The early evolution of vertebrates involved significant changes in development and repeated episodes of whole-genome duplication (Donoghue et al. 2008; Heimberg et al. 2010; Donoghue and Keating 2014; Marlétaz et al. 2018). The “new head” hypothesis posits that these transformations were connected to a prolonged ecological shift, transitioning from suspension feeding in a lancelet-like proto-vertebrate to progressively more active locomotion and predation-oriented lifestyles (Gans and Northcutt 1983; Northcutt and Gans 1983; Gans 1989, 1993; Northcutt 1996, 2005). This hypothesis draws upon insights from developmental biology, neurobiology, functional morphology, and systematics, complemented by fossil data, to construct a narrative that elucidates the sequence of events and the selective pressures that underpinned the emergence of vertebrates. Similar conceptual frameworks have also been advanced by certain researchers to elucidate the origin of gnathostomes (Mallatt 1984a,b, 1985, 1996). These proposals contend that the majority of adaptations driving the evolution of vertebrate jaws were geared toward enhanced ventilation, propelling the evolutionary journey of gnathostomes through intensified activity and predation-driven selection pressures. Ostracoderms, representing a series of paraphyletic groups of jawless stem-gnathostome vertebrates, play a crucial role in contrasting these hypotheses. Phylogenetically intermediate, they bridge the gap between living agnathans (lampreys and hagfishes) and gnathostomes (the rest of modern vertebrates; i.e., osteichthyans and chondrichthyans), showcasing the progressive acquisition of gnathostome traits across a continuum of body designs (Donoghue and Keating 2014). Ostracoderms have conventionally been interpreted as benthic mud-grubbers with poor swimming capabilities and low maneuverability (Aleyev and Novitskaya 1983; White and Toombs 1983; Belles-Isles 1987; Mark-Kurik 1992), and the origin of jaws is frequently presented as the key innovation that precipitated the ecological diversification and dominance of the group (Anderson et al. 2011; Friedman and Sallan 2012). Nonetheless, our comprehension of ostracoderm biology remains limited (Janvier 1996), and attempts to draw ecological conclusions about these groups have yielded primarily anecdotal evidence (Northcutt 2005) or, at best, contentious findings (Purnell 1995, 2001a,b, 2002; Ferrón and Botella 2017; Ferrón et al. 2018, 2020, 2021; Miyashita et al. 2021; Ferrón and Donoghue 2022). Consequently, clarifying the morphological adaptations of these extinct groups to their environment is key to the understanding of the mechanisms underlying the origin of jawed vertebrates.

Most ostracoderms were armored Paleozoic agnathans characterized by having a large rigid external shield enclosing the anterior portion of their bodies (Janvier 1996). These structures exhibited a wide morphological diversity, with shapes that varied from oblate to prolate, with conspicuous frontal, lateral, and/or dorsal processes often interpreted as being for stabilization, enhancing lift or reducing drag (Afanassieva 1992; Mark-Kurik 1992; Dec 2019a; Ferrón et al. 2020, 2021), deterring predators (Janvier 1977), housing sense organs (Voichyshyn 2006), substrate anchoring (Janvier 1985), or specialized feeding strategies (Dineley 1994). Notwithstanding this, the absence of modern analogues for most stem-gnathostomes makes it difficult to draw reliable functional inferences about these groups. As a result, our understanding of the role of ostracoderm headshields and, more specifically, headshield processes, remains restricted, despite their substantial potential to offer insights into the ecological scenario of early vertebrate evolution. To address this, we attempt to provide insight into the ecology of stem-gnathostomes and the evolutionary scenario of the “new head” hypothesis by evaluating the hydrodynamic performance and the functional implications of the headshield processes of ostracoderms. Here, we focus on Pteraspidomorphi, which represent the most diverse group of agnathan stem-gnathostomes both in terms of species and the various morphologies of their headshields, encompassing the complete spectrum of shapes displayed by ostracoderms. To achieve this, our approach involved two main steps: initially, we assessed the degree of homoplasy in these structures by employing geometric morphometrics and various quantitative convergence metrics. This evaluation aimed to determine whether headshield processes had independently evolved in distinct groups. Subsequently, we delved into their functional and ecological significance by employing computational fluid dynamics simulations and tools for biomechanical inference. These analyses were conducted on models that included intact headshield processes, as well as on models where these processes were digitally removed.

Materials and Methods

Three-Dimensional Virtual Modeling

Pteraspidomorphs, spanning from the Early Ordovician to the Late Devonian, represent a clade of stem gnathostomes characterized by the possession of leaf-shaped tubercles on the headshield. Usually, large median ventral and dorsal dermal plates constitute most of their head armor, and only the caudal fin is present. Pteraspidomorphs comprise four major clades (Arandaspida, Astraspida, Eriptychiida, and Heterostraci), and the majority of which were marine species, typically confined to nearshore areas such as lagoons. However, a few are regarded as freshwater forms that inhabited river systems or deltas (Janvier 1996). We created a total of 60 digital 3D models that encompass virtually all known pteraspidomorphs that preserve headshields from which 3D information can be gleaned (Supplementary Data S1, Supplementary Table S1). Our sample includes 57 heterostracans (including cyathaspids, amphiaspids, and pteraspidiformes) plus two arandaspids (Sacabambaspis, Arandaspis) and one astraspid (Astraspis). Digital models were built using Rhinoceros 7¹ by leveraging published reconstructions or photographic references of fossil specimens in multiple viewpoints (i.e., frontal, lateral, and dorsal). Research has demonstrated that digital models created through 3D virtual modeling serve as dependable instruments for evaluating function via computational analysis, yielding outcomes highly akin to those obtained through tomographic or surface-based methods (Rahman and Lautenschlager 2016). The majority of the chosen species are known by specimens with remarkably preserved fossil headshields, guaranteeing precise modeling of these elements. The postcranial region was accurately modeled for taxa in which it is known, while for the remaining species, the morphology of their closest phylogenetic relative was utilized. All models were repaired and converted into NURBS surfaces using Geomagic Studio 2012² and scaled to 100 mm total body length, oriented, centered with respect to their volume centroids, and converted into STEP format using Rhinoceros 7.

From the dataset, we handpicked 37 models that fulfilled one or more of the following criteria: the presence of prominent or extensively developed rostral plates, cornual plates, and/or spinal plates (Fig. 1A). These encompass frontal, lateral, and dorsal headshield processes that could potentially impact the hydrodynamics of the organism. To explore this further, we crafted modified versions of these models, each with a single headshield process removed, enabling subsequent computational fluid dynamics (CFD) analyses (Fig. 1 A–D, Supplementary Data S1). This allowed us to examine the individual impact of each process on their hydrodynamics. This subsample includes representatives from Arandaspida, Astraspida, and all major groups of Heterostraci, such as Cyathaspididae, Amphiaspididae, Traquairaspididae, Psammosteidae, and Pteraspidiformes, as well as several “tessellate” taxa.

Figure 1. Preprocessing of 3D models and landmark configurations. A, Visualization of the 3D model showcasing frontal, dorsal, and lateral headshield processes, highlighted in yellow, red, and orange, respectively. The original 3D models underwent preprocessing steps before computational fluid dynamics (CFD) simulations involving the removal of (B) frontal, (C) lateral, and (D) dorsal headshield processes. E, The landmark configurations utilized in the geometric morphometric analyses are depicted. As an illustrative example of a pteraspidomorph, the 3D model of Errivaspis waynensis is presented in lateral, frontal, and dorsal views, arranged from left to right.

Geometric Morphometrics

We used geometric morphometrics to numerically describe the morphological diversity of Pteraspidomorphi headshields. The analyses were conducted at the generic level, selecting one specimen per genus based on the best preservation. Three different landmark configurations were obtained for each genus in frontal, lateral, and dorsal views, respectively, using the projections of each 3D model in the space axes (x, y, and z). Our landmark configurations included a total of 70 landmarks of type II and III that were equally interpolated along the specimen’s outline using TpsUtil v. 1.58 and TpsDig v. 2.26 (Rohlf 2015; Fig. 1E). Generalized Procrustes analysis (GPA) was performed in R (R Core Team 2022) using the package geomorph v. 4.0.8 (Adams et al. 2021) in order to remove the variation in translation, rotation, and size from the original landmark configurations. We enriched our analysis by synergistically integrating morphological information from the different landmark configurations through the “combinland” approach (Profico et al. 2019), implemented in the R package Arothron v. 2.0.5 (Profico et al. 2021). This method integrates shape information from multiple 2D landmark configurations into a unified dataset, preserving the geometric accuracy of the original 3D shapes. Initially, each 2D landmark configuration, derived from distinct anatomical views, undergoes GPA to align the shapes. To ensure comparability across datasets with varying numbers of landmarks, the combinland approach applies a size correction, adjusting for the number of landmarks and dimensions. These corrected datasets are concatenated into a single matrix and analyzed using principal component analysis (PCA). This process effectively combines the shape data from different views, enabling the retrieval of 3D-quality morphological insights from 2D data. We favor this method over conventional 3D approaches due to its remarkable efficiency and speed, as it utilizes 2D images to quickly and accurately capture detailed contours with precision equivalent to that of 3D models, while empirical studies confirm that it effectively preserves morphological variation with minimal information loss (Profico et al. 2019), thus offering a highly reliable alternative to 3D digitization.

Morphological Homoplasy Analyses

To investigate the correlation between morphospace occupation and phylogeny, phylomorphospaces were constructed and visualized using the R packages ape v. 5.8 (Paradis and Schliep 2019), phytools v. 2.3.0 (Revell 2012), and ggplot2 v. 3.5.1 (Wickham 2023). Our tree is based on the phylogenetic hypothesis proposed by Randle et al. (2022), in which the Heterostraci, Amphiaspididae, Traquairaspididae, and Psammosteidae constitute monophyletic groups (the last two as part of Pteraspidiformes, which also includes pteraspids). Cyathaspididae, by contrast, is considered a polyphyletic group. This phylogeny was augmented by incorporating specific psammosteids, including Pycnosteus and Tartuosteus as detailed in Glinskiy (2018) and Guerichosteus from Dec (2019b), along with astraspids and arandaspids following Janvier (1996). Mesquite v. 3.81 (Maddison and Maddison 2023) was employed to construct the phylogenetic tree. This tree was then time-calibrated using the paleotree v. 3.4.7 (Bapst 2012) and ape v. 5.8 (Paradis and Schliep 2019) R packages, based on stratigraphic ranges published in Randle et al. (2022). We employed the bin_timePaleoPhy function from paleotree, utilizing the mbl method for minimum branch length calibration. This method scales all branches to be greater than or equal to vartime (set to 1 in our analysis) and subtracts the time added to later branches from earlier branches to maintain the temporal structure of events. We generated 1000 trees to account for uncertainty, with stochastic binning disabled. The analysis was set to randomly sample dates within observed stratigraphic ranges. Post-calibration, the trees were dichotomized using multi2di function from ape to resolve any remaining polytomies (Supplementary Data S2).

We employed two distinct methodologies to assess the degree of homoplasy within our dataset: the Mantel test and Stayton’s C1 and Ct1 metrics. Initially, we gauged the extent of dissociation between the distance matrices derived from both the morphometric and phylogenetic datasets. To achieve this, we conducted Mantel tests utilizing the R packages ape v. 5.8 (Paradis and Schliep 2019), stats v. 4.4.1 (Bolar 2019), and vegan v. 2.6.6.1 (Oksanen et al. 2022) using 1000 phylogenies that accounted for time calibration and phylogenetic uncertainty. To do this, polytomies in the original tree were randomly resolved, and each tree was calibrated by randomizing the tip age of each species within the chronostratigraphic unit corresponding to the species’ first appearance. The presence of homoplasy is indicated by enhanced decoupling in the distance matrices, leading to correspondingly reduced correlations. Next, we utilized Stayton’s C1 metric (Stayton 2015) to quantify the relationship between the phenotypic distances among pairs of potentially convergent taxa ( $ {D}_{\mathrm{tip}} $ ) and the utmost phenotypic distance observed between any pair of ancestors ( $ {D}_{\mathrm{max}} $ ) as

(1) $$ C1=1-\left(\frac{D_{\mathrm{tip}}}{D_{\mathrm{max}}}\right). $$

The values of this metric span from 0 to 1, representing the extent of convergence between pairs. We also used the recently proposed Ct1 metrics, which incorporate time when assessing convergence (Grossnickle et al. 2024). Unlike traditional C-measures, which compare morphological distances without considering timing, Ct-measures limit comparisons to synchronous time slices coinciding with internal phylogenetic nodes. This approach reduces the likelihood of misidentifying divergence as convergence. Additionally, while traditional C-measures might assign a score of zero to divergent lineages when $ {D}_{\mathrm{max}} $ equals $ {D}_{\mathrm{tip}} $ , Ct-measures address this issue by using $ {D}_{\mathrm{t}.\max } $ , which represents the maximum phenotypic distance observed across all time slices but excludes tip-to-tip distances. This ensures that $ {D}_{\mathrm{t}.\max } $ cannot equal $ {D}_{\mathrm{tip}} $ , resulting in negative scores for divergent lineages and providing a more accurate assessment of divergence.

To evaluate the statistical significance of the calculated C1 and Ct1 values, we executed 100 simulations for each pair using Brownian motion models. This was done to determine whether the observed C1 and Ct1 values were higher than what would be expected by random chance. We used the first 11 principal component scores for this analysis, which accounted for greater than 95% of the total variance in the sample. These analytical simulations were carried out using the convevol v. 1.0 (Stayton 2023) R package.

Biomechanical Analyses

CFD constitutes a computational technique aimed at simulating the flow of fluids (both liquids and gases) and their interactions with solid surfaces. This approach facilitates extensive comparative analyses, enabling the resolution of functional and ecological hypotheses concerning extinct taxa on a broad scale (Rahman 2017). We conducted simulations involving the flow around the 3D models (both with and without headshield processes) using ANSYS-Fluent 2020 R1 Academic.³ The methodology described by Gai et al. (2022) was employed for this purpose.

The computational domain was a rectangular prism, measuring 1820 mm in length and 610 mm in both width and height. The model was placed centrally within this domain, oriented horizontally. We conducted three different simulations with varying pitch and yaw angles: (1) pitch 0°, yaw 0°, (2) pitch 15° with a nose-up rotation, yaw 0°, and (3) pitch 0°, yaw 15° with clockwise rotation in the dorsal view. To accurately capture the wake, an additional refinement volume was included. This volume took the form of a rectangular prism with dimensions: length of 700 mm, height of 200 mm, and width of 120 mm. At the upstream end of the domain, an inlet boundary condition was set, specifying a normal inflow velocity and a turbulence intensity of 0.05. At the opposite end of the domain, a zero-pressure outlet boundary condition was applied. The sides of the rectangular solid were treated with slip symmetry boundary conditions, simulating a zero-shear wall effect. The walls of the model were given a no-slip boundary condition, fixing the fluid velocity at zero (Fig. 2A).

Figure 2. Experimental setup and computational domain for computational fluid dynamics (CFD) simulations. Schematic representations of (A) the computational domain and (B) the mesh, with detailed view highlighting the inflation layers encompassing the 3D model. Boundary conditions are designated as follows: if, inflation layers; in, inlet; ns, no-slip; ou, outlet; rf, refinement volume; ss, slip symmetry. Modified from Gai et al. (2022).

The ANSYS-Meshing software was employed to generate the mesh for the domain (Fig. 2B). Tetrahedral elements were utilized for both the interior and the walls of the rectangular prism. To capture the boundary layer phenomena, an inflation layer was introduced along the no-slip boundaries (i.e., the pteraspidomorph model’s walls). This layer consisted of up to 23 layers of prismatic elements, each with a growth rate of 1.1. The mesh structure was then transformed into a polyhedral arrangement. For details regarding the mesh used in the CFD simulations, please refer to Supplementary Table S2. Independence tests for mesh size, domain size, and refinement volume for this setting have been performed by Gai et al. (2022).

We conducted simulations of 3D, incompressible flow within the domain. A stationary solver with double precision was employed alongside a second-order discretization method to compute the steady-state flow. The Spalart–Allmaras model and a pressure-based segregated algorithm were utilized to solve the Reynolds-averaged Navier–Stokes equations. Convergence of the simulations was assessed based on stable numerical solutions for integrated drag and lift values, root-mean-square residual levels at 10⁻⁴, and a mass flow rate imbalance below 1%. For these evaluations, we adopted a flow velocity of 0.3 m s⁻¹, corresponding to a Reynolds number of approximately 30,000. This velocity was chosen to align with swimming speed records of similarly sized living fishes (Videler 1993).

We extracted drag (D), lift (L), lift to drag ratio (L/D) values, as well as pitch and yaw moments (M _pitch and M _yaw, respectively). D denotes the resistance force in the direction of fluid flow, L refers to the upward force perpendicular to the fluid flow, and L/D is a measure of an object’s efficiency in generating L compared with the D commonly used in aerodynamics, and M _pitch and M _yaw encompass the rotational effects around the z and y axes (i.e., the tendency to rotate or twist around specific axes). Positive M _pitch values indicate a tendency toward dorsal (nose-up) rotation, while negative values indicate ventral (nose-down) rotation. Similarly, positive M _yaw values indicate a tendency toward clockwise rotation, and negative values indicate counterclockwise rotation, both in the dorsal view. M _pitch and M _yaw were computed with respect to the volume centroid, which was taken as the center of mass, assuming a uniform density distribution in our models. Considering the multitude of scenarios and parameter combinations, our attention is directed toward specific parameters within particular scenarios that hold the potential to offer insights into particular hydrodynamic aspects of interest for the aims of the present study, following prior research (Barnard and Philpott 2010). Processes located dorsally can potentially impact yawing stability (we evaluated M _yaw under conditions of pitch at 0° and yaw at 15°); processes positioned frontally might affect drag reduction (we evaluated D under conditions of pitch at 0° and yaw at 0°); processes situated laterally could influence lift generation and lift-to-drag ratio with increasing pitch (we evaluated L and L/D under conditions of pitch at 15° and yaw at 0°), as well as pitching stability (we evaluated M _pitch under conditions of pitch at 15° and yaw at 0°) and yawing stability (we evaluated M _yaw under conditions of pitch at 0° and yaw at 15°). In this context, negative values for ΔM _pitch and ΔM _yaw indicate a stabilizing effect (i.e., rotating the model back toward the flow direction), while positive values indicate a destabilizing effect (i.e., rotating the model further away from the flow direction).

Moreover, to assess rolling stability and maneuverability, we determined the rotational inertia around the rolling axis (I_x), commonly referred to as moment of inertia. This relates to an object’s resistance to changes in its rotational motion. It is a measure of how an object’s mass is distributed with respect to the axis of rotation. Objects with larger rotational inertia are harder to rotate or stop from rotating compared with objects with smaller rotational inertia. Rotational inertia depends not only on the mass of the object but also on how that mass is distributed in relation to the axis of rotation. Objects with mass concentrated farther from the axis of rotation have larger rotational inertia. I_x was calculated under conditions of pitch at 0° and yaw at 0°, implementing the command “VolumeMoments” in Rhinoceros 7, assuming a uniform density distribution in our models.

We followed two different approaches to evaluate and visually depict the differences in these biomechanical parameters between the original models and the models with processes removed (ΔD, ΔL, ΔL/D, ΔM _pitch, and ΔM _yaw, representing the impact of these processes on pteraspidomorph hydrodynamics). Initially, we fit linear models between the parameter differences and the respective process areas. A significant effect is anticipated if these structures indeed influence the biomechanical parameter in question. We utilized R’s built-in linear modeling functions and showcased the outcomes using the package ggplot2 v. 3.5.1 (Wickham 2023). To address potential phylogenetic autocorrelation in the residual error terms of these linear models, we also conducted analyses using phylogenetic generalized least squares (PGLS). This approach was chosen to appropriately account for the phylogenetic relationships among taxa. Specifically, we used the gls function from the package nlme v. 3.1.164 (Pinheiro et al. 2024), incorporating Pagel’s lambda to account for the phylogenetic autocorrelation. These adjustments ensure that the models account for phylogenetic autocorrelation, reducing to ordinary least squares if no autocorrelation is present or approximating Felsenstein’s contrasts if phylogeny fully explains the autocorrelation. We incorporated methods to address the issue of unequal variances in the residuals, enhancing the robustness of our analysis. Alternatively, we constructed heat maps of parameter difference values across the morphospace using the R package akima v. 0.6.3.4 (Akima et al. 2022), effectively highlighting the correlations between morphology and biomechanics.

We utilized the color blind–friendly inferno palette from the viridis v. 0.6.5 package (Garnier et al. 2024) in all figures. Supplementary Data S2 contains all the R scripts and associated files, encompassing the analysis workflow as well as the code to generate the figures.

Results

Geometric Morphometrics and Homoplasy Tests

The first three principal components (PCs) account for 72.1% of the total variance, with PC 1, PC 2, and PC 3 explaining 40.5%, 21.2%, and 10.4%, respectively. In the morphospace (Fig. 3), the lowest PC 1 scores are associated with headshields resembling torpedoes, exhibiting elongation in dorsal view and lateral compression in cross section. These headshields feature well-developed rostral and dorsal processes. Conversely, the highest PC 1 scores correspond to headshields that are shorter and wider, displaying pronounced dorsoventral flattening and well-developed lateral processes. Minimum PC 2 scores are indicative of deltoid-winged headshields with a triangular cross section and extended lateral and dorsal processes. In contrast, the highest PC 2 scores align with headshields that are elongated and ellipsoidal in dorsal view, showcasing an inverted triangular cross section, along with smaller dorsal processes. At the lowest end of PC 3 scores, headshields show a rounded shape in both dorsal view and cross section, with slightly developed dorsal processes. Conversely, the maximum PC 3 scores are linked to headshields resembling triangles in dorsal view. These triangular headshields exhibit a dorsoventrally flattened rhomboidal cross section and frontal and dorsal processes.

Figure 3. Morphological disparity of headshields in Pteraspidomorphi. Outcomes of geometric morphometric analysis and phylomorphospace representation summarizing the morphological diversity within distinct groups of Pteraspidomorphi. The Pteraspidomorphi category encompasses Arandaspida, Astraspida, and Heterostraci taxa. Within Heterostraci, the subdivisions include Cyathaspididae (in yellow), Amphiaspididae (in orange), Pteraspidiformes (consisting of Pteraspididae, Psammosteidae, and Traquairaspididae, in red), and a few taxa with uncertain affiliations within the group (labeled “other Heterostraci”). This taxonomic classification follows Randle et al. (2022). The lower part of the figure highlights shape transformations associated with the minimum (min) and maximum (max) scores of the principal components (PCs).

Significant overlap is evident within the morphospace among different major groups of pteraspidomorphs (Fig. 3). Arandaspida and Astraspida are positioned with average, positive, and negative scores along PC 1, PC 2, and PC 3, respectively. Meanwhile, Amphiaspididae and Cyathaspididae share an analogous positioning to that of the latter but form two distinct clusters with minor overlap. These groups are characterized by rounded headshields with relatively underdeveloped processes. The majority of the morphospace is occupied by Pteraspidiformes. Among them, Pteraspididae stands out as the most widespread group, showcasing a diverse array of forms defined by the general development of cephalic processes. This includes well-formed rostral plates that give rise to elongated and pointed snouts, along with prominent cornual plates and conspicuous dorsal spines. In addition, Psammosteidae display positive, average, and negative scores along PC 1, PC 2, and PC 3, respectively. They feature laterally pointing-downward expanded headshields characterized by the absence of dorsal processes and limited development of the rostral plate. Finally, Traquairaspididae exhibit average scores along PC 1 and PC 2, and negative scores along PC 3. They display elongated rostral plates, cornuals, and even spinal plates, featuring a discreet dorsal crest that contributes to an almost pentagonal cross-sectional morphology.

The phylomorphospace displays a pattern of generalized convergence, with tree branches conspicuously entwined (Fig. 3). The outcomes of the Mantel tests indicate a lack of significant correlation between phenetic and phylogenetic distances. The Mantel statistic r values span from −0.028 to 0.095, while p-values range between 0.126 and 0.638 (Fig. 4). An examination of Stayton metrics further reveals that an important number of the analyzed pairs exhibit C1 and Ct1 values that significantly exceed what would be anticipated under a Brownian motion evolution scenario. Thus, 818, 375, and 44 pairs have C1 values above 0.25, 0.50, and 0.75, respectively, while 286, 169, and 49 pairs have associated p-values below 0.10, 0.05, and 0.01 (Fig. 5A); and 86, 51, and 44 pairs have Ct1 values above 0.25, 0.50, and 0.75, respectively, while 225, 128, and 46 pairs have associated p-values below 0.10, 0.05, and 0.01 (Fig. 5B). A considerable portion of these instances involve taxa showing elongated dorsal processes (e.g., Pteraspis–Podolaspis, Zascinaspis–Rachiaspis, Ulutitaspis–Europrotaspis), frontal processes (e.g., Podolaspis–Lamiaspis), and lateral processes (e.g., Unarkaspis–Canadapteraspis); in all these cases, the pairs exhibit p-values below 0.01 for both C1 and Ct1 metrics. These processes typically occur in tandem and are seldom observed as isolated traits. For instance, specimens exhibiting curved lateral processes pointing backward, along with blunt snouts and vertically oriented dorsal processes (e.g., Canadapteraspis, Larnovaspis, and Unarkaspis), converge toward negative scores of the PC 1 and PC 3. On the other hand, specimens featuring triangular lateral processes that are not particularly prominent laterally, pronounced and pointed snouts, and dorsal processes curving backward (e.g., Pteraspis, Errivaspis, Lamiaspis, Podolaspis, and Panamintaspis) converge toward negative scores of the PC 1 and positive scores of the PC 2 and PC 3. Meanwhile, specimens with a more rounded cephalic shape, where the curved backward lateral structures harmonize with the headshield’s contour, rounded snouts that blend with the overall head structure, and pronounced backward-curving dorsal processes (e.g., Rachiaspis, Zascinaspis, and Eucyclaspis), converge into the negative scores of the PC 1 and PC 3. Conversely, notable convergence (significance indicated by p-value < 0.01 and C1 and Ct1 values > 0.6) also exists among forms with limited or indistinct cephalic processes (e.g., Putonaraspis–Ariaspis, Nahaninaspis–Anglaspis).

Figure 4. Quantifying convergence in headshield processes of Pteraspidomorphi through Mantel test. Mantel statistic r values and their associated p-values were computed considering 1000 trees, accounting for uncertainties in time calibration and phylogenetic relationships.

Figure 5. Quantifying convergence in headshield processes of Pteraspidomorphi through Stayton’s (A) C1 and (B) Ct1 metrics. The Stayton’s C1 and Ct1 metrics, depicted in purple, provide a measure of convergence strength, while the accompanying p-values, shown in yellow-orange, indicate the statistical significance of the observed convergence. Note that instances where Ct1 metrics and p-values are displayed in white result from uninformative analyses, typically due to comparisons between sister taxa or taxa lacking contemporaneous ancestral nodes. The phylogenetic tree is based on Randle et al. (2022), Glinskiy (2018), Dec (2019b), and Janvier (1996).

Biomechanical Analyses

No statistically significant effects have been observed in relation to the frontal and dorsal process areas on ΔD and ΔM _yaw (slopes = 3.20*10⁻⁴ and 6.10*10⁻⁵; p-values = 0.10 and 0.11; respectively) (Fig. 6A,B, respectively). In both instances, a notable dispersion of data points is evident, with various taxa exhibiting both negative and positive ΔD and ΔM _yaw values (R ² = 0.08 and 0.06). Upon examination of the corresponding heat maps, distinctively heterogeneous patterns emerge, indicating the absence of an apparent correlation between morphology and biomechanics within the morphospace (Fig. 6A,B). Significant effects, along with relatively limited point dispersion, have been identified in the remaining instances. Notably, the area of the dorsal process has a positive effect on ΔI_x (slope = 0.13; p-value = 3.40*10⁻¹⁰; R ² = 0.74) (Fig. 6C). On the other side, the area of the lateral process has a positive effect on ΔL, ΔL/D, and ΔI_x (slopes = 1.90*10⁻², 6.40*10⁻⁴, and 0.66; p-values = 5.10*10⁻²¹, 4.60*10⁻¹¹, and 1.20*10⁻¹²; R ² = 0.95, 0.77, and 0.82, respectively), yet conversely has a negative effect on ΔM _pitch and ΔM _yaw (slopes = −8.80*10⁻⁵ and −4.20*10⁻⁵; p-values = 1.70*10⁻⁴ and 4.20*10⁻⁴; R ² = 0.37 and 0.33, respectively) (Fig. 6D–H, respectively). Aligned with these findings, a close examination of the associated heat maps reveals that the highest values for ΔI_x, ΔL, and ΔL/D tend to cluster in areas of the morphospace corresponding to significant development of either the dorsal or lateral processes (i.e., negative scores of the PC 1 and PC 2 and average scores of the PC 3; and positive scores of the PC 1 and PC 3 and negative scores of the PC 2; respectively) (Fig. 6D–F). Conversely, the regions of the morphospace associated with the most pronounced lateral process development (i.e., positive scores of the PC 1 and PC 3 and negative scores of the PC 2) consistently display the lowest values for ΔM _pitch and ΔM _yaw (Fig. 6G,H). These results remain consistent even when accounting for phylogenetic autocorrelation in the models, indicating that the observed relationships between process areas and biomechanical parameters are robust and not merely artifacts of phylogenetic relatedness (Supplementary Fig. S1).

Figure 6. Hydrodynamic impact of headshield processes of Pteraspidomorphi. Differences in various biomechanical parameters between the unaltered original models and those with specific processes removed are depicted through linear regression plots against the corresponding process area (depicted on the left), while complementary heat maps are overlaid onto the morphospaces (illustrated on the right). Impact of the frontal process on (A) the difference in drag; impact of dorsal process on the difference in (B) yawing moment and (C) rotational inertia; and impact of the lateral processes on the difference in (D) lift, (E) lift-to-drag ratio, (F) rotational inertia, (G) yawing moment, and (H) pitching moment. Negative values for ΔM _pitch and ΔM _yaw indicate a stabilizing effect, while positive values indicate a destabilizing effect.

Discussion

Our analyses unveil a widespread convergence in both the morphologies of pteraspidomorph headshields and the arrangement and extent of development exhibited by their processes. This is supported by qualitative phylomorphospace interpretation (Fig. 3) and quantitative metrics (Figs. 4, 5). The Mantel test indicates weak correlation, suggesting substantial convergence within the clade (Fig. 4). Stayton’s metrics analysis strengthens this by showing that numerous pteraspidomorphs evolved to closely resemble each other, defying stochastic evolutionary expectations (Fig. 5).

The exploration of biomechanical properties associated with these structures suggests that, in many cases, their repeated evolution may have been prompted by adaptation to specific hydrodynamic and locomotory conditions (Fig. 6). Our findings demonstrate that in the majority of species with lateral processes, these structures enhance lift and optimize the lift-to-drag ratio at low pitching angles, similar to the effect of aircraft wings (Barnard and Philpott 2010) (with positive effects on ΔL and ΔL/D; Fig. 6D,E). This high lift-to-drag ratio confers a hydrodynamic advantage by decreasing the energy required to maintain a specific position in the water column. Moreover, lateral processes offer stability in pitch and yaw (with negative effects on ΔM _pitch and ΔM _yaw; Fig. 6G,H). Dorsal processes contribute to yaw stability in certain species (when ΔM _yaw is below 0), despite the absence of a general effect (Fig. 6B). The positioning, orientation, and shape of lateral and, in some cases, dorsal processes facilitate stability along the yaw and pitch axes, functioning as directional stabilizers. This is particularly evident in species where these processes are situated posteriorly, likely behind the animal’s center of mass, which would have been shifted even farther anteriorly by the presence of thick bony headshields compared with our models that assume homogeneous density. In this position, their role as stabilizers is reinforced as they generate corrective torques and stabilizing moments that allow steady forward movement and counteract the yaw and pitch induced by the movement of the caudal pedicle, which could compromise hydrodynamic stability (Aleyev 1977; Webb et al. 1996; Fish 2002; Fletcher et al. 2014; Fish and Lauder 2017). This mechanism bears resemblance to an aircraft’s empennage, where the vertical stabilizer controls yaw and the horizontal stabilizers manage pitch, as influenced by the wings (Struett 2012; Ciliberti et al. 2017; Liu 2020). The combined contribution of both processes in stabilization is further evidenced by a positive correlation between the M _yaw of the dorsal process and the M _yaw and M _pitch of the lateral processes (slopes = 0.354 and 8.13*10⁻²; p-values = 4.70*10⁻³ and 7.83*10⁻²; R ² = 0.29 and 0.10, respectively). Both the dorsal and lateral processes contribute to boosting rolling stability by increasing rotational inertia (with a positive effect on ΔI_x; Fig. 6C,F), albeit at the expense of agility and rapid maneuvering. This is because higher rotational inertia demands increased energy for rolling and results in greater resistance to rotation. Frontal processes contribute to drag reduction in certain species (Fig. 6A), particularly when their straight and sharp orientation effectively streamlines the body shape, facilitating smoother flow over the surface and reducing flow separation (e.g., Phialaspis, Unarkaspis, Europrotaspis, Panamintaspis, Pteraspis, or Errivaspis). However, the lack of a consistent effect of frontal process area on ΔD and the wide dispersion in ΔD values (Fig. 6A) indicate that these structures could serve diverse functions across species, extending beyond hydrodynamics. These functions might encompass predator deterrence, accommodation of sensory organs (Voichyshyn 2006), or adaptations to specific feeding habits (Tarlo 1961, 1962; Janvier 1996). This could be especially applicable to species featuring upward-pointing rostral processes that result in a higher amount of drag (e.g., Eglonaspis, Angaraspis, Olbiaspis, Larnovaspis, or Kureykaspis).

Our discoveries underscore the intricate interplay between form and function throughout the evolution of headshield processes, making a substantial contribution to the ongoing debate regarding the ancestral role of ostracoderm headshields. Several non–mutually exclusive hypotheses have been proposed, including defense against predators and ectoparasites and abrasion resistance, osmoregulation, and hydrodynamics (Mark-Kurik 1992). This work not only highlights the primary role these processes played in hydrodynamics, an idea that has persisted untested for decades (Mark 1961; Aleyev 1976; Aleyev and Novitskaya 1983; Bendix-Almgreen 1986; Belles-Isles 1987; Afanassieva 1992; Mark-Kurik 1992; Novitskaya 2000; Moloshnikov 2001; Elliott et al. 2010), but also acknowledges instances where additional complementary or alternative functions might have been at play. Our findings yield further insights into the remarkable ecological diversity achieved by pteraspidomorphs, showcasing a wide spectrum of adaptations that facilitate passive flow control around the body during locomotion (Mark-Kurik 1992; Botella and Fariña 2008; Blieck and Elliott 2017; Ferrón et al. 2020, 2021). We have unveiled a diverse array of variations in the architecture and development of headshield processes, each aligned with distinct biomechanical performances. This strong correlation underscores the multifaceted functions associated with these variations, tailored to specific ecological roles. Of particular interest are the headshields featuring well-developed lateral and/or dorsal processes positioned at the rear. These adaptations appear finely tuned for maintaining a consistent swimming direction, hinting at a potential preference for continuous swimming—characteristics frequently observed in numerous pteraspids and psammosteids. In contrast, species lacking or possessing less pronounced processes seem to excel in agility and maneuverability, advantageous traits in more intricate habitats. These kinds of headshields are evident across astraspids, arandaspids, cyathaspids, amphiaspids, and traquairaspids. These outcomes dispel the notion that pteraspidomorphs and ostracoderms were confined to bottom-dwelling lifestyles with limited swimming abilities (Aleyev and Novitskaya 1983; White and Toombs 1983; Belles-Isles 1987; Mark-Kurik 1992; Ferrón et al. 2020). Instead, it is highly probable that many pteraspidomorphs actively interacted with the water column in different ways through the evolution of distinctive body plans, thus underpinning ecological diversification. Ultimately, these findings contribute to the growing body of evidence (Purnell 1995, 2001a,b, 2002; Ferrón and Botella 2017; Ferrón et al. 2018, 2020, 2021; Miyashita et al. 2021; Ferrón and Donoghue 2022) that challenges the conventional narrative of a unidirectional trend toward progressively active lifestyles and constrained ecological niches among ostracoderms that seeks to explain the origin and early diversification of jawed vertebrates within the framework of the new head hypothesis (Gans and Northcutt 1983; Northcutt and Gans 1983; Gans 1989, 1993; Northcutt 1996, 2005).

Conclusions

Through the application of geometric morphometrics, combined with qualitative assessments of phylomorphospaces and quantitative metrics for convergence, we present a compelling demonstration of the independent and repeated emergence of frontal, dorsal, and lateral headshield processes within Pteraspidomorphi. Employing CFD, we establish that these structures evolved in response to shared hydrodynamic challenges, addressing common functional demands. Lateral processes, for instance, are shown to enhance lift generation and hydrodynamic efficiency by increasing the lift-to-drag ratio. Meanwhile, dorsal processes potentially contribute to yawing stability, and lateral processes offer yawing and pitch stability. Often positioned backward to the animal’s center of mass, dorsal and lateral processes could have acted analogously to an aircraft’s empennage, serving as directional stabilizers. Additionally, these processes bolster roll stability by increasing the rotational inertia, albeit at the expense of maneuverability. In various cases, frontal processes are linked to drag reduction, but they may serve a broader array of functions. In summation, our outcomes underscore the substantial variation in the arrangement and development of distinct headshield processes across Pteraspidomorphi. This diversity is intrinsically tied to biomechanical considerations and, consequently, likely reflects functional and ecological diversity. The recurring instances of both morphological and functional convergence in the headshields and processes of heterostracans illuminate the intricate evolutionary pathways of lifestyles and ecologies exhibited by stem-gnathostomes. These findings collectively challenge the notion of a unidirectional trend toward increasingly active lifestyles and constrained ecologies among ostracoderms, as previously put forth to explain the origin and early diversification of jawed vertebrates within the framework of the new head hypothesis.