Morphological evolution during the last hurrah of the trilobites: morphometric analysis of the Devonian asteropyginid trilobites

Abstract

The Asteropyginae Delo, 1935 is a group of phacopid trilobites in the family Acastidae Delo, 1935 that has served as the focus for several studies due to their distinctive morphologies and diversity. However, despite an interest in these characteristic morphologies, there have been no studies that have examined this group using morphometric techniques. Our investigation utilized both geometric morphometric and elliptical Fourier methods to quantify the morphology of cephalic sclerites of asteropyginid specimens representing wide taxonomic sampling of the clade. We constructed a phylomorphospace that shows temporal and spatial patterns of phenotypic evolution within the framework of a novel tip-dated phylogenetic tree generated using Bayesian inference. We recovered similar patterns in disparity regardless of the morphometric approach. Both analyses illustrated a marked expansion into morphospace throughout the temporal range of the clade, peaking in disparity in the Emsian and with European taxa exhibiting the highest disparity in glabellar morphospace. Additionally, glabellar shape showed low phylogenetic signal and no major patterns in phylomorphospace. This study highlights the utility of employing different methodologies to quantitatively explore the disparity of fossil taxa. It also illustrates some of the patterns of morphological change occurring during one of the final and major evolutionary radiations within Phacopida.

Introduction

The Asteropyginae Delo, 1935 is a group of phacopid trilobites that famously exhibit some of the most remarkable morphologies of all Devonian trilobites (Bignon and Crônier 2013; Bignon et al. 2014). Members of this group show a high degree of phenotypic diversity, varying from the extraordinary pre-cephalic projections of Walliserops trifurcatus Morzadec, 2001 to the relatively simple forms seen in Greenops boothi (Green, 1837). It is important to recognize that disparity has been referred to and has been used in a variety of different ways by many different authors (e.g., Gould 1991; Smith and Lieberman 1999; Wills 2001; Crônier 2013; Hopkins and Gerber 2017; Guillerme et al. 2020). For the purpose of this study, disparity will be described as the measure of morphological variation among species within the Asteropyginae. The Asteropyginae itself is essentially cosmopolitan and composed of more than 250 species (Harrington et al. 1959; Lieberman and Kloc 1997; Bignon and Crônier 2013; Bignon et al. 2014), typically treated as ranging from the Lochkovian to the Frasnian (Feist 1991; Morzadec 1992), though at least one species, Asteropyge gdoumontensis Asselberghs, 1930, occurs in the Pridolian (Thomas et al. 1984; Van Viersen and Prescher 2009; Storey 2012). Thus, the Asteropyginae comprises a singular example of one of the final major evolutionary radiations within Trilobita. This radiation occurred during and appears to have been profoundly influenced by a time of major climatic oscillations, with concomitant effects on sea level, as well as substantial tectonic reorganization. It also was associated with substantial biogeographic shifts (Abe and Lieberman 2009; Holloway and Rustán 2012; Carbonaro et al. 2018; Dowding and Ebach 2019; Penn-Clarke 2019; Bault et al. 2022b).

The Asteropyginae has not previously been subjected to quantitative morphometric analysis despite broad interest in its distinctive morphologies and the fact that it has been incorporated in several phylogenetic studies, including Lieberman and Kloc (1997), Bignon and Crônier (2013), and Bignon et al. (2014), with their broader phylogenetic placement also considered by Edgecombe (1993). Here we utilize both landmark-based geometric morphometric (GM; for applications to trilobites, see Smith and Lieberman 1999; Sheets et al. 2004; Crônier et al. 2005, 2015; Crônier and Fortey 2006; Hopkins and Webster 2009; Abe and Lieberman 2012; Bignon and Crônier 2012; Hopkins and Pearson 2016; Álvaro et al. 2018) and outline-based (Crônier et al. 1998, 2005; López Carranza and Carlson 2021)/elliptical Fourier (EF) methods (for applications to trilobites and their close relatives, see Foote 1989; Crônier et al. 1998, 2005; Hopkins 2014; Jackson and Budd 2017) in order to consider and also compare and contrast different types of quantitative information on the morphology of cephalic sclerites. Cephala have been frequently used in landmark-based analyses of trilobite morphology (e.g., Hughes 1994; Smith and Lieberman 1999; Adrain 2005; Crônier et al. 2005, 2015; Crônier and Fortey 2006; Webster and Zelditch 2005; Webber and Hunda 2007; Hopkins 2011, 2017; Abe and Lieberman 2012; Monti 2018; Webster and Sundberg 2020), because they contain abundant character information often used to define species and higher taxonomic categories, are abundantly preserved, and are considered a valuable repository of information on overall trilobite shape. We also consider how cephalic disparity changes through time and varies by biogeographic region in order to assess the relationship between what are anecdotally considered to be times or regions possessing highly distinctive morphologies and actual calculated values. Further, we place this quantitative information on morphology into the context of a phylogenetic hypothesis for the group, which builds on previous perspectives, to view changes in morphology in relation to the diversification of the clade. To visualize the disparity of the group, we construct a phylomorphospace that shows temporal and spatial patterns of phenotypic evolution using an updated tip-dated phylogenetic tree that was generated using Bayesian inference.

