Introduction
Echinoderms show complex patterns of morphological disparity and taxonomic diversity throughout their evolutionary history (Sumrall Reference Sumrall1997; Deline Reference Deline, Zamora and Rábano2015; Deline and Thomka Reference Deline and Thomka2017; Deline et al. Reference Deline, Smith, Zamora, Rahman, Sheffield, Ausich, Kammer and Sumrall2020) and have been suggested to respond to long-term environmental and climatic patterns (Paul Reference Paul1968; Clausen Reference Clausen2004; Dickson Reference Dickson2004; Clausen and Smith Reference Clausen and Smith2005, Reference Clausen and Smith2008; Zamora and Smith Reference Zamora and Smith2008; Rahman and Zamora Reference Rahman and Zamora2009). The early Paleozoic experienced large shifts in geochemical cycles (Haq and Schutter Reference Haq and Schutter2008; Trotter et al. Reference Trotter, Williams, Barnes, Leécuyer and Nicoll2008; Rasmussen et al. Reference Rasmussen, Ullmann, Jakobsen, Lindskog, Hansen, Hansen, Eriksson, Dronov, Frei, Korte, Nielsen and Harper2016; Albanesi et al. Reference Albanesi, Barnes, Trotter, Williams and Bergström2019) that impacted microfossil and marine invertebrate groups (reviewed in Stigall et al. Reference Stigall, Edwards, Freeman and Rasmussen2019), suggesting that blastozoan echinoderms would have been impacted. However, there are limited studies that have analyzed the paleobiogeography of extinct echinoderms to infer how they evolved and dispersed in response to major climatic perturbations.
Of the blastozoan echinoderms, we primarily focus on the Diploporita in this study. Diploporitan echinoderms (Ordovician–Devonian) appeared in the fossil record in Lower Ordovician rocks in the Prague Basin, and quickly gained high biogeographic diversity and morphological disparity (Lefebvre et al. Reference Lefebvre, Sumrall, Shroat-Lewis, Reich, Webster, Hunter, Nardin, Rozhnov, Guensberg and Touzeau2013); all three large groupings within Diploporita (i.e., Asteroblastida, Sphaeronitida, and Glyptosphaeronitida) are found in similarly aged strata within the Prague Basin. Diploporitans were originally placed in the same group by the presence of diplopore (double-pore) respiratory structures, as all blastozoan groups are traditionally classified by their respiratory structures (Sprinkle Reference Sprinkle1973). However, it is suspected that many groups of blastozoans are not monophyletic and that respiratory structures might be convergent (Paul Reference Paul, Paul and Smith1988; Sumrall Reference Sumrall, Harris, Böttger, Walker and Lesser2010; Sheffield and Sumrall Reference Sheffield and Sumrall2017, Reference Sheffield and Sumrall2019a). Sheffield and Sumrall (Reference Sheffield and Sumrall2019a), the first study to place diplopore-bearing echinoderms within a phylogenetic framework, concluded that Diploporita is likely polyphyletic and diplopores have evolved more than once in echinoderms. As Diploporita is not a real evolutionary grouping, we refer to diplopore-bearing blastozoans as diploporans herein.
While biogeographic patterns of blastozoan echinoderms have been explored to a degree (e.g., Lefebvre Reference Lefebvre2007; Nardin and Lefebvre Reference Nardin and Lefebvre2010; Lefebvre et al. Reference Lefebvre, Sumrall, Shroat-Lewis, Reich, Webster, Hunter, Nardin, Rozhnov, Guensberg and Touzeau2013; Zamora et al. Reference Zamora, Lefebvre, J, Àlvaro, Clausen, Elicki, Fatka, Jell, Kouchinsky, Lin, Nardin, Parsley, Rozhnov, Sprinkle, Sumrall, Vizcaino and Smith2013), these studies have treated blastozoan groups as monophyletic, which can limit our understanding of biogeographic and evolutionary patterns. Further, these studies have largely been conducted by analyzing the counts of echinoderm taxa within stratigraphic locations, not through phylogenetic methods. Similarly, many studies of diversity curves, origination, and extinction rates of Paleozoic echinoderms have been conducted at the genus level, without the context of evolution (e.g., Nardin and Lefebvre Reference Nardin and Lefebvre2010). Interpretation of such patterns of evolution within a phylogenetic framework is important, because biogeographic processes are critical for understanding large-scale patterns in biodiversity and paleoecology through time, particularly through the Middle to Late Ordovician (e.g., Miller Reference Miller1997; Wright and Stigall Reference Wright and Stigall2013; Trubovitz and Stigall Reference Trubovitz and Stigall2016; Lamsdell et al. Reference Lamsdell, Congreve, Hopkins, Krug and Patzkowsky2017; Stigall et al. Reference Stigall, Bauer, Lam and Wright2017). The Ordovician is a critical time to investigate in order to understand patterns of evolution and dispersal, as two of the largest diversity booms and busts occur within this period: the great Ordovician biodiversification event (GOBE) and the Late Ordovician mass extinction (LOME). To date, there is no information regarding the speciation types and dispersal patterns of diploporan species, and few regarding echinoderms (a notable exception being Bauer Reference Bauer2020) within a rigorous phylogenetic and statistical framework through the Ordovician Period. As diploporans were likely reevolving respiratory structures through their evolutionary history, it is critical to assess why this pattern occurred. During the Ordovician, organisms like diploporans and other blastozoans were likely responding to changing conditions, such as increasing oxygen levels (e.g., Lenton et al. Reference Lenton, Daines and Mills2018), by evolving these respiratory structures and other morphological features.
In this contribution, we use the R package BioGeoBEARS (Matzke Reference Matzke2013) to infer the biogeographic history of the diploporans through the entire Ordovician Period. Further analysis of the patterns of dispersal and speciation within the group is investigated using biogeographic stochastic mapping (BSM; Dupin et al. Reference Dupin, Matzke, Särkinen, Knapp, Olmstead, Bohs and Smith2017). These analyses are conducted to (1) investigate the origin of the diploporan echinoderms; (2) assess the contribution of vicariance, dispersal, and founder-event speciation (or “jump dispersal”) that led to the early Paleozoic radiation of the diploporans; and (3) reconstruct the source and directionality of dispersal events before, during, and after the GOBE.
Background
The Ordovician Earth System
The Ordovician Period (485.4–443.8 Ma; Ogg et al. Reference Ogg, Ogg and Gradstein2016) was a time of increased tectonic activity, widely dispersed continents, and major climatic shifts. Throughout the entire study interval, volcanic island arcs fringed the major paleocontinents and terranes (e.g., Fitton and Hughes Reference Fitton and Hughes1970; Pedersen et al. Reference Pedersen, Bruton and Furnes1992; Harper et al. Reference Harper, Mac Niocaill and Williams1996; Glen et al. Reference Glen, Walshe, Barron and Watkins1998; Cocks and Torsvik Reference Cocks and Torsvik2011; Samygin Reference Samygin2019; Zhang et al. Reference Zhang, Buckman, Bennett and Nutman2019), especially within the Iapetus Ocean (Rasmussen and Harper Reference Rasmussen and Harper2011; Liljeroth et al. Reference Liljeroth, Harper, Carlisle and Nielsen2017). These tectonic and climatic factors undoubtedly had a large influence on the biogeographic patterns and evolutionary history of blastozoan echinoderms.
The Early Ordovician was characterized by relatively high sea level, with the Iapetus Ocean at its widest. Continents were at their most widely distributed due to increased rates of seafloor spreading and high continental drift rates for Baltica, Siberia, Gondwana, and Avalonia (Scotese and McKerrow Reference Scotese and McKerrow1991; van de Pluijm et al. Reference van de Pluijm, R, Torsvik, Hibbard, van Staal and Cawood1995; Harper et al. Reference Harper, Mac Niocaill and Williams1996; Cocks and Torsvik Reference Cocks and Torsvik2011). Throughout the Ordovician, the Iapetus Ocean between Laurentia and Baltica constricted, leading to the eventual collision of Laurentia, Baltica, and Avalonia in the latest Ordovician to early Silurian (Hutton and Murphy Reference Hutton and Murphy1987; Soper et al. Reference Soper, Strachan, Holdsworth, Gayer and Greiling1992; Cocks and Torsvik Reference Cocks and Torsvik2011; Zagorevski and Van Staal Reference Zagorevski, Van Staal, Brown and Ryan2011).
Several tectonic events characterized the Ordovician Period, with the one most impactful to the Laurentian-centric blastozoans being the Taconian orogeny. Within Laurentia, the initiation of the Taconian orogeny Blountian tectophase began along the southern margin of the paleocontinent during the Middle Ordovician (Ettensohn Reference Ettensohn2010). During this tectophase, the accretion of small island arcs and microcontinents led to the formation of the Sevier foreland basin that stretched from present-day Alabama to Virginia (Shanmugam and Lash Reference Shanmugam and Lash1982; Ettensohn Reference Ettensohn2010). This narrow basin was characterized by major clastic wedges (Ettensohn Reference Ettensohn2004) that did not interrupt tropical carbonate sedimentation across the southern margin of Laurentia (Holland and Patzkowsky Reference Holland and Patzkowsky1997). The Taconic tectophase began during the early Late Ordovician, with the locus of tectonic activity moving east toward present-day New York (Ettensohn et al. Reference Ettensohn, Barnes, Williams, Witman and Roy1991; Ettensohn Reference Ettensohn2010). Nearly synchronous with this event was the creation of the Sebree Trough, a narrow depression that stretched from the southwestern edge of Laurentia toward the midcontinent region, terminating near present-day Ohio (Kolata et al. Reference Kolata, Huff and Bergstro¨m2001). This trough funneled cooler Iapetus waters into the craton, which may have led to mixing of different water masses and a change in intracratonic circulation patterns. Cooler waters generated a density-stratified water column characterized by oxygen-poor and phosphate-rich conditions that lasted until the earliest Katian Age (Wilde Reference Wilde, Barnes and Williams1991; Kolata et al. Reference Kolata, Huff and Bergstro¨m2001; Young et al. Reference Young, Brett and McLaughlin2015; Quinton et al. Reference Quinton, Law, Macleod, Herrmann, Haynes and Leslie2017). With the initiation of cool oceanic waters into Laurentia and expansion of the Sebree Trough (Young et al. Reference Young, Brett and McLaughlin2015), warm carbonate platform deposition was replaced with mixed cool-water carbonate and clastic facies throughout eastern central Laurentia (Keith Reference Keith1989; Ettensohn et al. Reference Ettensohn, Hohman, Kulp and Rast2002). Warm-water carbonate facies resumed during the Late Ordovician Katian Age (Ettensohn et al. Reference Ettensohn, Hohman, Kulp and Rast2002).
Throughout the Ordovician, several climatic shifts occurred, leading to a transition from a greenhouse to an icehouse world. The Early Ordovician was characterized by very warm global average temperatures, with perhaps shorter-term cooling intervals (Trotter et al. Reference Trotter, Williams, Barnes, Leécuyer and Nicoll2008; Albanesi et al. Reference Albanesi, Barnes, Trotter, Williams and Bergström2019) and high global sea level, punctuated by shorter-term sea-level falls (Haq and Schutter Reference Haq and Schutter2008). A shift to a cooler climate began in the late Early Ordovician Floian Age, which continued into the Middle Ordovician and intensified in the Late Ordovician (Trotter et al. Reference Trotter, Williams, Barnes, Leécuyer and Nicoll2008; Finnegan et al. Reference Finnegan, Bergmann, Eiler, Jones, Fike, Eisenman, Hughes, Tripati and Fischer2011; Albanesi et al. Reference Albanesi, Barnes, Trotter, Williams and Bergström2019). Several studies based upon geochemistry, paleoecology, backstripping, and sea-level changes on a stable margin indicate that continental glaciations occurred during the Middle Ordovician, potentially as early as the Darriwilian Age (Vandenbroucke et al. Reference Vandenbroucke, Armstrong, Williams, Zalasiewicz and Sabbe2009; Turner et al. Reference Turner, Armstrong, Wilson and Makhlouf2012; Dabard et al. Reference Dabard, Loi, Paris, Ghienne, Pistis and Vidal2015; Amberg et al. Reference Amberg, Collart, Salenbien, Egger, Munnecke, Nielsen, Monnet, Hammer and Vandenbroucke2016; Pohl et al. Reference Pohl, Donnadieu, Le Hir, Ladant, Dumas, Alvarez-Solas and Vandenbroucke2016a; Rasmussen et al. Reference Rasmussen, Ullmann, Jakobsen, Lindskog, Hansen, Hansen, Eriksson, Dronov, Frei, Korte, Nielsen and Harper2016). Cooling continued into the Late Ordovician with well-developed continental ice sheets over west Gondwana responding to orbital forcing (Beuf et al. Reference Beuf, Biju-Duval, Charpal, Gariel and Bennacef1971; Hambrey Reference Hambrey1985; Crowley and Baum Reference Crowley and Baum1991; Eyles Reference Eyles1993; Crowell Reference Crowell1999; Scotese et al. Reference Scotese, Boucot and McKerrow1999; Sutcliffe et al. Reference Sutcliffe, Dowdeswell, Whittington, Theron and Craig2000; Sinnesael et al. Reference Sinnesael, Desrochers, Rasmussen, Claeys and Vandenbroucke2019).
Modeling experiments under different atmospheric carbon dioxide (CO2) conditions have reconstructed wind-driven surface circulation throughout the Ordovician Period, allowing comparisons of reconstructed dispersal paths to current and gyre flow directions. Five major ocean gyre systems operated within the Ordovician ocean (Poussart et al. Reference Poussart, Weaver and Barnes1999; Herrmann et al. Reference Herrmann, Haupt, Patzkowsky, Seidov and Slingerland2004; Pohl et al. Reference Pohl, Nardin, Vandenbroucke and Donnadieu2016b; Fig. 1). In addition to major gyre systems, several currents flowed along the margins of major continental land masses, such as the Iapetus Current in the Iapetus Ocean. These currents likely caused increased coastal upwelling and thus increased primary productivity along the margins of large paleocontinents, volcanic island arcs, and in the midlatitude regions (Pope and Steffen Reference Pope and Steffen2003; Pohl et al. Reference Pohl, Harper, Donnadieu, Le Hir, Nardin and Servais2018).
Paleogeographic map of the Middle Ordovician Darriwilian Age with major ocean currents (solid lines) and gyre systems (dotted lines). Geographic basins defined in this study are denoted by letters as follows: G, Gondwana; B, Baltica; S, southern Laurentia; P, northern Appalachian Basin; A, southern Appalachian Basin; C, Cincinnati Basin; W, western midcontinent; and R, north of the Transcontinental Arch. Basin letters are the same as those in Fig. 2. Ocean gyres are numbered as follows: 1, north Panthalassic convergence; 2, south Panthalassic convergence; 3, south Paleo-Tethys convergence; 4, north Paleo-Tethys convergence; and 5, Rheic convergence. Ocean currents are abbreviated as follows: NEC, north equatorial current SEC, south equatorial current; IC, Iapetus Current; SLC, southern Laurentia Current; AC, Antarctica Current; SW, southern westerlies; and GC, Gondwanan westerlies. Paleogeographic map is modified from Torsvik and Cocks (Reference Torsvik and Cocks2013), with the approximate locations of currents and gyre systems from Pohl et al. (Reference Pohl, Nardin, Vandenbroucke and Donnadieu2016b).
Biotic Events of the Ordovician
The Ordovician Period not only included substantial changes in geochemical cycles and major climatic shifts, but also a major boom (GOBE) and bust (LOME) in microfossil and marine invertebrate diversity.
The GOBE was a time in Earth's history when there were unparalleled amounts of diversity across multiple taxonomic groups, primarily at the superfamily to species levels (Webby et al. Reference Webby, Paris, Droser and Percival2004; Harper Reference Harper2006; Stigall Reference Stigall2018; Stigall et al. Reference Stigall, Edwards, Freeman and Rasmussen2019). At the family level, diversity increased on average threefold, whereas generic diversity increased fourfold. In a recent review by Stigall et al. (Reference Stigall, Edwards, Freeman and Rasmussen2019), biota belonging to planktic, benthic, and nektonic groups all exhibit peak diversification rates within the Middle Ordovician Darriwilian Age. Thus, the use of the term “GOBE” is now restricted only to the Darriwilian to indicate the rather short period of time when statistically elevated diversification rates occurred (Kröger and Lintulaakso Reference Kröger and Lintulaakso2017).
This expansion of diversification has been linked to extrinsic factors of the Ordovician Earth system, of which many have been proposed, including: oceanic cooling, increased atmospheric oxygen levels, sea-level transgressions, increased weathering and terrestrial delivery of nutrients to the oceans, and spin-up of ocean gyres and currents in response to cooling (Fortey Reference Fortey and Bruton1984; Trotter et al. Reference Trotter, Williams, Barnes, Leécuyer and Nicoll2008; Saltzman et al. Reference Saltzman, Edwards, Adrain and Westrop2015; Pohl et al. Reference Pohl, Donnadieu, Le Hir, Ladant, Dumas, Alvarez-Solas and Vandenbroucke2016a; Rasmussen et al. Reference Rasmussen, Ullmann, Jakobsen, Lindskog, Hansen, Hansen, Eriksson, Dronov, Frei, Korte, Nielsen and Harper2016; Trubovitz and Stigall Reference Trubovitz and Stigall2016; Young et al. Reference Young, Gill, Edwards, Saltzman and Leslie2016; Edwards et al. Reference Edwards, Saltzman, Royer and Fike2017; Lam et al. Reference Lam, Stigall and Matzke2018). Increased atmospheric oxygen levels could have contributed to increased metabolic processes for marine invertebrates and allowed for increased calcite precipitation (Pruss et al. Reference Pruss, Finnegan, Fischer and Knoll2010). Tectonic activity created several volcanic island arcs and microterranes surrounding paleocontinents, which provided the means for species to accomplish long-distance dispersal among areas (Harper et al. Reference Harper, Mac Niocaill and Williams1996; Liljeroth et al. Reference Liljeroth, Harper, Carlisle and Nielsen2017; Lam et al. Reference Lam, Stigall and Matzke2018). All of these factors primed the Ordovician Earth system and its biota for the GOBE, as the evolutionary mechanism that led to an increase in biodiversity was likely produced in part by a series of dispersal and vicariance events, termed “biotic immigration events” (BIMES; Stigall et al. Reference Stigall, Bauer, Lam and Wright2017).
Concurrent with the latest Ordovician was the LOME, which occurred in two major pulses (Brenchley et al. Reference Brenchley, Marshall, Carden, Robertson, Long, Meidla, Hints and Anderson1994; Brenchley Reference Brenchley1995; Brenchley and Marshall Reference Brenchley and Marshall1999), the first occurring at the Katian/Hirnantian boundary and the second at the end of the Hirnantian Age. The first pulse of the extinction is closely linked with the rapid growth of continental ice sheets on Gondwana causing glacioeustatic sea-level and temperature decreases (Finnegan et al. Reference Finnegan, Bergmann, Eiler, Jones, Fike, Eisenman, Hughes, Tripati and Fischer2011). Genera that mainly occupied epicontinental seas were hit hardest by this first extinction pulse as these seaways drained, reducing the amount of available niche space (Finnegan et al. Reference Finnegan, Heim, Peters and Fischer2012, Reference Finnegan, Rasmussen and Harper2016). Modeling studies indicate the combined effects of Ordovician paleogeography along with relatively fast cooling and eustatic sea-level fall are factors that help explain the severity of the first LOME pulse (Saupe et al. Reference Saupe, Qiao, Donnadieu, Farnsworth, Kennedy-Asser, Ladant, Lunt, Pohl, Valdes and Finnegan2020). The second pulse of the LOME corresponds to receding glaciers and warming that may have caused other environmental perturbations, such as warming oceans and anoxic waters transgressing into shelf areas (Brenchley et al. Reference Brenchley, Marshall, Carden, Robertson, Long, Meidla, Hints and Anderson1994; Sheehan Reference Sheehan2001; Melchin et al. Reference Melchin, Mitchell, Holmden and Štorch2013). More recent studies hint that large igneous province (LIP) emplacement in the Late Ordovician, as inferred from mercury-enrichment episodes, could have been a primary driver of the LOME (Jones et al. Reference Jones, Martini, Fike and Kaiho2017).
Methods
Geographic Framework and Species Occurrences
We defined eight biogeographic areas, with six of them within Laurentia (Fig. 1) based on thermal (such as the Sebree Trough) and physical (such as the Transcontinental Arch in western Laurentia) barriers that existed between basins. Laurentian basins include the southern Appalachian Basin; northern Appalachian Basin; Cincinnati Basin, southern Laurentia; north of the Transcontinental Arch; and the western midcontinent (Fig. 1). The Laurentian basins follow marine invertebrate and microfossil provinces (e.g., conodonts, ostracods, corals, brachiopods; Amsden Reference Amsden1974, Reference Amsden1986; Thompson and Satterfield Reference Thompson and Satterfield1975; Elias Reference Elias1983; Barrick Reference Barrick1986; Mohibullah et al. Reference Mohibullah, Williams, Vandenbroucke, Sabbe and Zalasiewicz2012). Laurentia is split into several areas due to increased sampling and species occurrences in North America, which allows for finer-scale paleobiogeographic analyses. Additional basins include the paleocontinents of Gondwana and Baltica. We have grouped Russia, located on the paleo-terrane of Siberia, with Gondwana in this study (Fig. 1), as only one species occurs in this area, and that occurrence is questionable (Supplementary Table 1). The basins in this study follow those defined in Lam et al. (Reference Lam, Stigall and Matzke2018) so a direct comparison among dispersal paths could be made for the blastozoan echinoderms in this study and brachiopods and trilobites studied in Lam et al. (Reference Lam, Stigall and Matzke2018).
The fossil record of the blastozoans, particularly in groups like the diploporans, which do not often preserve in high numbers, is still lacking, but progress over recent years (e.g., Sumrall et al. Reference Sumrall, Heredia, Rodríguez and Mestre2013, Reference Sumrall, Deline, Colmenar, Sheffield, Zamora, Zamora and Rábano2015; Zamora et al. Reference Zamora, Sumrall, Zhu and Lefebvre2016; Sheffield et al. Reference Sheffield, Ausich and Sumrall2017; Sumrall and Zamora Reference Sumrall and Zamora2018) has shed light on the evolutionary history of this group, especially in relatively unexplored areas of the globe. Here, we utilized the phylogeny developed by Sheffield and Sumrall (Reference Sheffield and Sumrall2019a), which tested the hypothesis that Diploporita was not a monophyletic group of echinoderms. To test this hypothesis, other groups of echinoderms (e.g., crinoids, eocrinoids, rhombiferans) were included in the analysis. This tree ranges from the Cambrian to Devonian (Fig. 2); however, because the majority of speciation events occur in the Ordovician, we focus on reconstructing patterns of dispersal and speciation within this period only.
Maximum-likelihood ancestral range estimation of the diploporan and other related blastozoan species within the DIVALIKE + j model. A, Pie charts indicating the probabilities of ancestral ranges at the nodes of the phylogeny. B, Most likely areas as inferred from the pie charts in A. Using these estimates, dispersal (range expansion) and vicariance events were interpreted through the tree. The three time slices used in the biogeographic stochastic mapping (BSM) analysis are denoted beside global chronostratigraphy in the gray boxes. Area letter denotations and colors are the same as in Figs. 1 and 4. Chronostratigraphy and age from the Geologic Time Scale 2016 (Ogg et al. Reference Ogg, Ogg and Gradstein2016). Refer to Supplementary Table 1 for species’ groupings within Echinodermata. Fossil image of diploporan Eumorphocystis multiporata Branson and Peck Reference Branson and Peck1940, SUI 97599, modified from Sheffield and Sumrall (Reference Sheffield and Sumrall2019b).
For species utilized in this study, their stratigraphic and geographic occurrence data were compiled from published literature sources, the Paleobiology Database, and FossilID (https://fossiilid.info; Supplementary Table 1). Age-level occurrence data from published literature sources were then updated to the Geologic Time Scale 2016 (Ogg et al. Reference Ogg, Ogg and Gradstein2016) using graptolite and conodont zones to gain an estimate of species ranges for use in time-calibrating the blastozoan phylogeny.
Time-Calibrated Phylogenetic Hypothesis
Sheffield and Sumrall (Reference Sheffield and Sumrall2019a) constructed a morphological character matrix of 61 characters (of which 41 were parsimony informative). Twenty-eight taxa across Diploporita and other stemmed echinoderms were included in the analysis; this study primarily used taxa that the authors were able to study in person. Parsimony-uninformative characters were removed from the analysis. A heuristic search of most parsimonious trees was run using a tree–bisection–reconnection branch-swapping algorithm with a reconnection limit of eight. The matrix was polarized using Gogia spiralis Robison, Reference Robison1965 as the out-group. The analysis recovered 18 optimal trees. The analysis presented here uses the strict consensus of the most optimal trees recovered in Sheffield and Sumrall (Reference Sheffield and Sumrall2019a).
Biogeographic models require a fully bifurcating, time-calibrated phylogenetic tree. The polytomies present in Sheffield and Sumrall (Reference Sheffield and Sumrall2019a) were resolved in the Bayesian analysis presented herein. Overall, the topologies of both trees are largely very similar, although a few relationships changed. The monophyletic grouping uncovered in Sheffield and Sumrall (Reference Sheffield and Sumrall2019a), Sphaeronitida (which also includes the monophyletic grouping Holocystitidae), is present in the fully bifurcating tree, though relationships within sphaeronitids have changed (e.g., Holocystites salmoensis is the sister taxon to other holocystitids, though this group is still monophyletic as well). The other two traditional hierarchies within Diploporita, Glyptosphaeritida and Asteroblastida, share similarities in both phylogenetic reconstructions. Glyptosphaeritida is polyphyletic in both reconstructions and, notably, the taxon Eumorphocystis multiporata is sister taxon to Hybocrinu nitidus, a crinoid, something uncovered in another recent phylogenetic reconstruction (Sheffield and Sumrall Reference Sheffield and Sumrall2019b). Asteroblastid (Asteroblastus stellatus) groups away from other diploporans. In Sheffield and Sumrall (Reference Sheffield and Sumrall2019a), this species is the sister taxon to edrioblastoid Eurekablastus rozhnovi; in this analysis, it is sister to E. multiporata and Hybocrinus, while E. rozhnovi is the sister taxon. Rhombiferan echinoderms are uncovered as polyphyletic in both analyses. Eocrinoid Rhopalocystis destombesi was uncovered in both analyses as the taxon most distantly related to the other taxa in this analysis and groups most closely with G. spiralis, the out-group.
As our goal was to obtain a fully bifurcating phylogeny, with more realistic divergence dates than obtained from minimum ages alone, we reanalyzed the character matrix using Bayesian tip-dating methods (Matzke and Wright Reference Matzke and Wright2016) in Beast2 (Bouckaert et al. Reference Bouckaert, Heled, Kühnert, Vaughan, Wu, Xie, Suchard, Rambaut and Drummond2014). BEASTmasteR (Matzke Reference Matzke2019) was then used to convert the data matrix into the XML format that is required for Beast2.
As fossil specimen temporal ranges are not known exactly, and the characters of each species occur across the time range of that species, the date of each tip species was allowed to vary according to a uniform prior stretching from the bottom of the time bin of the first occurrence to the top of the time bin of the last occurrence. The exception was Tristomiocystis globosus, the most recently occurring taxon in the dataset, which was given a fixed age of 382.7 Ma in order to provide a fixed time point for the other time priors. To account for the fact that invariant characters are not collected in parsimony analyses, the Mk v model (Lewis Reference Lewis2001) was implemented in the Beast2 XML via BEASTmasteR (Matzke and Irmis Reference Matzke and Irmis2018). Among-site rate variation was modeled with a discretized gamma distribution with four categories. A relaxed clock model, with uncorrelated lognormally distributed rates, was used to model variation in the rate of morphological evolution among branches. The clock rate had a flat uniform (0.00001, 2) prior, and the clock standard deviation had a flat uniform (0, 10) prior; this combination covers all reasonable values for the average rate and its variation among branches. A fossilized birth–death prior (using Beast2's BDSKY model) was used for the tree prior, as it allows sampling of fossils through time. See Lam et al. (Reference Lam, Stigall and Matzke2018) for an expanded discussion of the BEASTmasteR setup. The Bayesian analysis was run for 1 billion generations, sampling every 500,000 generations, producing 2000 samples. After 10% burn-in, all parameters had estimated sample-size values above 419, larger than the recommended minimum of 200. TreeAnnotator was used to summarize the post-burn-in posterior distribution of dated phylogenies by generating a maximum clade credibility tree, used for downstream analyses.
Ancestral Range Estimation
Probabilistic inference of ancestral geographic ranges within the phylogeny was conducted in the R package BioGeoBEARS (Matzke Reference Matzke2013). BioGeoBEARS implements, in a common likelihood framework, the key assumptions from three programs popular in historical biogeography: the dispersal–extinction–cladogenesis (DEC) model of Lagrange (Ree and Smith Reference Ree and Smith2008); dispersal–vicariance–analysis (DIVA; Ronquist Reference Ronquist1997); and BayArea (Landis et al. Reference Landis, Matzke, Moore and Huelsenbeck2013). In BioGeoBEARS, DIVA and BayArea are converted into a likelihood framework directly comparable to DEC, and are thus referred to as DIVALIKE and BAYAREALIKE within the package. The models allow for geographic ranges to evolve along branches in the phylogenetic hypothesis using range expansion (dispersal) and range loss (“extinction,” more accurately termed “local extirpation”) processes, modeled with the rate parameters d and e, respectively. Each model makes different assumptions about which cladogenetic processes are allowed (Matzke Reference Matzke2013), and each model can be modified with an additional + j parameter to model the relative probability of founder-event speciation (or jump dispersal) during cladogenesis (Matzke Reference Matzke2014). This mode of speciation has been found to be important for many extant clades (Matzke Reference Matzke2014; Dupin et al. Reference Dupin, Matzke, Särkinen, Knapp, Olmstead, Bohs and Smith2017), as well as Ordovician brachiopods and trilobites (Lam et al. Reference Lam, Stigall and Matzke2018; Censullo and Stigall Reference Censullo and Stigall2019; Congreve et al. Reference Congreve, Krug and Patzkowsky2019).
Although a recent critique argues that DEC and DEC + j cannot be statistically compared (Ree and Sanmartín Reference Ree and Sanmartín2018), that critique can be rejected on several grounds (Klaus and Matzke Reference Klaus and Matzke2019). For example, the DEC/DEC + j statistical comparison was validated with simulation-inference tests in Matzke (Reference Matzke2014), but these tests were not rebutted, or even acknowledged, by Ree and Sanmartín (Reference Ree and Sanmartín2018). Instead, the critique relies on peculiar behavior found with two tiny datasets (of 2 or 4 taxa), but maximum-likelihood inferences are not expected to be reliable for such small datasets where the number of data is close to the number of inferred parameters (2–3). Similarly, the claimed pathological inferences (jump dispersal at all nodes) are not observed on typical empirical datasets (McDonald-Spicer et al. Reference McDonald-Spicer, Knerr, Encinas-Viso and Schmidt-Lebuhn2019).
The BioGeoBEARS analysis was conducted without constraining the direction or timing of speciation events. The maximum number of areas that species could occupy (‘max_range_size’) was set to 2, which is the maximum number of areas any species occupied in the phylogeny (Supplementary File 2). This limits the size of the state space, speeding up the computational analysis. Relative model fit was measured with sample-size corrected Akaike information criterion (AICc; Burnham and Anderson Reference Burnham and Anderson2002). Following O'Meara et al. (Reference O'Meara, Ané, Sanderson and Wainwright2006), AICc values were calculated using the total number of taxa in the phylogeny as the sample size.
Estimation of Type and Number of Dispersal Events
The number and percentage of dispersal events (both anagenic range expansion and jump dispersal) was estimated using BSM implemented within BioGeoBEARS (Dupin et al. Reference Dupin, Matzke, Särkinen, Knapp, Olmstead, Bohs and Smith2017). Stochastic mapping techniques were first developed for discrete characters like morphology and DNA (Nielsen Reference Nielsen2002; Huelsenbeck et al. Reference Huelsenbeck, Nielsen and Bollback2003; Bollback Reference Bollback2006). Instead of calculating the probability of each character state at each node (traditional ancestral state probabilities; Felsenstein Reference Felsenstein2003), stochastic mapping simulates an exact history of character states across an entire tree (often represented by “painting” or “mapping” the character state history onto the phylogeny). Every “stochastic map” is a probabilistic sample of the possible histories of the character, conditional on the phylogeny, the model parameters, and the observed tip data. BSM implements the same method, using biogeographic models. The advantage of BSM over standard ancestral range probabilities is that exact event histories are sampled, allowing an estimation of the number and directionality of different anagenetic and cladogenetic events, along with uncertainties in these counts. The oft-used approach of “eyeballing” the number and directionality of dispersal events from a plot of ancestral range probabilities, while correlating with the true counts, will often miss rarer or uncertain events, resulting in underestimation of the number of events. In addition, eyeball counts produce no estimate of uncertainty. The accuracy of BSM can easily be checked, as averaging many BSMs will asymptotically approach the analytically calculated ancestral range probabilities (Dupin et al. Reference Dupin, Matzke, Särkinen, Knapp, Olmstead, Bohs and Smith2017).
Although the entire phylogeny was subject to ancestral range estimation in BioGeoBEARS, we only focused the BSM analysis on the Ordovician sections of the tree. We time-sliced the BSM analysis to determine the amount of dispersal (range expansion) and jump dispersal within three time slices and to determine whether dispersal paths changed during the Ordovician Period, especially within the GOBE interval (Supplementary File 2). The first time slice encompasses the Early to Middle Ordovician Tremadocian to Dapinigian ages (485.4–467.3 Ma). The second time slice includes just the Middle Ordovician Darriwilian Age (467.3–458.4 Ma), the time of the GOBE. The third time slice covers the Late Ordovician Sandbian to Hirnantian ages (458.4–443.8 Ma), including the LOME (Fig. 2).
Results
Ancestral Range Estimation
In all cases, the addition of the + j parameter (founder-event speciation) improved model fit (Table 1), suggesting that anagenetic dispersal (range expansion) and sympatry (within-area speciation) alone were not sufficient to fully explain the biogeographic patterns within the clade. An improved model fit with the addition of the + j parameter was also observed for Ordovician brachiopod and trilobite species and genera (Lam et al. Reference Lam, Stigall and Matzke2018; Censullo and Stigall Reference Censullo and Stigall2019; Congreve et al. Reference Congreve, Krug and Patzkowsky2019).
Comparison of model fits for the Blastozoan BioGeoBEARS analysis, with the best-fit model DIVALIKE + j bolded. lnL, log-likelihood values; AIC, Akaike information criterion; AICc, corrected AIC; ΔAIC, difference between the best-fit model and other models; ωi, Akaike weight, the relative likelihood of each model; d, dispersal rate (range expansion) along each branch within the phylogeny; e, extinction rate (range contraction) along each branch within the phylogeny; j, relative weight of jump dispersal in each model.
The best-fit model for the blastozoan tree presented here was DIVALIKE + j, with the DEC + j model having a slightly weaker fit. The results from these analyses are very similar, with only slight differences (e.g., in the DEC + j analysis, the southern Appalachian Basin has a higher probability of being a source region for the ancestor of E. multiporata and H. nitidus; Supplementary File 2, pie charts of DEC + j and DIVALIKE + j). However, the most probable ranges from BioGeoBEARS for the blastozoans are identical (Supplementary File 2, most probably ancestral states of DEC + j and DIVALIKE + j). Therefore, we have focused our discussion on the results from the DIVALIKE + j model.
Ancestral range estimation using DIVALIKE + j indicates that there are potentially several source regions for the species in this analysis, as there was not a very high probability of any one area as ancestral (Fig. 2A). Deep nodes in the phylogeny do indicate that Baltica, Gondwana, and western Laurentia (north of the Transcontinental Arch) could be source regions for the earliest species. The geographic ranges of the ancestors of diploporan blastozoans (Fig. 2A) are highly uncertain, so these inferences should be viewed as hypotheses that can be tested by collecting specimens and geographic range data near in time to the nodes in question; species with smaller ranges would tend to reduce the estimate range size of nearby nodes, as well as the uncertainty.
For comparison with previous studies, biogeographic events were also inferred from the most probable ranges (Fig. 2B) under the DIVALIKE + j model. Dispersal was inferred by a descendant most-probable range occupying a new area than that of its ancestor, and vicariance interpreted when the descendant species occupied subsets of the ancestor. From the tree with most likely areas (Fig. 2B), five dispersal events were recovered in the first time slice, four in the second time slice, and six in the third time slice. Throughout the tree, the majority of speciation is dominated by dispersal, with only one vicariance event occurring within the Cambrian. Interestingly, this clade does not display patterns of alternating vicariance and dispersal events (BIMES) that have been recorded in Ordovician brachiopod and trilobite clades (Stigall et al. Reference Stigall, Bauer, Lam and Wright2017; Lam et al. Reference Lam, Stigall and Matzke2018; Censullo and Stigall Reference Censullo and Stigall2019; Congreve et al. Reference Congreve, Krug and Patzkowsky2019).
Biogeographic Dispersal Type, Numbers, and Directionality
Results from the Ordovician BSM analysis on the DIVALIKE + j models indicate that range expansion (dispersal), founder-event speciation (jump dispersal), and within-area speciation were the dominant biogeographic events for these blastozoans throughout the Ordovician Period (Table 2). As the phylogeny contains only 26 species within the Ordovician Period, the total number of events is relatively low. Across all time slices, vicariance estimates were very low, with “extinction” (meaning anagenetic range loss) estimated at 0. Within-area speciation events are expected to be higher in this study due to the large geographic areas defined outside Laurentia (e.g., Baltica and Gondwana). Among the dispersal types, founder events were more than twice as common as range expansions throughout the Ordovician (Table 2). For this reason, we focus the following discussion on dispersal patterns of the species in this analysis (both founder events and range expansions) through time.
Biogeographic stochastic mapping results for the diploporan and other blastozoan species across the three Ordovician times. Mean values, standard deviations (SDs), and the percent of each speciation mode are shown. There were no extinction or range contractions estimated for any of the times. Time 1 is the Early to early Middle Ordovician (Tremadocian to Dapingian ages, 485.4–467.3 Ma); time 2 encompasses the Middle Ordovician Darriwilian Age, or GOBE interval (467.3–458.4 Ma); and time 3 includes the Late Ordovician Sandbian to Hirnantian ages (458.4–443.8 Ma).
Throughout the Early to early Middle Ordovician, or Tremadocian to Dapingian ages, founder events comprised nearly half of all speciation events, and range expansions made up 16% of the estimated dispersal events. Vicariance is lowest within this interval of time, and within-area speciation (sympatry) comprises 32% of total estimated events. Focusing on dispersal events, it is clear that Baltica was a major source region for dispersals, as this area comprises more than a third of the estimated dispersal events (Fig. 3). The western midcontinent region within Laurentia was the second-highest source region for dispersals, with the majority of movement occurring into the southern Appalachian Basin. Gondwana was the third-highest contributor to estimated dispersal events during the Early to early Middle Ordovician, with an estimated 16.8% of events occurring from Gondwana. The majority of Gondwanan dispersals occurred into the western midcontinent and north of the Transcontinental Arch within Laurentia. Areas that acted as major sinks, or areas where species dispersed into, include the western midcontinent, Gondwana, and southern Laurentia, respectively. Only 10% of estimated dispersal events occur into Baltica, indicating this region was mainly a place out of which species dispersed. In this time slice, the BSM results indicate that there was increased faunal exchange between the southern Appalachian Basin and western midcontinent, and between Baltica and southern Laurentia.
Number of dispersal (range expansion and jump dispersal) events for the three time slices (Fig. 2) from the DIVALIKE + j model. Counts of events were averaged across the 100 biogeographic stochastic mappings (BSMs) with standard deviations in parentheses. Rows represent source areas, columns represent destinations (or sinks). Darker shades indicate a higher frequency of dispersal events. The sum and percent of events in each row and column are given on the margins. Areas are abbreviated as follows: S. Laur., southern Laurentia; W. Mid., western midcontinent; Gond., Gondwana; Cincin., Cincinnati Basin; N. App, northern Appalachian Basin; S. App, southern Appalachian Basin; N. TCA, north of the Transcontinental Arch.
Speciation type changed very little during the Middle Ordovician Darriwilian Age, which encompasses the GOBE (Table 2), but dispersal source and sink regions changed compared with the first time slice. Founder-event speciation was still the dominant speciation mode, with a slight increase in the number of estimated range expansions. Vicariance increased slightly during this time, but was still generally very low. Within-area speciation dropped, indicating sympatric speciation became less important during this time due to species evolving in different regions. During this time, the source regions for dispersals switched to Laurentian basins, with several events still occurring from Baltica (over 45%; Fig. 3). Namely, the Cincinnati Basin and western midcontinent were the largest source areas for dispersal events. Gondwana was still important, but its importance as a source for species dispersal was reduced during the Middle Ordovician (11.3% in time 2 compared with 16.8% in time 1). Instead, it acted as a sink for species (Fig. 3), with 22.2% of estimated dispersals occurring into Gondwana. Basins within Laurentia, namely the Cincinnati Basin and the western midcontinent, were also major regions for species to disperse into. During this time, there was a strong exchange of taxa between the southern Appalachian Basin and western midcontinent, and the Cincinnati Basin and Baltica (Fig. 3).
The third time slice encompasses the Late Ordovician Sandbian to Hirnantian ages, which include the LOME. In this time slice, founder-event speciation events increased to their highest estimates of the three time slices, with the number of range expansions remaining unchanged from the Middle Ordovician GOBE time slice (Table 2). Vicariance and within-area speciation both decreased slightly compared with the Middle Ordovician. Areas that acted as sources for dispersal events changed little from the second time slice, but areas into which species dispersed differed greatly. Baltica remained one of the dominant areas from which the most estimated dispersals occurred, with most dispersals occurring out of the western midcontinent of Laurentia. Together, dispersal events from these two areas comprise 77% of all estimated dispersal events in the third time slice. The three areas that were most important as sink regions include the Cincinnati Basin, Gondwana, and southern Laurentia, respectively. Within the Late Ordovician, estimated dispersal events between Baltica and the Cincinnati Basin increased from the previous two time slices, indicating an increasing number of events between these regions through the entire Ordovician (Fig. 3). There are also increased dispersal events between the western midcontinent and southern Laurentia and between the western midcontinent and the southern Appalachian Basin during this time compared with earlier time slices.
Discussion
Our analyses indicate that there are still several questions remaining as to where these groups originated. The most probable ranges estimated from the BioGeoBEARS analysis indicate the earliest likely ancestors of the species represented, based on their respiratory structures (G. spiralis; Fig. 2), originated in Gondwana and western North American (north of the Transcontinental Arch), with the earliest diploporan blastozoans occurring within the Ordovician. Ordovician species likely also originated within Gondwana, with deep nodes in the Early Ordovician, indicating Baltica may also be a place of major evolutionary radiations for this group. Additional occurrence data across Blastozoa will allow for increased division of areas to more finely determine the origins for this clade. The first documented fossil occurrences of all of the major groups within the traditional group of Diploporita (i.e., Asteroblastida, Sphaeronitida, and Glyptosphaeritida) appear in the same age strata in the Ordovician Prague Basin of Gondwana. Because the morphology of each of these groupings is disparate and diploporans are polyphyletic, we hypothesize that it is possible that diploporans have definitive origins earlier than the Ordovician. We suspect that these diplopore respiratory structures likely evolved in the late Cambrian; this analysis may give an indication as to the area in which these organisms might be found, as data indicate that we see strong dispersal patterns from Gondwana.
An understanding of the limits of echinoderm dispersal (e.g., thermal barriers) would add immensely to this study. However, there is limited knowledge about fossil echinoderm larvae, especially within diplopore-bearing echinoderms, for which juvenile specimens are rarely found (Sheffield et al. Reference Sheffield, Ausich and Sumrall2017). Presumably, Paleozoic echinoderms dispersed during larval stages, as modern echinoderms do today, but there is little information available about the larval stages of Paleozoic echinoderms. Some studies using more modern echinoderms, such as ophiuroids, indicate that long-distance dispersals are possible (Richards et al. Reference Richards, DeBiasse and Shivii2015), while other studies question whether long-distance dispersal in Paleozoic echinoderms was possible (e.g., Rozhnov Reference Rozhnov2013); however, echinoderms have extremely varied larval morphologies, and we cannot currently estimate whether or not diplopore-bearing groups would have dispersed similarly.
Abiotic Drivers of Dispersal and Speciation
It is apparent from the reconstructed ancestral areas and dispersal paths that the evolutionary history of the blastozoans was influenced by major ocean currents and gyre systems, and these paths changed through the Ordovician time slices defined in this study. Recovered dispersal paths from the most likely areas reconstructed from BioGeoBEARS (Fig. 2B) and the BSM analysis (Fig. 3) were compared with the surface ocean circulation model of Pohl et al. (Reference Pohl, Nardin, Vandenbroucke and Donnadieu2016b) for the Middle Ordovician (460 Ma; Fig. 1) under 8 times preindustrial atmospheric CO2 levels (PAL). From the Pohl et al. (Reference Pohl, Nardin, Vandenbroucke and Donnadieu2016b) analysis, there is little change in ocean circulation throughout the Ordovician; thus we use this reconstruction for the entire period, although we acknowledge that changing tectonic arrangements of the continents and terranes likely influenced changes in upwelling zones and circulation patterns.
Interestingly, the majority of dispersal events occurred during the Early to early Middle Ordovician time slice, when the continents were most widely distributed (Scotese and McKerrow Reference Scotese and McKerrow1991; Fig. 4). Namely, there were several dispersal events between Gondwana and Laurentia, as well as between Baltica and Laurentia (Fig. 4). This finding indicates that distance among terranes and continents was not an important factor in the dispersal and evolutionary history of the diploporans. Because our best-fit model incorporated jump dispersal (Table 1), species were likely dispersing among areas defined in this study through the use of microterranes and volcanic island arcs that were surrounding major continents (Cocks and Torsvik Reference Cocks and Torsvik2011).
Dispersal directions for the three time slices inferred from the most likely areas occupied by ancestors in the phylogeny (bold lines, Fig. 2B) and from dispersal events recovered within the biogeographic stochastic mapping (BSM) analysis (thin lines, Fig. 3). The majority of long-distance dispersal events occur within the first and third time slices. Dispersal events were mainly constrained between Baltica and Laurentia and between Laurentian basins in the second time slice. Basin and area shading are the same as in Figs. 1 and 2.
The Early Ordovician, which was characterized by low atmospheric oxygen levels, generally warm seas (Fig. 5), and a lower equator to pole temperature gradient due to no or few glaciers compared with the Middle and Late Ordovician (Trotter et al. Reference Trotter, Williams, Barnes, Leécuyer and Nicoll2008; Pohl et al. Reference Pohl, Donnadieu, Le Hir, Ladant, Dumas, Alvarez-Solas and Vandenbroucke2016a; Edwards et al. Reference Edwards, Saltzman, Royer and Fike2017; Quinton et al. Reference Quinton, Speir, Miller, Ethington and MacLeod2018), was likely a time of few thermal ocean barriers for blastozoan dispersal. Thus, with warm seas and plenty of island arcs surrounding terranes and continents, dispersal occurred relatively unimpeded during this time. Dispersal was likely mediated among areas by the almost zonal flow of surface currents from Gondwana to Baltica, and subsequently to Laurentia via the Iapetus Current (Fig. 1; see also Pohl et al. Reference Pohl, Nardin, Vandenbroucke and Donnadieu2016b: Fig. 5b). Significant dispersal events between Laurentia and Baltica (Figs. 3, 4) were likely influenced by the Iapetus Current and the anticyclonic surface currents between the two continents in the Iapetus Ocean that are prominent under models with 16 and 8 PAL (Pohl et al. Reference Pohl, Nardin, Vandenbroucke and Donnadieu2016b).
Geochemical and evolutionary data for the Ordovician Period, with the great Ordovician biodiversification event (GOBE) indicated by the gray horizontal bar. Blastozoan speciation events are scattered evenly through the Darriwilian to Katian ages, making interpretations of abiotic drivers of evolution in this dataset difficult. Sea-level data are from Haq and Schutter (Reference Haq and Schutter2008); generalized δ13C carbonate curve from Bergström et al. (Reference Bergström, Chen, Gutierrez-Marco and Dronov2009), plotted with mercury (Hg) over total organic carbon (TOC) from Monitor Range, Nevada (gray line) from Jones et al. (Reference Jones, Martini, Fike and Kaiho2017); δ18O curve from Albanesi et al. (Reference Albanesi, Barnes, Trotter, Williams and Bergström2019); red atmospheric CO2 curve from the GEOCARBSULF model of Krause et al. (Reference Krause, Mills, Zhang, Planavsky, Lenton and Poulton2018); black curve of oxygen-dependent atmospheric CO2 model from Edwards et al. (Reference Edwards, Saltzman, Royer and Fike2017); and blue curve of atmospheric CO2 from the COPSE biogeochemical model of Lenton et al. (Reference Lenton, Daines and Mills2018). Brachiopod and trilobite speciation events from Lam et al. (Reference Lam, Stigall and Matzke2018). Abbreviations for carbon isotope excursions are as follows: ICE, isotope excursion; MDICE, mid-Darriwilian excursion; GICE, Guttenberg excursion; and HICE, Hirnantian excursion.
Among Laurentian basins, the highest number of dispersal events as recovered from the BSM analysis (Fig. 3) occur within the Early to early Middle Ordovician time slice. Our finding corroborates other statistical studies of blastozoan speciation that indicate endemism among genera was depressed during the Early Ordovician (Nardin and Lefebvre Reference Nardin and Lefebvre2010), likely due to increased dispersal events. To date there are no models that have captured surface water flow within the Laurentian epicontinental seas, but based on reconstructed wind vectors, surface currents likely flowed from the northeast to southwest (Ettensohn Reference Ettensohn2010), with major hurricanes sweeping across the paleocontinent from 10° to 30° north and south of the equator (Jin et al. Reference Jin, Harper, Cocks, McCausland, Rasmussen and Sheehan2013). As concluded in other studies that have focused on intracratonic Laurentian dispersal events of marine invertebrates (Wright and Stigall Reference Wright and Stigall2013; Bauer and Stigall Reference Bauer and Stigall2014; Lam and Stigall Reference Lam and Stigall2015; Lam et al. Reference Lam, Stigall and Matzke2018), we assume these mechanisms were the primary drivers of dispersal among basins. The lack of tectonic and thermal barriers that later resulted from the Taconian orogeny (e.g., the Sebree Trough) in the Early Ordovician, combined with strong wind and storm activity, likely led to increased blastozoan dispersal events among areas and decreased endemism during this time.
The second time slice in our analysis corresponds to the Darriwilian Age and the GOBE interval (Fig. 5), a time when oceans cooled with the growth of continental land glaciers over Gondwana (Quinton et al. Reference Quinton, Percival, Zhen and Macleod2015; Pohl et al. Reference Pohl, Donnadieu, Le Hir, Ladant, Dumas, Alvarez-Solas and Vandenbroucke2016a; Edwards et al. Reference Edwards, Jones and Fike2018; Albanesi et al. Reference Albanesi, Barnes, Trotter, Williams and Bergström2019). During this time, ocean circulation may have vigorously increased as the temperature gradients from the equator to poles increased and thus aided in dispersal dynamics of benthic marine invertebrates during their larval phases (Rasmussen et al. Reference Rasmussen, Ullmann, Jakobsen, Lindskog, Hansen, Hansen, Eriksson, Dronov, Frei, Korte, Nielsen and Harper2016; Lam et al. Reference Lam, Stigall and Matzke2018; Stigall et al. Reference Stigall, Edwards, Freeman and Rasmussen2019), in part leading to the continued low endemicity of blastozoans through the Darriwilian Age (Nardin and Lefebvre Reference Nardin and Lefebvre2010). This interval of our study is characterized by a decreased number of dispersal events as interpreted from the most likely areas reconstructed within BioGeoBEARS (Fig. 2) and the BSM analysis (Fig. 3), but only by 0.79 compared with the first time slice and 0.36 compared with the third time slice. The reduced number of dispersal events may be due to the second time slice being of a much shorter duration, as compared with the first and third time slices, but the data also indicate that the fewer dispersal paths recovered during the Darriwilian Age are of higher counts (Fig. 3).
Compared with the Early to early Middle Ordovician time slice, there are several patterns that differ during the Darriwilian Age. Within the Darriwilian, it is clear that there were several dispersal events between Laurentia and Baltica, namely to and from the Cincinnati Basin (Figs. 3, 4). Fewer dispersal events occur among Laurentia, Baltica, and Gondwana, and even fewer among Laurentian basins (Fig. 4). The increased number of dispersal events between Baltica and the Cincinnati Basin in the Darriwilian Age compared with the first time slice (Fig. 3) is likely due to the movement of Baltica closer to Laurentia (Cocks and Torsvik Reference Cocks and Torsvik2011). With increased cooling (Fig. 5) and intensification of currents and upwelling due to reduced temperatures and increased latitudinal temperature gradients (Pohl et al. Reference Pohl, Nardin, Vandenbroucke and Donnadieu2016b; Rasmussen et al. Reference Rasmussen, Ullmann, Jakobsen, Lindskog, Hansen, Hansen, Eriksson, Dronov, Frei, Korte, Nielsen and Harper2016; Lam et al. Reference Lam, Stigall and Matzke2018), the Iapetus Current (Fig. 1) was likely a major factor facilitating trans-Iapetus dispersal. On the southern margin of Laurentia, the Taconian orogeny Blountian tectophase was in full swing during this time (Ettensohn Reference Ettensohn2010). Our data indicate that limited dispersal among Laurentian basins may have been caused by increased tectonic activity, which would have led to increased cool-water influx into the craton through the Sebree Trough (Young et al. Reference Young, Brett and McLaughlin2015) and the addition of physical barriers blocking species dispersal.
The triggers of the GOBE have been discussed in great detail in previous works, including increased atmospheric oxygen levels, a cooling climate and ocean, and cessation of anoxic events, coinciding in part with increased diversity among microfossil and invertebrate clades (reviewed in Stigall et al. Reference Stigall, Edwards, Freeman and Rasmussen2019). From our new analyses on this blastozoan echinoderm phylogeny, it is clear that these Earth system changes may not have affected diploporans and other blastozoans as they did other shelly marine invertebrates, or perhaps global cooling and the growth of Gondwanan glaciers (Pohl et al. Reference Pohl, Donnadieu, Le Hir, Ladant, Dumas, Alvarez-Solas and Vandenbroucke2016a) negatively affected them. Our data indicate slightly reduced and limited dispersal events and paths during the GOBE interval, compared with the Early Ordovician, along with low endemism (Nardin and Lefebvre Reference Nardin and Lefebvre2010). Diversity curves for the blastozoans at the genus level in other studies do suggest an increase in diversity during the Darriwilian (Nardin and Lefebvre Reference Nardin and Lefebvre2010). Because of increased diversity and low endemism, we would expect to see increased dispersal events among several basins during this time. We suggest additional analyses of blastozoan diversity and endemism at the species level be conducted before trying to further explore patterns and drivers of blastozoan speciation across the GOBE interval. Other echinoderm studies indicate that members of this large and diverse group did obtain increased richness and diversity during the GOBE (e.g., Wright and Toom Reference Wright and Toom2017; Cole Reference Cole2019), hinting that the blastozoans may have followed a similar trend.
During the Late Ordovician time slice, there is again an increased number of dispersal events recovered from the BioGeoBEARS and BSM analyses. However, dispersals from the Late Ordovician BSM analysis are still less than those recovered for the Early to early Middle Ordovician time slice (Fig. 3). During this interval, species dispersals among Laurentia, Gondwana, and Baltica are prevalent, with several dispersal events apparent between Gondwana and Laurentia (Fig. 4). The highest number of dispersals between Baltica and the Cincinnati Basin occur during this time (Fig. 3), indicating increased connectivity between these regions through the entire Ordovician Period as these continents moved closer together. Interestingly, the majority of dispersals into Laurentian basins occur into areas that are located on the southern margin. This may be due to the deposition of warm-water facies that resumed in the Katian Age (Ettensohn et al. Reference Ettensohn, Hohman, Kulp and Rast2002), providing more niche space within these regions for species.
Several Earth system changes occurred during the Late Ordovician, namely the possible emplacement of a LIP (Jones et al. Reference Jones, Martini, Fike and Kaiho2017) in the late Katian to Hirnantian ages (Fig. 5), increased global cooling and the growth of glaciers in the Katian, then subsequent glacial melting (Saltzman and Young Reference Saltzman and Young2005; Trotter et al. Reference Trotter, Williams, Barnes, Leécuyer and Nicoll2008; Finnegan et al. Reference Finnegan, Bergmann, Eiler, Jones, Fike, Eisenman, Hughes, Tripati and Fischer2011; Melchin et al. Reference Melchin, Mitchell, Holmden and Štorch2013; Pohl et al. Reference Pohl, Donnadieu, Le Hir, Ladant, Dumas, Alvarez-Solas and Vandenbroucke2016a) that all culminated in the LOME in the latest Katian and Hirnantian ages. Throughout this time slice, especially in the Sandbian and Katian, blastozoan species exhibit speciation events regardless of these environmental perturbations, until the time before the LOME (Fig. 5). Interestingly, there is one speciation event that occurs between the two major Hg/total organic carbon (TOC) enrichments (Jones et al. Reference Jones, Martini, Fike and Kaiho2017; Fig. 5). Modeling studies have indicated that emplacement and subsequent weathering of a LIP could lead to 5°C cooling within 2 Myr of emplacement within the Late Ordovician Katian Age (Lefebvre et al. Reference Lefebvre, Servais, François and Averbuch2010). If gradual cooling does lead to increased speciation, as other evolutionary studies during times of Cenozoic global cooling have proposed (e.g., Fraass et al. Reference Fraass, Kelly and Peters2015; Davis et al. Reference Davis, Hill, Astrop and Wills2016; Baarli et al. Reference Baarli, Malay, Santos, Johnson, Silva, Meco, Cachão and Mayoral2017; Lam and Leckie Reference Lam and Leckie2020; A. R. Lam and R. M. Leckie unpublished data), we would predict speciation events occurring within 2 Myr of LIP activity. Indeed, we do see this occur in this blastozoan phylogeny and for the brachiopods that exhibit vicariance events (Fig. 5) after the first Mg/TOC pulse in the latest Katian. We caution that we have limited data on which to base this hypothesis, but propose that additional data and studies with more refined stratigraphic control on species temporal ranges and occurrences could reveal more interesting evolutionary patterns in response to the Late Ordovician LIP events.
It is clear from the phylogenetic hypothesis of blastozoans presented here that this group was hit hard during the first pulse of the LOME at the Katian/Hirnantian boundary, with several species going extinct just before or at the boundary (Fig. 2). However, one grouping of diploporans, the holocystitids, continues across the LOME into the Silurian, with reconstructed ancestral ranges indicating strong affinities to the Cincinnati Basin (Fig. 2A). Further work should explore the trends and drivers of extinction selectivity across the LOME in this particular group.
Comparison of Dispersal Paths and Speciation Types
Few other studies have utilized BioGeoBEARS to reconstruct ancestral relationships and infer patterns of speciation and dispersal, but here we compare our results with those of Lam et al. (Reference Lam, Stigall and Matzke2018) and Congreve et al. (Reference Congreve, Krug and Patzkowsky2019). Our results indicate that speciation events in the Ordovician were dominated by dispersal, a pattern that is similar to brachiopod clades in the Lam et al. (Reference Lam, Stigall and Matzke2018) study. Speciation events for the diploporans and other blastozoan species are scattered throughout the Ordovician, whereas brachiopod speciation events are clustered within the Late Ordovician, and trilobite speciation events are clustered within the Floian Age and Dapingian to Darriwilian ages, which are times of global cooling (Fig. 5).
The three groups (blastozoans, brachiopods, and trilobites) converge in their similarity of type of speciation through the Ordovician. Among all the models for the groups, models that incorporated the + j parameter improved model fit, suggesting that jump dispersal was a speciation mode utilized by several marine invertebrates. Thus, the finding from our study and those of Lam et al. (Reference Lam, Stigall and Matzke2018) and Congreve et al. (Reference Congreve, Krug and Patzkowsky2019) support the hypothesis of Stigall et al. (Reference Stigall, Edwards, Freeman and Rasmussen2019) that several volcanic island arcs allowed for increased long-distance dispersal events and likely primed the GOBE event that occurred in the Darriwilian Age. Additional BioGeoBEARS results from Bauer (Reference Bauer2020) for echinoderm taxa also indicate models with the + j parameter fit the data best. Therefore, jump dispersal may be a prevalent mode of speciation during several times in the Paleozoic.
Lam et al. (Reference Lam, Stigall and Matzke2018) and Congreve et al. (Reference Congreve, Krug and Patzkowsky2019) also recovered several dispersal events between Laurentia and Baltica from the Dapingian to Sandbian ages for trilobites and brachiopods. Interestingly, the majority of the Lam et al. (Reference Lam, Stigall and Matzke2018) dispersal paths are from Laurentia to Baltica, whereas those reconstructed in this analysis are bidirectional between these areas (Figs. 3, 4). For brachiopods and trilobites, there are a number of dispersal events, mainly from Laurentia to Gondwana, during the Middle Ordovician (Lam et al. Reference Lam, Stigall and Matzke2018) and between the two areas (Congreve et al. Reference Congreve, Krug and Patzkowsky2019). This study indicates dispersals of these blastozoan species between Laurentia and Gondwana were more restricted during this interval but still prevalent (Fig. 4).
There are several differences and similarities when Ordovician dispersal paths are compared within Laurentia across clades. In the Middle Ordovician, trilobite dispersal events occur among most Laurentian basins, with no preference for direction or any one basin (Lam et al. Reference Lam, Stigall and Matzke2018). The Early to early Middle Ordovician dispersal paths for the blastozoans are very similar, with rampant paths recovered among all Laurentian basins with no perceivable preference for direction. In Lam et al. (Reference Lam, Stigall and Matzke2018), Late Ordovician Laurentian dispersal and vicariance events for brachiopods and trilobites retain a random pattern among the Laurentian southern and western basins, with a directionality preference from interior basins into northern and southern basins. The Late Ordovician blastozoan dispersal paths are similar, except dispersal into the southern Laurentian Basin (Fig. 4) occurs bidirectionally between basins.
In summary, there are several differences and similarities among dispersal paths across clades during the Ordovician. Comparisons are made difficult across studies by the inconsistent delineation of basins and areas as restricted by the occurrences of species and the time slicing of dispersal paths discussed. However, it is clear that jump dispersal improves biogeographic model fit for invertebrate clades, and the exchange of taxa between Baltica and Laurentia was dominant throughout the Ordovician Period.
Conclusions and Future Work
This study is the first to infer Ordovician diploporan blastozoan speciation and dispersal events within a rigorous statistically informed phylogenetic framework. Probabilistic inference of ancestral ranges was conducted within the R package BioGeoBEARS to determine the origin for these organisms and most likely mode of speciation through time. The number and percent of dispersal events was estimated using BSM to determine how dispersal paths changed across three time slices through the Ordovician. The first time slice includes the Early to early Middle Ordovician; the Darriwilian Age time slice includes the GOBE; and the Late Ordovician time slice includes several Earth system perturbations and the LOME.
The results of the BioGeoBEARS analysis indicate that all models that included the + j parameter for jump dispersal improved model fit, with the best-fit model being DIVALIKE + j. This suggests that jump dispersal was an important speciation mode for these species. From the DIVALIKE + j model, there are several potential source regions for the earliest ancestors of the diploporans, with true diploporans appearing in the Early Ordovician and likely originating in either Baltica or Gondwana, potentially during the late Cambrian. Results from the BSM analysis indicate that founder-event speciation accounted for nearly half of all speciation events within all three time slices.
Directionality of dispersal events was inferred from both range expansion and jump-dispersal speciation. Within the Early to early Middle Ordovician time slice, the highest number of dispersal events was recovered, with several occurring among Laurentia, Baltica, and Gondwana. Several events were also recovered among Laurentian basins. The Middle Ordovician time slice, characterized by the GOBE, includes the least amount of recovered dispersal events, but with important dispersal events between Baltica and the Cincinnati Basin. The Late Ordovician again shows an increased number of dispersals among Laurentia, Baltica, and Gondwana, with dispersals most prevalent on the southern edge of the Laurentian paleocontinent. These findings indicate that distance did not impede dispersal, and species likely utilized the major current and gyre systems that flowed between continents. Intracontinental Laurentian dispersal events were likely mediated by surface currents and strong storm activity, and perhaps the lack of tectonically emplaced thermal and physical barriers that did not develop until the Middle Ordovician.
Compared with geochemical data and model results for the Ordovician, blastozoan dispersal and jump dispersal do not cluster during times of global cooling or atmospheric oxygenation, as has been found for brachiopods and trilobites. This indicates that the blastozoans may have had different environmental tolerances than other marine invertebrates and did not respond similarly to abiotic factors. Further analyses may also help to understand how biogeographic patterns and oxygen levels may have played different roles in blastozoan evolutionary patterns, as compared with other invertebrate groups. A larger dataset of blastozoan taxa with additional and broader occurrence data and improved stratigraphic and temporal ranges will help reveal further drivers of speciation and dispersal during the early Paleozoic.
We acknowledge that a geographic bias in sampling toward Baltica and Laurentia exists across many invertebrate taxa, with Gondwana and other paleocontinents being less broadly explored (Miller Reference Miller2000; Tarver et al. Reference Tarver, Braddy and Benton2007; Sumrall and Zamora Reference Sumrall and Zamora2011; Lefebvre et al. Reference Lefebvre, Sumrall, Shroat-Lewis, Reich, Webster, Hunter, Nardin, Rozhnov, Guensberg and Touzeau2013). This biases overall interpretations of paleoecological and biogeographic trends, as Laurentia and South China were primarily situated in more tropical, carbonate regimes, whereas Gondwanan depositional environments were more siliciclastic during the Ordovician. More accurate paleobiogeographic trends will be assessed with better global sampling, as new fossil species incorporated into phylogenetically informed paleobiogeographic studies can drastically change interpretations of diversity and biogeography (Lefebvre et al. Reference Lefebvre, Ghobadipour and Nardin2005; Sumrall et al. Reference Sumrall, Deline, Colmenar, Sheffield, Zamora, Zamora and Rábano2015; Zamora et al. Reference Zamora, Sumrall, Zhu and Lefebvre2016). Future work should focus on sampling more global occurrences of blastozoan echinoderms to be included in future biogeographic studies, so that we might be able to better assess the abiotic factors driving Paleozoic echinoderm speciation and biogeography.
Acknowledgments
We would like to thank J. E. Bauer and M. R. Limbeck for reading earlier versions of this article. We thank B. Deline and an anonymous reviewer for critical comments and suggestions that greatly strengthened this contribution. A.R.L. was supported by a Cushman Foundation for Foraminiferal Research Johanna M. Resig Foraminiferal Research Fellowship and the Binghamton University Presidential Diversity Postdoctoral Fellowship. S.L.S. was supported by the Paleontological Society Arthur James Boucot Research Grant. N.J.M.'s development of BioGeoBEARS and BEASTmasteR were supported by the National Institute for Mathematical and Biological Synthesis, an institute sponsored by the National Science Foundation (NSF), the U.S. Department of Homeland Security, and the U.S. Department of Agriculture through NSF award no. EFJ0832858, with additional support from the University of Tennessee, Knoxville. N.J.M. was also funded by the Australian Research Council's Discovery Early Career Researcher award no. DE150101773, and by the Australian National University. He is currently supported by the University of Auckland and New Zealand Marsden Grants 16–UOA–277 and 18–UOA–034. Funding for publication was kindly provided by the University of South Florida Libraries. Supplementary Files 1 and 2 can be found at https://doi.org/10.5061/dryad.4tmpg4f6j.
Open access
