Introduction

Echinoids are a highly diverse group of marine invertebrates known from a plethora of modern and ancient habitats (Smith, 1984; Smith and Jeffery, 2000; Mooi, 2001; Kroh and Mooi, 2018; Arachchige et al., 2019). The group experienced a sharp diversity decline at the end of the Paleozoic, when the majority of stem-group echinoids went extinct after the end-Permian mass extinction (Thompson et al., 2018). However, their diversity rebounded in the Mesozoic era with the appearance of irregular echinoids (Saucède et al., 2007). Echinoids were abundant denizens of the shallow benthic communities of the southern parts of the Western Interior Seaway (WIS), in what is today the southwestern United States (e.g., Texas and New Mexico) and northern Mexico (Sonora and Coahuila). This is where the WIS opened into the much larger North Atlantic/West Tethys for much of its duration (Blakey, 2014). Many of the first echinoids found in this region were described in various nineteenth-century studies (e.g., Roemer, 1849, 1852; Clark, 1893). Subsequent works include descriptions by Clark (1915) and Cooke (1946, 1953, 1955) and publications focused on taxa from Texas (Adkins, 1928; Whitney and Kellum, 1966; Smith and Rader, 2009; Morgan, 2016; Thompson, 2016). Echinoids have been frequent subjects of phylogenetic analysis (relatively recent publications in this area include Kroh and Smith, 2010; Mooi et al., 2014; Koch et al., 2018; Thompson et al., 2018).

Among the irregular echinoids, the most diverse and abundant are the heart urchins of the order Spatangoida, which inhabit all oceans and a wide array of habitats (Stockley et al., 2005). Spatangoids are deposit-feeding burrowers that radiated in the Cretaceous and survive into the present (Mortensen, 1950; Fischer, 1966; Smith, 1984; Mooi, 2001; Smith and Wright, 2008, 2012; Smith and Kroh, 2011; Kroh and Mooi, 2018). They have also been incorporated into phylogenetic analyses that considered both living and fossil taxa (Villier et al., 2004; Stockley et al., 2005; Saucède et al., 2007; Kroh and Smith, 2010; Kroh et al., 2014). The goal of the present study is to focus on species-level phylogenetic relationships of Cretaceous spatangoids, with special emphasis on taxa from the WIS, although it was not possible to consider to any great extent the diverse range of taxa present in the Mississippi Embayment, an important and related biogeographic region (see Zachos, 2017), as taxa in this region primarily radiated too late to be relevant to the present analysis. The present study also uses the phylogenetic perspective to consider aspects of the biogeographic relationship between the WIS and other parts of the Cretaceous world.

Table 1.

Taxa considered in phylogenetic analysis, relevant material examined, and geographic occurrence; KUMIP = Division of Invertebrate Paleontology, Biodiversity Institute, University of Kansas, Lawrence, USA; NMMNH = New Mexico Museum of Natural History, Albuquerque, USA; STIB = Universität Bonn, Germany; USNM = United States National Museum of Natural History, Washington, D.C., USA; UT = Non-Vertebrate Paleontology Laboratory, University of Texas, Austin, USA. For geographic occurrence, 0 = Western Europe; 1 = North Africa; 2 = South America; 3 = the WIS.

Materials and methods

Taxonomic sampling.—Twenty-one taxa were selected for phylogenetic analysis, comprising representatives of the “Toxasteridae,” Hemiasteridae, and Micrasteridae, the largest and most morphologically diverse clades within the Spatangoida (Smith and Kroh, 2011). Due to the large number of spatangoid species, not all could be considered in the phylogenetic analysis. Thus, only Cretaceous species, and not Cenozoic and modern forms, were analyzed. All well-preserved Cretaceous WIS species for which at least 75% of the characters used in the phylogenetic analysis could be coded were considered. Additional Cretaceous spatangoid species from outside the WIS were selected using the criteria that: (1) type specimens or quality photographs of these needed to be accessible, (2) species needed to be defined such that at least 75% of characters could be coded, and (3) all major Cretaceous biogeographic regions should be represented. This meant that species that were the subject of more recent systematic treatments were better represented in the analysis. Taxa subjected to phylogenetic analysis and relevant material examined are listed in Table 1. We follow the classification scheme of Smith and Kroh (2011) as it has been widely utilized, although an updated classification has recently been presented by Kroh (2020).