Materials and Methods

Specimens

Specimens were included in order to represent a wide taxonomic sampling of the Asteropyginae. Landmark data were captured from photographs of previously published cephala in dorsal view belonging to 64 species across 37 genera within the Asteropyginae (Supplementary Table 1). Of the photographs used, 41 of the images were primary type specimens (i.e., holotype, paratype, neotype) and six species are represented by multiple specimens (see Supplementary Table 1) to test whether individuals of the same species cluster in morphospace. Only specimens that were sufficiently complete for all landmarks considered, that were lacking in visible signs of deformation, and that appeared to be in the proper and correct orientation were included in the analysis as described in Shaw (1957). Any images of specimens that obviously differed from such an orientation were excluded from the present study. However, it is possible that various photographers have slightly different conceptions or perceptions of what comprises standard orientation (or even that a single photographer's perception could vary over time or from one specimen to another). This could cause subtle differences between species to be introduced artifactually. Given the range of morphological differences observed across the clade, this is unlikely to be playing a large role in the patterns observed. Glabella were the specific focus of this study, because other important parts of the cephalon such as lateral, posterior, and anterior borders, facial sutures, and precise position of the ocular lobes were not consistently preserved and not always discernible across specimens.

Note that while the pygidial sclerites of members of this group can be highly recognizable and do contain a wealth of diagnostic information, detailed investigation and analysis using landmark and outline approaches conducted on asteropyginids revealed a degree of taphonomic distortion in several individuals that skewed results of GM analyses. Although current retrodeformation methods exist that resolve problems with deformed specimens (Motani 1997), the images that served as a basis for measurements in this study were sourced from various digital resources and thus lacked a common scale for deformation. Additionally, we assessed images of pygidia that contained more variation in orientation than those of the cephala. For these reasons, the results using pygidia are not presented. Instead, cephala were the focus of the analysis.

Geometric Morphometrics

Geometric morphometric approaches were used to describe the variation of asteropyginid trilobite dorsal glabellar morphology. Landmarks were digitally placed on previously published photographs of the dorsal side of 62 species across 36 genera. Each glabella was initially delineated by nine homologous landmarks and 20 sliding semilandmarks (Fig. 1, blue) situated along five curves (Table 1). Landmarks and curves were placed preferentially on the left side of the glabella, because this side was the most consistently preserved for the available specimens. However, if the left side was damaged or incompletely preserved, then the image was flipped horizontally, and the reflected right side was used. These homologous landmarks and curves were digitally placed in R (R Core Team 2020), following similar positions used in other morphometric studies of trilobites (e.g., Crônier et al. 2004, 2015; Sheets et al. 2004; Webster and Sheets 2010; Abe and Lieberman 2012) using the package Stereomorph v. 1.6.4 (Olsen and Westneat 2015). To fully describe the bilaterally symmetrical shape of the glabella, these landmarks and curves, with the exception of midplane landmarks, were vertically mirrored to produce a complete configuration (see Cardini 2016). These reflected landmarks and semilandmarks are shown in yellow in Figure 1. Descriptions of the landmarks and semilandmarks used in this study are presented in Table 1. A general Procrustes analysis (GPA) was performed on the complete configuration (left side, midplane, and mirrored landmarks and semilandmarks) using the R package geomorph v. 4.0.0 (Adams et al. 2021). Because both landmarks and semilandmarks were analyzed, minimum bending energy was used to superimpose semilandmarks. Bending energy optimization minimizes the thin-plate spline bending energy between specimens, as opposed to Procrustes distance optimization, which minimizes the Euclidian distance between specimen shape coordinates.

Figure 1. Landmark and semilandmark scheme used in this study. In blue, landmarks (dark blue) and semilandmarks (light blue) initially digitized on specimen images. In yellow, landmarks (dark yellow) and semilandmarks (light yellow) generated by mirroring initial configuration. See Table 1 for further descriptions.

Table 1. Landmark and semilandmark descriptions.

Outline Analysis

Glabellar outlines were manually digitized from dorsal photographs belonging to 46 species across 33 genera using Adobe Illustrator 2020. To neutralize confounding variables related to asymmetry, points were digitized for each specimen from half of the glabella, depending on which half was better preserved in each specimen. Coordinate configurations (Fig. 1) were then mirrored across the medial plane so that all outline configurations represent the whole glabella (i.e., had left and right halves). Outline coordinates were then automatically extracted using the R package Momocs v. 1.3.3 (Bonhomme et al. 2014). To normalize (i.e., eliminate variation due to size, rotation, and translation) the outline coordinates before the elliptical Fourier analysis (EFA), a GPA was performed using three control landmarks: (1) anterior-most point of the glabella in the sagittal line, (2) interior-most point of right S3, and (3) interior-most point of left S3. Once outline coordinates were aligned, an EFA was performed using 99% of harmonic power (Bonhomme et al. 2014; López Carranza and Carlson 2021). The resulting morphological variables (EF coefficients) were then used as input for a principal component analysis (PCA).

Calculating Variation and Patterns of Disparity by Time and Geography

The PCAs for each morphometric analysis were plotted using the gm.prcomp function in geomorph (Adams et al. 2021). To visualize patterns of shape change among asteropyginids by geography and time, taxa in the PCA were colored based on geographic locality or time of first appearance. Geographic and stratigraphic data were obtained via a comprehensive review of the literature and examination of data housed in the Paleobiology Database (https://paleobiodb.org) and Integrated Digitized Biocollections (https://www.idigbio.org). Stratigraphic resolution was to the level of stage. Some studies, such as that of Bignon and Crônier (2013), were able to use greater stratigraphic resolution, especially for the Emsian interval, but because we were considering taxa distributed across several continents, we were not able to subdivide the Emsian or other stages to a finer scale. Geographic areas used generally correspond to broad areas of endemism and biogeographic regions identified in numerous previous studies of Devonian paleogeography (e.g., Lieberman and Eldredge 1996; Scotese et al. 1999; Rode and Lieberman 2005; Carrera and Rustán 2015; Dowding and Ebach 2018, 2019; Bault et al. 2022b). However, we refer to the region “Europe” only for heuristic purposes and in reference to the modern day (involving fossils distributed in what are today parts of western and central Europe, including Belgium, France, Germany, and Spain), as in the Siluro-Devonian these comprised a set of different terranes. Trilobite taxa may have moved within and between distinct regions at this time due to changes in sea level, tectonic collisions, or via dispersal (Lieberman and Eldredge 1996; Bault et al. 2022b). Data for species’ geography and first appearance are presented in Supplementary Table 1. To compare differences in the amount of glabellar variation between taxa by geologic time and by geography, we calculated Procrustes variances using the function morphol.disparity from the geomorph package (Adams et al. 2021) using a procD.lm model following procedures similar to those used in Friedman et al. (2019) and Martin et al. (2022). Unless otherwise noted, p-values were based on 10,000-iteration permutations.

Phylogeny

The character–taxon matrix developed by Bignon and Crônier (2013) for phylogenetic analysis of the Asteropyginae was used, augmented with the two additional taxa incorporated by Bignon et. al (2014). To confirm the results of previous studies, a heuristic search in PAUP*4.0a build 169 (Swofford 2003) was performed to find the most parsimonious cladogram, with branch swapping performed using tree bisection reconnection and taxa added by random sequence addition with 100 replicates. The recovered topology of the resultant most parsimonious tree aligned very closely with the results of Bignon and Crônier (2013) and Bignon et al. (2014). The matrix was then loaded into BEAUti v. 2.6.4 (Bouckaert et al. 2014) for parameterization, then used in a Bayesian tree search conducted using BEAST v. 2.6.4 (Bouckaert et al. 2014) employing the Mk model (Lewis 2001) of morphological trait evolution. Characters were subject to default partitioning based on the number of states. Bayes factors (Kass and Raftery 1995) were calculated for each of four different clock models, with marginal-likelihood estimates generated from path sampling. In this case, the uncorrelated exponential relaxed clock was found to be the preferred model and was used in the final tree search. This model allows rates to vary independently across branches with values drawn from an exponential distribution, the variance of which is estimated from the data and for each partition of character states.

Tip dates were specified based on stratigraphic occurrence data collected from the primary literature for each taxon. Stages were correlated with the 2012 timescales for the Silurian and Devonian (Gradstein et al. 2012) and subsequent updates to the International Chronostratigraphic Chart (Cohen et al. 2013) in order to obtain numerical tip dates that were used in the Bayesian tree search. This method of tip dating allows resolution to the lowest international stage boundary for each taxon and has been used in previous studies (e.g., Paterson et al. 2019).

The tree prior chosen for this analysis is the fossilized birth–death (FBD) model. This model describes the probability of the tree topology and fossils given a set of parameters. It is a stochastic branching model with parameters that describe speciation rate (λ), extinction rate (μ), fossil recovery rate (ψ), and the probability of sampling extant species (ρ) (Stadler 2010; Heath et al. 2014). Gavryushkina et al. (2014) modified this model in the SA (“Sampled Ancestors”) package in BEAST2, allowing a version of FBD to be used as a tree prior with new parameters for Markov chain Monte Carlo (MCMC) optimization: net diversification rate (λ − μ), turnover (μ/λ), and sampling proportion (ψ/(μ + ψ)). The ρ parameter was excluded because there are no extant representatives of this clade. A constraint in-group Asteropyginae monophyly was applied, which has been confirmed by multiple studies (e.g., Lieberman and Kloc 1997; Bignon and Crônier 2013). MCMC analyses consisted of independent runs sampling every 2500 generations for 50 million generations per run with a burn-in of 25%. A total of four runs was sufficient for convergence, which was assessed by visual confirmation of log-likelihood plots in Tracer v. 1.7.1 (Rambaut et al. 2018) and effective sample sizes greater than 200. Tree Annotator v. 2.6.4 (Bouckaert et al. 2014) was used to generate a maximum clade credibility (MCC) tree. Each clade within the tree was given a score based on its frequency within the sampled posterior trees, and the product of these scores within a tree is its score. The tree with the highest score is the MCC tree—a fully resolved tree that summarizes the posterior distribution of tree topologies. Parsimony and Bayesian approaches to phylogenetic analysis use different methods of optimization, and the latter is inherently a type of statistical phylogenetic analysis; for additional discussion of this, the reader is referred to Wiley and Lieberman (2011) and other relevant references on the topic. The characters used in the phylogenetic analysis include an array of cephalic and pygidial characters that are not identical to nor do they correspond in a one-to-one way with the landmark and outline data and are not considered in the morphological portion of this study. However, both are considered aspects of trilobite morphology and thus are only quasi-independent.

Analysis of Phylogenetic Signal

To evaluate patterns of glabellar shape change across the phylogenetic history of Asteropyginae, a phylomorphospace (Sidlauskas 2008) was created using the PCA from each morphometric analysis and the R package phytools (Revell 2012). Phylomorphospaces were plotted using a calculated average location of each genus in morphospace. Paraphyletic genera remained unaveraged, and for these genera, each species in the genus was incorporated. Tips were associated with species present in both the tree and morphometric datasets. Phylomorphospaces were generated and the trilobite phylogeny inferred in this study. Then, phylogenetic signal was analyzed for the landmark (GM) data and our tree using the function physignal from the package geomorph. The resulting K-statistic was compared with a null distribution generated from permutation tests using the average shapes of species or genera. Species or genera not present in the GM study were trimmed from the tree, and any species present in the GM study but not in the tree were removed from the dataset.

Results

Principal Component Analyses of Landmarks and Outlines

For the landmark-based (GM) analysis (Fig. 2), principal component (PC) 1 and PC 2 account for 38.98% and 29.97% of the total variance, respectively (Fig. 2A). The consensus landmark configuration (Fig. 2B) depicts the average glabellar morphology as black points. The gray points in Figure 2B represent the Procrustes-aligned individual landmark configurations. Variation along PC 1 mainly relates to the lateral narrowing of the anterior margin. Shape change along PC axes is shown in Figure 2C. Negative PC 1 scores correspond to rounder anterior margins, while positive PC 1 scores correspond to more laterally constrained, pinched margins. Similarly, negative PC 1 scores correlate with more curved lateral glabellar margins compared with straighter margins in positive PC 1 scores. Shape change along PC 2 is primarily characterized by an outward curving of the lateral margin of L3, as observed in positive PC 2 scores. Furthermore, positive PC 2 scores correspond to a narrower anterior lobe. PC 3 accounts for 7.53% of shape variation and shows an axis of variation associated with an overall narrowing of the glabella (Supplementary Fig. 2).

Figure 2. A, Asteropyginid specimens shown in tangent space of principal component (PC) 1 and PC 2 resulting from the geometric morphometric analysis, with colors corresponding to genera. B, Consensus configuration from the Procrustes analysis. Black dots depict positions of landmarks and semilandmarks; gray dots depict the variation in specimen landmark locations around the average. C, Deformation grids showing glabellar shape at the extremes of PC 1 and PC 2.

Glabellar morphology within some genera can be diverse (e.g., Alcaldops Arbizu, 1979, Bellacartwrightia Lieberman and Kloc, 1997, Bradocryphaeus Haas and Mensink, 1969, Neocalmonia Pillet, 1969), with species moderately spread across shape space. Some genera, however, cluster in morphospace (e.g., Hollardops Morzadec, 1977 or Rheingoldium Basse, 2003), indicating less-variable morphologies (Fig. 2A, see also Supplementary Fig. 1).

For the outline-based (EF) analysis (Fig. 3), PC 1 and PC 2 account for 53.65% and 16.51% of the total variance, respectively (Fig. 3A). The mean outline shape is depicted in Figure 3B. Morphological variation along PC 1 corresponds to a medial narrowing, resulting in an overall lateral broadening. The anterior glabellar lobe widens from negative to positive on the PC 1 axis, and the lateral margins become narrower, which is also related to how S3 intersects with the axial furrow. The anterior margin of the glabella also becomes more tapered (Fig. 3C). Along PC 2, shape change is dominated by an overall narrowing of the glabella. Negative PC 2 scores correspond to wider glabellae with rounder frontal lobes, both anteriorly and laterally (Fig. 3C). Positive PC 2 scores relate to narrower, more elongated glabellar shapes with more pointed anterior lobes. PC 3 accounts for 13.15% of the total variance, and shape variation along this axis is associated with a narrowing of the anterior lobe and an overall posterior widening (L1–L3 and posterior margin; Supplementary Fig. 4).

Figure 3. Principal component analysis (PCA) of glabellar outlines described by elliptical Fourier coefficients and changes in outline shape corresponding to PCs. A, PC 1–PC 2; B, mean outline shape; and C, thin-plate spline deformation grids associated with PC 1 and PC 2.

Overall, there is minimal overlap among specimens, including specimens of the same genus (Fig. 3A). Similar to the results from landmark analysis, some genera show moderate spread in morphospace (e.g., Alcaldops, Bellacartwrightia, Minicryphaeus Bignon and Crônier, 2013). However, there does not appear to be significant clustering of species within genera (Fig. 3A; see also Supplementary Fig. 3).

Patterns of Disparity by Time and Geography

Results from the landmark dataset analyzing disparity in glabellar shape by time (Fig. 4A, Table 2) show Emsian taxa with the highest disparity in glabellar shape (0.0065532) and Givetian taxa with the lowest (0.0041207). Emsian taxa were close to significantly higher in their relative glabellar shape disparity compared with Givetian taxa (p = 0.062) and the upper Silurian taxon (disparity 0, p = 0.086; Table 2). Note though, the upper Silurian is only represented by one out-group taxon and that it is not necessarily valid to think of a taxon's disparity if it consists of a single individual. Figure 4A shows centralized distributions by time, with taxa from the Givetian additionally present in the positive extreme of PC 1 and negative extreme of PC 2 and taxa from the Frasnian generally localized around the negative values of PC 1. Although there are differences in location within morphospace, the spread across morphospace of taxa by each time is approximately equal. Please refer to Supplementary Figure 5 for landmark PCA panels without overlapping-time convex hulls.

Figure 4. Principal component analysis (PCA) of landmark data from the geometric morphometric analysis. Specimens and convex hulls colored by time (A) and geography (B).

Table 2. Disparity and p-values for both geography and time resulting from the shape disparity analysis on the landmark data.

Glabellar shape analysis by outline and time shows a pattern similar to that of the landmark data, with Emsian taxa having the highest disparity across morphospace (Fig. 5A). Groups by time are more separated in morphospace compared with their locations using landmark data. In addition to the Frasnian taxa, Eifelian taxa are similarly spread toward the negative end of PC 1 (Fig. 5A). The outline data also suggest a larger spread across PC 1 with slightly more restriction across PC 2, with the exception of taxa from the Eifelian. Please refer to Supplementary Figure 6 for outline PCA panels without overlapping-time convex hulls.

Figure 5. Principal component analysis (PCA) of outline data from the elliptical Fourier analysis. Specimens and convex hulls colored by time (A) and geography (B).

Assessing disparity of glabellar shape by geography on the landmark dataset (Fig. 4B; Table 2) shows that European taxa have the highest disparity (0.0082629) and eastern North American taxa have the lowest disparity (0.0034837) relative to taxa found in other regions in this dataset. European taxa were significantly higher in their glabellar disparity compared with North American taxa (p = 0.013), but there were no other significant differences in disparity between geographic groups. There is no distinct pattern in glabellar shape disparity by geography and geographic region, except that eastern North American and western Asian taxa are more constrained in morphospace and have limited spread toward the positive end of PC 2. Additionally, eastern North American taxa do trend slightly toward the positive axes of PC 1 (Fig. 4B).

Similar to the landmark dataset, the glabellar shape analysis using outlines also shows high disparity in shape for European taxa (Fig. 5B, light blue area) and similar spread across morphospace by geographic region. Each geographic group is equally spread across PC 1 compared with PC 2, where they are much more constrained, except European taxa (Fig. 5B).

It is important to note that the nature of biogeographic regions does change through time, and we are not implying that all these regions have entirely monophyletic area relationships throughout the study interval.

Phylogeny

The results of the phylogenetic analysis (Fig. 6) largely agree with those of previous analyses on this and related datasets (e.g., Lieberman and Kloc 1997; Bignon and Crônier 2013; Bignon et al. 2014). While similar patterns of constituent clades are recovered, there are two distinct and moderately well-supported groups similar to those seen in Lieberman and Kloc (1997), Bignon and Crônier (2013), and Bignon et al. (2014). The first of these groups is recovered at a node with a posterior probability of 0.57 and contains several well-supported clades. The genera Greenops Delo, 1935, Stummiana Lieberman and Kloc, 1997, Deloops Lieberman and Kloc, 1997, and Breizhops Morzadec, 1983 are sister to the Minicryphaeus Bignon and Crônier, 2013–Destombesina Morzadec, 1997 clade with a posterior probability of 0.41. Sister to that group is a clade composed of the genera Kayserops Delo, 1935 and representative species of Rhenops Richter and Richter, 1943, Paracryphaeus Gandl, 1972, Gandlops Bignon and Crônier, 2013, Rheingoldium, and Braunops Lieberman and Kloc, 1997, joined at a node that is moderately well supported (posterior probability of 0.57). The second major grouping is a well-supported clade (posterior probability 0.64) of the remaining asteropyginids.

Figure 6. Bayesian maximum clade credibility (MCC) tree. Node symbols indicate posterior support, with darker values indicative of stronger support and lighter values showing weaker support. Posterior probability values provided in the Supplementary Material.

Relationships similar to those reported in Bignon et al. (2014) are recovered in the clade of Minicryphaeus; however, in contrast to the results from Bignon et al. (2014), Destombesina is not recovered as a sister to the Asteropyginae. Conversely, it is well supported (posterior probability of 0.96) as a sister to the Minicryphaeus clade and includes the species Ganetops ebbae (Richter and Richter, 1954) as a basal taxon (Bignon et al. 2014).

Bignon and Crônier (2013) retrieved a monophyletic group that they described as reflecting the “progressive type” of Struve (1959), and it was also recovered in the Bayesian analysis (Fig. 6), with strong support (posterior probability of 0.97). Instead of a grade topology recovered by previous analyses, the Alcadops clade is resolved sister to a Radiopyge Farsan, 1981–Neocalmonia Pillet, 1969 clade with a posterior probability of 0.99.

Although there appears to be reasonable support for important aspects of the topology and broad congruence with previous results, information regarding the timing of diversification events (Fig. 6) differs from the conclusions of some previous studies. For instance, some of the early evolutionary divergence within Asteropyginae is projected to have occurred well back into the Ordovician and Silurian, whereas several previous interpretations (Delo 1935; Struve 1959; Gandl 1972; Arbizu 1979; Morzadec 1983; Lieberman and Kloc 1997; Bignon and Crônier 2013; Bignon et al. 2014) have largely treated these diversification events as occurring within the Devonian and perhaps the late Silurian. However, by the same token, the results also suggest that more work may be needed on understanding the classification and origins of the subfamily. This aligns with the conclusions of Ramsköld and Edgecombe (1993) and Chatterton et al. (2006) regarding the status of some of the subfamilies within the Acastidae.

Phylomorphospace and Phylogenetic Signal

The phylomorphospace plots (Figs. 7, 8) show an array of overlapping branches with no strong pattern of separation by clade. The phylogenetic signal K-statistic on the landmark data using the phylogeny inferred here was K = 0.466, p < 0.001. A K-value this low suggests extremely low phylogenetic signal in the glabella landmark data.

Figure 7. Phylomorphospace plot of principal component (PC) 1 and PC 2 from the landmark data resulting from the geometric morphometric analysis. Circle positions represent the average location of genera in morphospace, with paraphyletic genera not averaged. Please refer to Fig. 6 for our hypothesized tree topology.

Figure 8. Phylomorphospace plot of principal component (PC) 1 and PC 2 from the outline data resulting from the elliptical Fourier analysis. Circle positions represent the average location of genera in morphospace, with paraphyletic genera not averaged. Please refer to Fig. 6 for our hypothesized tree topology.

Discussion

This study focused on using and comparing quantitative morphometric analyses to explore the phenotypic diversity of the Asteropyginae with the incorporation of phylogenetic, geographic, and temporal information. The asteropyginids were undergoing evolutionary radiation toward the end of the evolutionary history of one of the most diverse and distinctive orders of Trilobita, the Phacopida. As such, the Asteropyginae comprises one of the last gasps or hurrahs of diversification within Trilobita as a whole. We included 64 of the more than 250 species and 38 genera, recovering a range of variation in asteropygine glabellar shape and furrow position. For decades, geometric morphometrics and other landmark-based methods have been used to better understand and quantify shape variation across a wide range of study systems. Newer analytical techniques that incorporate curve fitting and additional shape information into analyses have become increasingly popular to address questions regarding shape variation, especially in organisms lacking an abundance of homologous landmarks. With the increasing use of EFA, there are a growing number of studies that perform a comparison of quantitative methods on the same dataset (e.g., Loy et al. 2000; Russo et al. 2007; Van Bocxlaer and Schultheiß 2010; Dujardin et al. 2014).

Our study compares the analytical outcomes of geometric morphometrics and EFA and recovers similar patterns of glabellar variation between morphospaces resulting from these PCAs. Similar to previous studies comparing these methods, we find that regardless of the morphometric approach, there are similar and distinct axes of variation (Figs. 2C, 3C). From negative to positive values along the PC 1 axis, there is a trend of lateral narrowing of the anterior margins in both analyses. Differences in PC 2 between analyses are more apparent, as GM methods pick up a distinct outward curving of the lateral margin of L3 (29.97% of shape variation), whereas EFA methods highlight the narrowing of the glabella (16.51% of shape variation). Analysis of PC 3 from the GM methods suggests axes of variation similar to that of PC 2 in the EFA methods (Supplementary Fig. 2), which accounts for 9.66% of overall shape variation. Results from the EFA methods also show this trend within PC 3 (13.15% of shape variation), reflecting an axis of shape variation similar to that produced using GM methods PC 2 (Supplementary Fig. 4). The slight discrepancies between these analyses may highlight the differences of including interior landmarking in the GM analysis in addition to using semilandmark curves versus just using glabellar outlines in the EFA analysis. Overall, there are similar trends across taxa among specimens in the EFA analysis compared with the GM analysis and specimens of the same genus, possibly indicating that these analyses are similar in distinguishing subtle morphological differences in the glabella, at least in this clade (Fig. 3). In essence, both methods seem to be highly effective at quantifying morphology, as Crônier et al. (2005) previously demonstrated for a different trilobite clade. Further, this suggests that even in a system such as trilobites, where there are numerous homologous features that can be utilized for GM analysis, incorporating information on outlines using EFA can prove beneficial.

This study suggests the Emsian represented the time of maximal asteropyginid disparity, which was also the apogee of the clade's diversity (Bignon and Crônier 2013). This study also suggests that although taxa occurring in the Emsian are most disparate in their glabellar morphology, there are occurrences of taxa in novel areas of morphospace (e.g., Eifelian, Frasnian, Pragian) not occupied by Emsian taxa (Figs. 4, 5). Equivalently, the European biogeographic region encompassed the clade's maximal disparity and diversity, with taxa in other biogeographic regions not occupying novel areas in morphospace.

The overall patterns of increasing disparity do not align well with what might be predicted during an adaptive radiation (Simões et al. 2016). In particular, there are no signs of an early burst or of clustering or partitioning of the different parts of morphospace, nor do certain morphological regions of the glabella act as attractors. Instead, morphospace occupied seems to broadly diffuse outward. The morphometric data are also not displaying any prominent phylogenetic signal. Several other clades of radiating Devonian trilobites show related patterns (Eldredge and Cracraft 1980; Lieberman 1993; Abe and Lieberman 2012; Carbonaro et al. 2018), though more data are needed to discern the extent to which the general model for Devonian trilobites in general, and asteropyginids in particular, is evolutionary radiation without adaptive radiation. Previous work by Fortey and Owens (1999) postulated how cephalon and glabellar disparity is tied to feeding habits in trilobites. They suggest that “advanced predators” trend toward expanded anterior glabellar lobes. The disparity in the anterior glabella we see in this study could be associated with changes in feeding strategy during trilobite evolution that may not be as tightly associated with phylogenetic history.

In total, the asteropyginids, in spite of their distinctive morphologies, seem to represent a clade that is drifting through glabellar morphospace, showing no concerted patterns, other than a steady increase in glabellar disparity as diversity increased, followed by a steady decline and then the blinking out of glabellar disparity as the diversity of the subfamily and family declined, with the entire order it belongs to eventually disappearing at the end of the Devonian (Feist 1991; Bignon and Crônier 2013; Crônier 2013; Van Viersen 2013; Bignon et al. 2014; Van Viersen and Bignon 2018; Van Viersen and Vanherle 2018; Bault et al. 2022a). A key question in any clade's history is the extent to which the changing patterns of disparity through time are related to patterns of speciation as opposed to extinction (Guillerme et al. 2020). One of the distinctive aspects of the Late Devonian biodiversity crisis in general is it seems to be produced not so much by an increase in extinction rate as a decline in speciation rate post-Eifelian or Givetian, with extinction rate remaining relatively constant throughout the Devonian (Bambach et al. 2004; Rode and Lieberman 2004). However, Morzadec (1992) argued that there was prominent diversification in the group during the Givetian–Frasnian. The initial expansion of disparity in asteropyginids may be associated with no biases in the production (speciation) or elimination (extinction) of morphology within the quantified morphospace, although it is worth noting that Bignon and Crônier (2014) highlighted how there was some specialization of asteropyginids into peri-reefal environments. In addition, the subsequent decline in disparity, if asteropyginids fit the general pattern of the Late Devonian biodiversity crisis, is likely attributable to continued extinction with limited diversification, and again there do not seem to be any particular biases in the morphology of forms that went extinct. To test this hypothesis in greater detail though, quantitative information on speciation and extinction rates in asteropyginids is needed. This will in turn provide further detail on how biogeographic and ecological factors contributed to the initial radiation, as well as the subsequent crisis (McGhee et al. 2004; Rode and Lieberman 2004; Abe and Lieberman 2009, 2012; Bault et al. 2022a).

There were a variety of interesting research questions that our study could not consider. For instance, several highly useful studies, including Smith (1998), Crônier et al. (1998, 2004, 2005, 2015), Crônier and Fortey (2006), Hunda and Hughes (2007), Webber and Hunda (2007), Hopkins (2011), and Bignon and Crônier (2012), have employed geometric morphometrics to consider and quantify levels of variation within species of trilobites, albeit from other clades. By contrast, a specific focus of our study was quantifying morphological change across the Asteropyginae and therefore, in most cases, morphometric data were not collected from multiple specimens within individual species. This meant that the degree of variation within individual species was not considered in detail in the present study, and thus it is possible that important patterns were missed. However, various phylogenetic studies, including Lieberman and Kloc (1997), Bignon and Crônier (2013), and Bignon et al. (2014), have suggested that levels of variation within individual asteropyginid species are low, and they have indicated low levels of polymorphism within species.