Estimating dispersal and evolutionary dynamics in diploporan blastozoans (Echinodermata) across the great Ordovician biodiversification event

Abstract

Echinoderms make up a substantial component of Ordovician marine invertebrates, yet their speciation and dispersal history as inferred within a rigorous phylogenetic and statistical framework is lacking. We use biogeographic stochastic mapping (BSM; implemented in the R package BioGeoBEARS) to infer ancestral area relationships and the number and type of dispersal events through the Ordovician for diploporan blastozoans and related species. The BSM analysis was divided into three time slices to analyze how dispersal paths changed before and during the great Ordovician biodiversification event (GOBE) and within the Late Ordovician mass extinction intervals. The best-fit biogeographic model incorporated jump dispersal, indicating this was an important speciation strategy. Reconstructed areas within the phylogeny indicate the first diploporan blastozoans likely originated within Baltica or Gondwana. Dispersal, jump dispersal, and sympatry dominated the BSM inference through the Ordovician, while dispersal paths varied in time. Long-distance dispersal events in the Early Ordovician indicate distance was not a significant predictor of dispersal, whereas increased dispersal events between Baltica and Laurentia are apparent during the GOBE, indicating these areas were important to blastozoan speciation. During the Late Ordovician, there is an increase in dispersal events among all paleocontinents. The drivers of dispersal are attributed to oceanic and epicontinental currents. Speciation events plotted against geochemical data indicate that blastozoans may not have responded to climate cooling events and other geochemical perturbations, but additional data will continue to shed light on the drivers of early Paleozoic blastozoan speciation and dispersal patterns.

Introduction

Echinoderms show complex patterns of morphological disparity and taxonomic diversity throughout their evolutionary history (Sumrall 1997; Deline 2015; Deline and Thomka 2017; Deline et al. 2020) and have been suggested to respond to long-term environmental and climatic patterns (Paul 1968; Clausen 2004; Dickson 2004; Clausen and Smith 2005, 2008; Zamora and Smith 2008; Rahman and Zamora 2009). The early Paleozoic experienced large shifts in geochemical cycles (Haq and Schutter 2008; Trotter et al. 2008; Rasmussen et al. 2016; Albanesi et al. 2019) that impacted microfossil and marine invertebrate groups (reviewed in Stigall et al. 2019), 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. 2013); 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 1973). However, it is suspected that many groups of blastozoans are not monophyletic and that respiratory structures might be convergent (Paul 1988; Sumrall 2010; Sheffield and Sumrall 2017, 2019a). Sheffield and Sumrall (2019a), 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 2007; Nardin and Lefebvre 2010; Lefebvre et al. 2013; Zamora et al. 2013), 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 2010). 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 1997; Wright and Stigall 2013; Trubovitz and Stigall 2016; Lamsdell et al. 2017; Stigall et al. 2017). 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 2020) 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. 2018), by evolving these respiratory structures and other morphological features.

In this contribution, we use the R package BioGeoBEARS (Matzke 2013) 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. 2017). 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. 2016) 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 1970; Pedersen et al. 1992; Harper et al. 1996; Glen et al. 1998; Cocks and Torsvik 2011; Samygin 2019; Zhang et al. 2019), especially within the Iapetus Ocean (Rasmussen and Harper 2011; Liljeroth et al. 2017). 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 1991; van de Pluijm et al. 1995; Harper et al. 1996; Cocks and Torsvik 2011). 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 1987; Soper et al. 1992; Cocks and Torsvik 2011; Zagorevski and Van Staal 2011).

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 2010). 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 1982; Ettensohn 2010). This narrow basin was characterized by major clastic wedges (Ettensohn 2004) that did not interrupt tropical carbonate sedimentation across the southern margin of Laurentia (Holland and Patzkowsky 1997). 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. 1991; Ettensohn 2010). 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. 2001). 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 1991; Kolata et al. 2001; Young et al. 2015; Quinton et al. 2017). With the initiation of cool oceanic waters into Laurentia and expansion of the Sebree Trough (Young et al. 2015), warm carbonate platform deposition was replaced with mixed cool-water carbonate and clastic facies throughout eastern central Laurentia (Keith 1989; Ettensohn et al. 2002). Warm-water carbonate facies resumed during the Late Ordovician Katian Age (Ettensohn et al. 2002).

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. 2008; Albanesi et al. 2019) and high global sea level, punctuated by shorter-term sea-level falls (Haq and Schutter 2008). 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. 2008; Finnegan et al. 2011; Albanesi et al. 2019). 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. 2009; Turner et al. 2012; Dabard et al. 2015; Amberg et al. 2016; Pohl et al. 2016a; Rasmussen et al. 2016). Cooling continued into the Late Ordovician with well-developed continental ice sheets over west Gondwana responding to orbital forcing (Beuf et al. 1971; Hambrey 1985; Crowley and Baum 1991; Eyles 1993; Crowell 1999; Scotese et al. 1999; Sutcliffe et al. 2000; Sinnesael et al. 2019).

