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 Reference Bignon and Crônier2013; Bignon et al. Reference Bignon, Corbacho and López-Soriano2014). Members of this group show a high degree of phenotypic diversity, varying from the extraordinary pre-cephalic projections of Walliserops trifurcatus Morzadec, Reference Morzadec2001 to the relatively simple forms seen in Greenops boothi (Green, Reference Green1837). 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 Reference Gould1991; Smith and Lieberman Reference Smith and Lieberman1999; Wills Reference Wills, Adrain, Edgecombe and Lieberman2001; Crônier Reference Crônier2013; Hopkins and Gerber Reference Hopkins, Gerber, Nuno de la Ruso and Müller2017; Guillerme et al. Reference Guillerme, Cooper, Brusatte, Davis, Jackson, Gerber, Goswami, Healy, Hopkins, Jones, Lloyd, O'Reilly, Pate, Puttick, Rayfield, Saupe, Sherratt, Slater, Weisbecker, Thomas and Donoghue2020). 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. Reference Harrington, Henningsmoen, Howell, Jaanusson, Lochman-Balk, Moore, Poulsen, Rasetti, Richter, Richter, Schmidt, Sdzuy, Struve, Stormer, Stubblefield, Tripp, Weller and Whittington1959; Lieberman and Kloc Reference Lieberman and Kloc1997; Bignon and Crônier Reference Bignon and Crônier2013; Bignon et al. Reference Bignon, Corbacho and López-Soriano2014), typically treated as ranging from the Lochkovian to the Frasnian (Feist Reference Feist1991; Morzadec Reference Morzadec1992), though at least one species, Asteropyge gdoumontensis Asselberghs, Reference Asselberghs1930, occurs in the Pridolian (Thomas et al. Reference Thomas, Owens and Rushton1984; Van Viersen and Prescher Reference Van Viersen and Prescher2009; Storey Reference Storey2012). 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 Reference Abe and Lieberman2009; Holloway and Rustán Reference Holloway and Rustán2012; Carbonaro et al. Reference Carbonaro, Langer, Nihei, de Souza Ferreira and Ghilardi2018; Dowding and Ebach Reference Dowding and Ebach2019; Penn-Clarke Reference Penn-Clarke2019; Bault et al. Reference Bault, Crônier and Bignon2022b).
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 (Reference Lieberman and Kloc1997), Bignon and Crônier (Reference Bignon and Crônier2013), and Bignon et al. (Reference Bignon, Corbacho and López-Soriano2014), with their broader phylogenetic placement also considered by Edgecombe (Reference Edgecombe1993). Here we utilize both landmark-based geometric morphometric (GM; for applications to trilobites, see Smith and Lieberman Reference Smith and Lieberman1999; Sheets et al. Reference Sheets, Kim, Mitchell and Elewa2004; Crônier et al. Reference Crônier, Auffray and Courville2005, Reference Crônier, Budil, Fatka and Laibl2015; Crônier and Fortey Reference Crônier and Fortey2006; Hopkins and Webster Reference Hopkins and Webster2009; Abe and Lieberman Reference Abe and Lieberman2012; Bignon and Crônier Reference Bignon and Crônier2012; Hopkins and Pearson Reference Hopkins and Pearson2016; Álvaro et al. Reference Álvaro, Esteve and Zamora2018) and outline-based (Crônier et al. Reference Crônier, Renaud, Feist and Auffray1998, Reference Crônier, Auffray and Courville2005; López Carranza and Carlson Reference López Carranza and Carlson2021)/elliptical Fourier (EF) methods (for applications to trilobites and their close relatives, see Foote Reference Foote1989; Crônier et al. Reference Crônier, Renaud, Feist and Auffray1998, Reference Crônier, Auffray and Courville2005; Hopkins Reference Hopkins2014; Jackson and Budd Reference Jackson and Budd2017) 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 Reference Hughes1994; Smith and Lieberman Reference Smith and Lieberman1999; Adrain Reference Adrain2005; Crônier et al. Reference Crônier, Auffray and Courville2005, Reference Crônier, Budil, Fatka and Laibl2015; Crônier and Fortey Reference Crônier and Fortey2006; Webster and Zelditch Reference Webster and Zelditch2005; Webber and Hunda Reference Webber and Hunda2007; Hopkins Reference Hopkins2011, Reference Hopkins2017; Abe and Lieberman Reference Abe and Lieberman2012; Monti Reference Monti2018; Webster and Sundberg Reference Webster and Sundberg2020), 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 (Reference Shaw1957). 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 Reference Motani1997), 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. Reference Crônier, Feist and Auffray2004, Reference Crônier, Budil, Fatka and Laibl2015; Sheets et al. Reference Sheets, Kim, Mitchell and Elewa2004; Webster and Sheets Reference Webster and Sheets2010; Abe and Lieberman Reference Abe and Lieberman2012) using the package Stereomorph v. 1.6.4 (Olsen and Westneat Reference Olsen and Westneat2015). 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 Reference Cardini2016). 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. Reference Adams, Collyer, Kaliontzopoulou and Baken2021). 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. Reference Bonhomme, Picq, Gaucherel and Claude2014). 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. Reference Bonhomme, Picq, Gaucherel and Claude2014; López Carranza and Carlson Reference López Carranza and Carlson2021). 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. Reference Adams, Collyer, Kaliontzopoulou and Baken2021). 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 (Reference Bignon and Crônier2013), 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 Reference Lieberman and Eldredge1996; Scotese et al. Reference Scotese, Boucot and McKerrow1999; Rode and Lieberman Reference Rode and Lieberman2005; Carrera and Rustán Reference Carrera and Rustán2015; Dowding and Ebach Reference Dowding and Ebach2018, Reference Dowding and Ebach2019; Bault et al. Reference Bault, Crônier and Bignon2022b). 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 Reference Lieberman and Eldredge1996; Bault et al. Reference Bault, Crônier and Bignon2022b). 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. Reference Adams, Collyer, Kaliontzopoulou and Baken2021) using a procD.lm model following procedures similar to those used in Friedman et al. (Reference Friedman, Martinez, Price and Wainwright2019) and Martin et al. (Reference Martin, Davis and Smith2022). Unless otherwise noted, p-values were based on 10,000-iteration permutations.
Phylogeny
The character–taxon matrix developed by Bignon and Crônier (Reference Bignon and Crônier2013) 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 Reference Swofford2003) 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 (Reference Bignon and Crônier2013) and Bignon et al. (Reference Bignon, Corbacho and López-Soriano2014). The matrix was then loaded into BEAUti v. 2.6.4 (Bouckaert et al. Reference Bouckaert, Heled, Kühnert, Vaughan, Wu, Xie, Suchard, Rambaut and Drummond2014) for parameterization, then used in a Bayesian tree search conducted using BEAST v. 2.6.4 (Bouckaert et al. Reference Bouckaert, Heled, Kühnert, Vaughan, Wu, Xie, Suchard, Rambaut and Drummond2014) employing the Mk model (Lewis Reference Lewis2001) of morphological trait evolution. Characters were subject to default partitioning based on the number of states. Bayes factors (Kass and Raftery Reference Kass and Raftery1995) 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. Reference Gradstein, Ogg and Hilgen2012) and subsequent updates to the International Chronostratigraphic Chart (Cohen et al. Reference Cohen, Finney, Gibbard and Fan2013) 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. Reference Paterson, Edgecombe and Lee2019).
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 Reference Stadler2010; Heath et al. Reference Heath, Huelsenbeck and Stadler2014). Gavryushkina et al. (Reference Gavryushkina, Welch, Stadler and Drummond2014) 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 Reference Lieberman and Kloc1997; Bignon and Crônier Reference Bignon and Crônier2013). 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. Reference Rambaut, Drummond, Xie, Baele and Suchard2018) and effective sample sizes greater than 200. Tree Annotator v. 2.6.4 (Bouckaert et al. Reference Bouckaert, Heled, Kühnert, Vaughan, Wu, Xie, Suchard, Rambaut and Drummond2014) 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 (Reference Wiley and Lieberman2011) 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 Reference Sidlauskas2008) was created using the PCA from each morphometric analysis and the R package phytools (Revell Reference Revell2012). 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, Reference Arbizu1979, Bellacartwrightia Lieberman and Kloc, Reference Lieberman and Kloc1997, Bradocryphaeus Haas and Mensink, Reference Haas and Mensink1969, Neocalmonia Pillet, Reference Pillet1969), with species moderately spread across shape space. Some genera, however, cluster in morphospace (e.g., Hollardops Morzadec, 1977 or Rheingoldium Basse, Reference Basse2003), 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, Reference Bignon and Crônier2013). 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 Reference Lieberman and Kloc1997; Bignon and Crônier Reference Bignon and Crônier2013; Bignon et al. Reference Bignon, Corbacho and López-Soriano2014). 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 (Reference Lieberman and Kloc1997), Bignon and Crônier (Reference Bignon and Crônier2013), and Bignon et al. (Reference Bignon, Corbacho and López-Soriano2014). 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, Reference Delo1935, Stummiana Lieberman and Kloc, Reference Lieberman and Kloc1997, Deloops Lieberman and Kloc, Reference Lieberman and Kloc1997, and Breizhops Morzadec, Reference Morzadec1983 are sister to the Minicryphaeus Bignon and Crônier, Reference Bignon and Crônier2013–Destombesina Morzadec, Reference Morzadec1997 clade with a posterior probability of 0.41. Sister to that group is a clade composed of the genera Kayserops Delo, Reference Delo1935 and representative species of Rhenops Richter and Richter, Reference Richter and Richter1943, Paracryphaeus Gandl, Reference Gandl1972, Gandlops Bignon and Crônier, Reference Bignon and Crônier2013, Rheingoldium, and Braunops Lieberman and Kloc, Reference Lieberman and Kloc1997, 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. (Reference Bignon, Corbacho and López-Soriano2014) are recovered in the clade of Minicryphaeus; however, in contrast to the results from Bignon et al. (Reference Bignon, Corbacho and López-Soriano2014), 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, Reference Richter and Richter1954) as a basal taxon (Bignon et al. Reference Bignon, Corbacho and López-Soriano2014).
Bignon and Crônier (Reference Bignon and Crônier2013) retrieved a monophyletic group that they described as reflecting the “progressive type” of Struve (Reference Struve1959), 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, Reference Farsan1981–Neocalmonia Pillet, Reference Pillet1969 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 Reference Delo1935; Struve Reference Struve1959; Gandl Reference Gandl1972; Arbizu Reference Arbizu1979; Morzadec Reference Morzadec1983; Lieberman and Kloc Reference Lieberman and Kloc1997; Bignon and Crônier Reference Bignon and Crônier2013; Bignon et al. Reference Bignon, Corbacho and López-Soriano2014) 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 (Reference Ramsköld and Edgecombe1993) and Chatterton et al. (Reference Chatterton, Fortey, Brett, Gibb and McKellar2006) 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. Reference Loy, Busilacchi, Costa, Ferlin and Cataudella2000; Russo et al. Reference Russo, Costa and Cataudella2007; Van Bocxlaer and Schultheiß Reference Van Bocxlaer and Schultheiß2010; Dujardin et al. Reference Dujardin, Kaba, Solano, Dupraz, McCoy and Jaramillo-O2014).
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. (Reference Crônier, Auffray and Courville2005) 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 Reference Bignon and Crônier2013). 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. Reference Simões, Alvarado, Breitkreuz, Baca, Cooper, Heins, Herzog and Lieberman2016). 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 Reference Eldredge and Cracraft1980; Lieberman Reference Lieberman1993; Abe and Lieberman Reference Abe and Lieberman2012; Carbonaro et al. Reference Carbonaro, Langer, Nihei, de Souza Ferreira and Ghilardi2018), 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 (Reference Fortey and Owens1999) 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 Reference Feist1991; Bignon and Crônier Reference Bignon and Crônier2013; Crônier Reference Crônier2013; Van Viersen Reference Viersen and P2013; Bignon et al. Reference Bignon, Corbacho and López-Soriano2014; Van Viersen and Bignon Reference Van Viersen and Vanherle2018; Van Viersen and Vanherle Reference Van Viersen and Bignon2018; Bault et al. Reference Bault, Balseiro, Monnet and Crônier2022a). 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. Reference Guillerme, Cooper, Brusatte, Davis, Jackson, Gerber, Goswami, Healy, Hopkins, Jones, Lloyd, O'Reilly, Pate, Puttick, Rayfield, Saupe, Sherratt, Slater, Weisbecker, Thomas and Donoghue2020). 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. Reference Bambach, Knoll and Wang2004; Rode and Lieberman Reference Rode and Lieberman2004). 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 (Reference Bignon and Crônier2014) 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. Reference McGhee, Sheehan, Bottjer and Droser2004; Rode and Lieberman Reference Rode and Lieberman2004; Abe and Lieberman Reference Abe and Lieberman2009, Reference Abe and Lieberman2012; Bault et al. Reference Bault, Balseiro, Monnet and Crônier2022a).
There were a variety of interesting research questions that our study could not consider. For instance, several highly useful studies, including Smith (Reference Smith1998), Crônier et al. (Reference Crônier, Renaud, Feist and Auffray1998, Reference Crônier, Feist and Auffray2004, Reference Crônier, Auffray and Courville2005, Reference Crônier, Budil, Fatka and Laibl2015), Crônier and Fortey (Reference Crônier and Fortey2006), Hunda and Hughes (Reference Hunda and Hughes2007), Webber and Hunda (Reference Webber and Hunda2007), Hopkins (Reference Hopkins2011), and Bignon and Crônier (Reference Bignon and Crônier2012), 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 (Reference Lieberman and Kloc1997), Bignon and Crônier (Reference Bignon and Crônier2013), and Bignon et al. (Reference Bignon, Corbacho and López-Soriano2014), have suggested that levels of variation within individual asteropyginid species are low, and they have indicated low levels of polymorphism within species.
Acknowledgments
We thank associate editor M. Hopkins, A. van Viersen, and one anonymous reviewer for their very helpful comments and suggestions on previous versions of this article. We additionally thank L. Boucher, University of Texas, Austin, for assistance obtaining images of specimens. Financial support for this project was provided by NSF DBI 1602067 and the Biodiversity Institute, University of Kansas. Open access publishing of this article supported by the David Henry Wenrich Memorial endowment fund at the University of Kansas.
Declaration of Competing Interests
The authors declare no competing interests.
Data Availability Statement
All relevant data files and scripts for the analyses presented here are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.9w0vt4bj8.
Open access