Characters.—Character and character state descriptions were developed via detailed comparative examination of test features. The preexisting literature was also considered to identify characters that had previously been used to diagnose groups and infer relationships among groups. Villier et al. (2004), Smith and Stockley (2005), and Kroh and Smith (2010) were especially helpful in this regard. In total, 32 characters were identified (Appendix), including a combination of qualitative and quantitative characters. Quantitative characters were measured for coding using the program ImageJ; character states are noncontinuous relative to each other, and multistate characters were treated as unordered (Table 2). A variety of means for coding characters exist, and we have followed the precepts and principles presented by Wiley and Lieberman (2011) for identifying hypotheses of homologies and coding them. No explicit models of character weighting were used, but different morphological systems varied in the number of characters they contained. Further, some characters may display more homoplasy than others; thus, different characters may provide varying degrees of phylogenetic resolution.

Phylogenetic methods.—A parsimony analysis was implemented using the matrix in Table 2 via PAUP 4.0 (Swofford, 2003), with Toxaster retusus (Lamarck, 1816) employed as the outgroup; this taxon has been previously treated as an earliest representative of the group by Stockley et al. (2005). (Nexus files utilized are provided as Supplemental Data.) Data were subjected to a branch-and-bound search (Hendy and Penny, 1982) to find the most parsimonious tree(s). A bootstrap analysis was conducted using a heuristic search in PAUP 4.0 to determine support for each node appearing in the most parsimonious tree(s) using 100 bootstrap replicates in a stepwise search that employed five random replicated addition sequences per bootstrap replication with one tree held. Jackknife analysis was also used, again using 100 replicates in a stepwise search that employed five random addition sequences per replication, with the percentage of characters removed value chosen at 10%, which equates to three characters, and one tree held. For both support tests, groups were retained that were compatible with the 50% majority rule consensus tree. In addition, a Bremer branch support analysis (Bremer, 1994) was conducted to calculate node support using the difference between the length of the most parsimonious tree and the largest tree in which a node is maintained.

Table 2.

Character state distributions for taxa used in phylogenetic analysis. ? = missing data.

Finally, a Bayesian analysis was conducted using MrBayes 3.2.7, again using Table 2 and Toxaster retusus as an outgroup, to explore how well supported the various nodes in the parsimony tree were. The original file was exported from the simplified nexus file provided as Supplemental Data. No partitioning of the data was implemented as there were no a priori reasons for differentiating among character types or morphological systems. A standard MKV model was employed without a particular starting tree and performed on the CIPRES cluster. The analysis ran for 10 million generations and sampled every 1,000 generations. The posterior probability values for nodes retrieved in the analysis were then compared with those in the parsimony tree (see Cole, 2019). Note that there has been extensive discussion about whether a Bayesian or a parsimony approach is better justified for the analysis of morphological character data, and consideration of this topic goes substantially beyond the scope of this contribution. We apply these different methods here only for purposes of considering commonality of results achieved across different methods. The interested reader is referred to discussions by Wiley and Lieberman (2011) and Wright (2019) for greater amplification of the various issues involved with applying different phylogenetic methods to morphological data.