Modeling experiments under different atmospheric carbon dioxide (CO₂) 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. 1999; Herrmann et al. 2004; Pohl et al. 2016b; 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 2003; Pohl et al. 2018).

Figure 1. 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 (2013), with the approximate locations of currents and gyre systems from Pohl et al. (2016b).

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. 2004; Harper 2006; Stigall 2018; Stigall et al. 2019). At the family level, diversity increased on average threefold, whereas generic diversity increased fourfold. In a recent review by Stigall et al. (2019), 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 2017).

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 1984; Trotter et al. 2008; Saltzman et al. 2015; Pohl et al. 2016a; Rasmussen et al. 2016; Trubovitz and Stigall 2016; Young et al. 2016; Edwards et al. 2017; Lam et al. 2018). Increased atmospheric oxygen levels could have contributed to increased metabolic processes for marine invertebrates and allowed for increased calcite precipitation (Pruss et al. 2010). 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. 1996; Liljeroth et al. 2017; Lam et al. 2018). 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. 2017).

Concurrent with the latest Ordovician was the LOME, which occurred in two major pulses (Brenchley et al. 1994; Brenchley 1995; Brenchley and Marshall 1999), 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. 2011). 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. 2012, 2016). 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. 2020). 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. 1994; Sheehan 2001; Melchin et al. 2013). 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. 2017).

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 1974, 1986; Thompson and Satterfield 1975; Elias 1983; Barrick 1986; Mohibullah et al. 2012). 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. (2018) 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. (2018).

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. 2013, 2015; Zamora et al. 2016; Sheffield et al. 2017; Sumrall and Zamora 2018) 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 (2019a), 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.

Figure 2. 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. 2016). Refer to Supplementary Table 1 for species’ groupings within Echinodermata. Fossil image of diploporan Eumorphocystis multiporata Branson and Peck 1940, SUI 97599, modified from Sheffield and Sumrall (2019b).

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. 2016) 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 (2019a) 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, 1965 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 (2019a).

Biogeographic models require a fully bifurcating, time-calibrated phylogenetic tree. The polytomies present in Sheffield and Sumrall (2019a) 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 (2019a), 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 2019b). Asteroblastid (Asteroblastus stellatus) groups away from other diploporans. In Sheffield and Sumrall (2019a), 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 2016) in Beast2 (Bouckaert et al. 2014). BEASTmasteR (Matzke 2019) 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 2001) was implemented in the Beast2 XML via BEASTmasteR (Matzke and Irmis 2018). 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. (2018) 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 2013). 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 2008); dispersal–vicariance–analysis (DIVA; Ronquist 1997); and BayArea (Landis et al. 2013). 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 2013), 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 2014). This mode of speciation has been found to be important for many extant clades (Matzke 2014; Dupin et al. 2017), as well as Ordovician brachiopods and trilobites (Lam et al. 2018; Censullo and Stigall 2019; Congreve et al. 2019).