Biogeographic analysis.—Biogeographic patterns were inferred by first substituting the area of occurrence of the species for a given species name. The areas utilized were Western Europe, North Africa, South America, and the WIS (Fig. 1), and these represent major areas of endemism for echinoids during the Late Cretaceous (see e.g., Smith, 1984, 1992; Rosen and Smith, 1988; Néraudeau and Mathey, 2000). Biogeographic states were then mapped to ancestral nodes using the modified version of Fitch parsimony optimization for multistate characters, which treats all transformations between different character states as unordered (described by Lieberman, 2000). This method minimizes the assumptions about the way biogeographic change happens between different regions. Transitions between ancestral and descendant nodes or ancestral and descendant terminals that involved range contraction were treated as potential instances of vicariance, and equivalent transitions involving range expansion were treated as potential instances of dispersal following Lieberman (2000).

Repositories and institutional abbreviations.—Type, figured, and other specimens examined in this study are deposited in the following institutions: Division of Invertebrate Paleontology, Biodiversity Institute, University of Kansas (KUMIP); New Mexico Museum of Natural History (NMMNH); United States National Museum of Natural History (USNM); and Non-Vertebrate Paleontology Laboratory, University of Texas, Austin (UT and UT BEG). High-quality images of type specimens were also provided by the Universität Bonn (STIB).

Systematic paleontology

Figure 1.

Paleogeographic reconstruction for Late Cretaceous ∼90 Ma; © 2016 Colorado Plateau Geosystems Inc.

Figure 2.

“Hemiaster” texanus: (1) STIB Texas 156; and “Heteraster” texanus (STIB Texas 157) (2). Scale bars = 1 cm.

Figure 3.

(1) “Diplodetus” americanus (USNM PAL 76285). (2) Proraster humphreysanus (USNM PAL 331). Scale bars = 1 cm.

Results

Phylogenetic analysis.—Two separate parsimony analyses were run. The first included all taxa except for Linthia variabilis Slocom, 1909, which was coded as missing for 25% of the character states. This analysis resulted in a strict consensus of four most parsimonious trees of length 111, consistency index (CI) of 0.41, and retention index (RI) of 0.60 (Fig. 4), with one character parsimony-uninformative. Then, Linthia variabilis was added and a branch-and-bound analysis performed again. (Both nexus files utilized are provided as Supplemental Data.) This resulted in a strict consensus of 83 trees of length 114, a CI of 0.40, and an RI of 0.60 (Fig. 5). Overall, the topology resembles that of Figure 4, but resolution decreased, possibly due to the high amount of missing data for Linthia variabilis. In terms of tree support statistics, four nodes of the cladogram in Figure 4 were supported by bootstrap values greater than 60, with 15 nodes supported by jackknife values greater than 60 (and seven of these values were above 80). Bremer support was found for four nodes in the tree (Fig. 4). The tree in Figure 4, with the characters used in phylogenetic analysis mapped to the tree, is provided as a figure in Supplemental Data. The Bayesian analysis (not shown) of the matrix without Linthia variabilis contained five nodes with probabilities greater than 60, generally aligning with the results of the parsimony analysis but containing poor resolution of later nodes.

Figure 4.

Strict consensus of four most parsimonious trees of length 111 recovered in analysis, excluding relatively incomplete taxon L. variabilis. Numbers above node in bold are bootstrap support values, numbers below node in italics are jackknife support values, and numbers in normal font below node are Bremer support values. *Species occurring outside of the WIS.

The analyses suggest several notable relationships relative to traditional groupings and previous analyses. Regarding traditional groupings, the genera Hemiaster, Heteraster, and Proraster are paraphyletic or even polyphyletic (see also Smith and Bengtson, 1991; Néraudeau, 1994). To make the classification more compatible with the phylogeny, several taxonomic changes were suggested (Fig. 4). For example, Proraster is treated herein as containing Proraster dalli, Proraster humphreysanus, and Proraster bexari (Clark, 1915) (Fig. 4). In addition, some taxa were assigned to paraphyletic genera using the within-quotes convention of Wiley (1979). These are “Heteraster” texanus; “Hemiaster” batnensis; “Hemiaster” texanus; “Hemiaster” estepi; and “Diplodetus” americanus (Fig. 4). At this time, this is preferable to establishing several new monotypic genera as not all relevant species could be considered in the present analysis. The phylogeny indicates that Mecaster fourneli (Aggasiz and Desor, 1847) and Palhemiaster Lambert, 1916 are closely related, with Proraster the sister group of these; a paraphyletic set of species referable to Hemiaster and “Hemiaster” represent several successive lineages sister to these. In addition, the following pairs of genera are closely related: Diplodetus plus Micraster Agassiz, 1836; Washitaster plus Heteraster; and Pliotoxaster Fourteau, 1907 plus Macraster Roemer, 1888. Further, Washitaster plus Heteraster is the sister group to all spatangoids included here (excluding Toxaster Agassiz, 1840).

Figure 5.

Strict consensus of 83 most parsimonious trees of length 114 recovered in analysis, including relatively incomplete L. variabilis. Numbers above node in bold are bootstrap support values, numbers below node in italics are jackknife support values, and numbers in normal font below node are Bremer support values. *Species occurring outside of the WIS.

Biogeographic analysis.—Biogeographic patterns provide evidence for vicariance within the spatangoids (Fig. 6). For example, following an initial range expansion by the group from being distributed in Western Europe and North Africa to being distributed in these areas and the WIS during the Valanginian, there was subsequent vicariance between northern Africa and the WIS plus Western Europe during the Albian comprising Heteraster and Washitaster. There is further vicariance between Western Europe and the WIS within Heteraster and Washitaster. Another vicariance event occured in the Campanian associated with the divergence of the WIS “Diplodetus” americanus and the Western European Micraster and Diplodetus parvistella. The rest of the diversification within the clade represents speciation occurring within the WIS, although at this time it could not be determined whether this involved smaller scale vicariant events within the WIS or actual sympatric differentiation.

There also appears to have been range expansion by spatangoids from the WIS into Europe, North Africa, and South America. For example, at some time between the Albian and the Cenomanian, there was a range expansion by the lineage leading to Hemiaster bufo (Fig. 6). There was also a range expansion from the WIS into South America during the Cenomanian associated with the origin of “Hemiaster” batnensis. Finally, sometime between the Albian and the Turonian, there was a range expansion from the WIS into North Africa and South America associated with the origin of Mecaster fourneli.

Discussion

Phylogenetic analysis.—In general, the phylogenetic results (Fig. 4) match various aspects of previous interpretations, and genera could be readily assigned to families as per the familial level classification presented by Smith and Kroh (2011). For example, Micraster and Diplodetus had been grouped into a monophyletic Micrasteridae due to the shared presence of an ethmophract apical disc and similarly patterned subanal fasciole (Smith and Stockley, 2005), although Diplodetus was found to be paraphyletic herein. Further, various authors (e.g., Smith and Bengtson, 1991; Stockley et al., 2005; Oliveira et al., 2013; Manso and Souza-Lima, 2017) have, either directly or on the basis of their taxonomic assignments, posited a relatively close relationship between Mecaster and Hemiaster. Some aspects of the phylogeny (Fig. 4) are also in agreement with the phylogeny presented by Stockley et al. (2005), although they included a much broader range of taxa. They also found Mecaster and Hemiaster to be closely related, although they did not retrieve a monophyletic Hemiasteridae. They also held Toxaster to be the earliest member of the spatangoids.

Figure 6.

Biogeographic distribution and ages of species mapped to tree shown in Figure 4 along with biogeographic distribution of ancestral nodes, calculated using modified version of Fitch parsimony described by Lieberman (2000), and minimum divergence age of nodes, determined using ghost lineage method described by Norell et al. (1992) and recently used by Guensburg et al. (2019). 0 = Western Europe; 1 = North Africa; 2 = South America; 3 = the WIS.