Although a recent critique argues that DEC and DEC + j cannot be statistically compared (Ree and Sanmartín 2018), that critique can be rejected on several grounds (Klaus and Matzke 2019). For example, the DEC/DEC + j statistical comparison was validated with simulation-inference tests in Matzke (2014), but these tests were not rebutted, or even acknowledged, by Ree and Sanmartín (2018). 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. 2019).

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 2002). Following O'Meara et al. (2006), 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. 2017). Stochastic mapping techniques were first developed for discrete characters like morphology and DNA (Nielsen 2002; Huelsenbeck et al. 2003; Bollback 2006). Instead of calculating the probability of each character state at each node (traditional ancestral state probabilities; Felsenstein 2003), 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. 2017).

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. 2018; Censullo and Stigall 2019; Congreve et al. 2019).

Table 1. 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. 2017; Lam et al. 2018; Censullo and Stigall 2019; Congreve et al. 2019).

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.

Table 2. 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.

Figure 3. 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. 2017). 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. 2015), while other studies question whether long-distance dispersal in Paleozoic echinoderms was possible (e.g., Rozhnov 2013); 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. (2016b) for the Middle Ordovician (460 Ma; Fig. 1) under 8 times preindustrial atmospheric CO₂ levels (PAL). From the Pohl et al. (2016b) 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 1991; 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 2011).

Figure 4. 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. 2008; Pohl et al. 2016a; Edwards et al. 2017; Quinton et al. 2018), 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. 2016b: 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. 2016b).

Figure 5. 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 (2008); generalized δ¹³C carbonate curve from Bergström et al. (2009), plotted with mercury (Hg) over total organic carbon (TOC) from Monitor Range, Nevada (gray line) from Jones et al. (2017); δ¹⁸O curve from Albanesi et al. (2019); red atmospheric CO₂ curve from the GEOCARBSULF model of Krause et al. (2018); black curve of oxygen-dependent atmospheric CO₂ model from Edwards et al. (2017); and blue curve of atmospheric CO₂ from the COPSE biogeochemical model of Lenton et al. (2018). Brachiopod and trilobite speciation events from Lam et al. (2018). 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 2010), 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 2010), with major hurricanes sweeping across the paleocontinent from 10° to 30° north and south of the equator (Jin et al. 2013). As concluded in other studies that have focused on intracratonic Laurentian dispersal events of marine invertebrates (Wright and Stigall 2013; Bauer and Stigall 2014; Lam and Stigall 2015; Lam et al. 2018), 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. 2015; Pohl et al. 2016a; Edwards et al. 2018; Albanesi et al. 2019). 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. 2016; Lam et al. 2018; Stigall et al. 2019), in part leading to the continued low endemicity of blastozoans through the Darriwilian Age (Nardin and Lefebvre 2010). 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 2011). With increased cooling (Fig. 5) and intensification of currents and upwelling due to reduced temperatures and increased latitudinal temperature gradients (Pohl et al. 2016b; Rasmussen et al. 2016; Lam et al. 2018), 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 2010). 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. 2015) 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. 2019). 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. 2016a) 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 2010). Diversity curves for the blastozoans at the genus level in other studies do suggest an increase in diversity during the Darriwilian (Nardin and Lefebvre 2010). 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 2017; Cole 2019), 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. 2002), 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. 2017) 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 2005; Trotter et al. 2008; Finnegan et al. 2011; Melchin et al. 2013; Pohl et al. 2016a) 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. 2017; 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. 2010). 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. 2015; Davis et al. 2016; Baarli et al. 2017; Lam and Leckie 2020; 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. (2018) and Congreve et al. (2019). 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. (2018) 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. (2018) and Congreve et al. (2019) support the hypothesis of Stigall et al. (2019) 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 (2020) 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. (2018) and Congreve et al. (2019) 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. (2018) 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. 2018) and between the two areas (Congreve et al. 2019). 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. 2018). 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. (2018), 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 2000; Tarver et al. 2007; Sumrall and Zamora 2011; Lefebvre et al. 2013). 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. 2005; Sumrall et al. 2015; Zamora et al. 2016). 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.