While the results from the analysis of Villier et al. (2004) cannot be directly compared with the analysis performed herein, as they focused on Early Cretaceous taxa and had fewer WIS representatives, as well as more representatives from outside of the WIS, some similarities emerge. In particular, Villier et al. (2004) also retrieved a clade consisting of Heteraster oblongus, “Heteraster” texanus, and Washitaster riovistae (Adkins, 1920) and found that a clade containing Macraster elegans was sister to two clades, one including Micrasteridae and the other including Hemiasteridae. They did not, however, find Macraster to be closely related to Pliotoxaster comanchei (Clark, 1915) (treated as Palhemiaster in their study).

Biogeographic patterns.—Smith (1984) summarized general patterns of echinoid biogeography and argued that during the Cretaceous, the present-day Mediterranean region (then Tethys) was a biodiversity hotspot for echinoids. The pattern of substantial within-region speciation in the WIS recovered herein suggests that the WIS might also have served as a biodiversity hotspot, at least among spatangoids. The recovery of episodes of range expansion between the WIS and South America matches patterns identified by Néraudeau and Mathey (2000) in their study of the biogeography of the Cretaceous of the present-day South Atlantic region. However, it is important to note that not all species of spatangoids could be considered herein. Changing the species sampled, or adding species, could potentially change the biogeographic character optimizations, which might alter the nature or number of episodes of range expansion, or vicariance; moreover, many species that undoubtedly existed are not preserved in the fossil record. Each of these will affect our ability to reconstruct biogeographic patterns (Turner et al., 2009). Thus, the biogeographic conclusions presented herein focus only on the available data.

Another relevant biogeographic study is that of Rosen and Smith (1988), who used parsimony analysis of endemicity (PAE) on Campanian–Maastrichtian echinoids and found the WIS (in particular, a joint region composed of Texas and Mexico) to be the sister region to North Africa plus Europe. This matches the episodes of vicariance and range expansion between the WIS and Western Europe recovered herein. Further, Smith (1992) analyzed biogeographic patterns in Cenomanian echinoids using distance-based methods, cladistic biogeography, and PAE. The PAE indicated Brazil and Angola were closely related biogeographically and sister to the Texas region. Smith's (1992) cladistic biogeographic analysis indicated that Texas could be the sister area of Western Europe and North Africa. Again, this matches some of the episodes posited herein for vicariance among spatangoids.

It is clear that sea-level changes would have had a major impact on the biogeography and evolution of echinoids (Smith and Rader, 2009). Indeed, the WIS evolved significantly during the time interval during which the taxa considered occurred. Multiple studies have been published on Cretaceous sea-level changes (e.g., Hallam, 1992; Haq, 2014). Prominent sea-level changes are associated with the following transgressive–regressive cycles seen in the North American Cretaceous carbonate platform series: the Coahuilan (Valanginian–mid-Aptian); the Comanchean (late Aptian–early Cenomanian); and the Gulfian (mid-Cenomanian–Maastrichtian) (Scott, 1993; Smith and Rader, 2009). The results from the biogeographic analysis presented herein suggest a possible connection between specific transgressive–regressive cycles and patterns of echinoid vicariance and range expansion. The regressions during the Coahuilan may coincide with the vicariance event that led to the diversification of the WIS “toxasterids,” including Heteraster and Washitaster (Fig. 6); the regression during the Comanchean may have led to the vicariance associated with the divergence at the base of all spatangoids, excluding the “toxasterids” (Fig. 6); another sea-level drop during the Gulfian may coincide with Campanian vicariance in the Micrasteridae (Fig. 6). In addition, sea-level rise associated with the Comanchean could account for the range expansions associated with Mecaster fourneli, Hemiaster bufo, and “Hemiaster” batnensis (Fig. 6). However, a broader range of marine taxa will need to be considered to determine whether these are patterns that consistently emerge in other groups, thereby indicating biogeographic congruence and an overarching earth history control (e.g., Lieberman, 2000). If, instead, they do not, the patterns may be related to specific aspects of spatangoid ecology. At this time, it is not possible to determine the precise paleogeographic features that served as the specific geographic barriers that caused the episodes of vicariance when sea levels fell.