Dissecting the paleocontinental and paleoenvironmental dynamics of the great Ordovician biodiversification

Abstract

The Ordovician was a time of drastic biological and geological change. Previous work has suggested that there was a dramatic increase in global diversity during this time, but also has indicated that regional dynamics and dynamics in specific environments might have been different. Here, we contrast two paleocontinents that have different geological histories through the Ordovician, namely Laurentia and Baltica. The first was situated close to the equator throughout the whole Ordovician, while the latter has traversed tens of latitudes during the same time. We predict that Baltica, which was under long-term environmental change, would show greater average and interval-to-interval origination and extinction rates than Laurentia. In addition, we are interested in the role of the environment in which taxa originated, specifically, the patterns of onshore–offshore dynamics of diversification, where onshore and offshore areas represent high-energy and low-energy environments, respectively. Here, we predict that high-energy environments might be more conducive for originations.

Our new analyses show that the global Ordovician spike in genus richness from the Dapingian to the Darriwilian Stage resulted from a very high origination rate at the Dapingian/Darriwilian boundary, while the extinction rate remained low. We found substantial interval-to-interval variation in the origination and extinction rates in Baltica and Laurentia, but the probabilities of origination and extinction are somewhat higher in Baltica than Laurentia. Onshore and offshore areas have largely indistinguishable origination and extinction rates, in contradiction to our predictions. The global spike in origination rates at the Dapingian/Darriwilian boundary is apparent in Baltica, Laurentia, and onshore and offshore areas, and abundant variability in diversification rates is apparent over other time intervals for these paleocontinents and paleoenvironments. This observation hints at global mechanisms for the spike in origination rates at the Dapingian/Darriwilian boundary but a domination of more regional and local mechanisms over other time intervals in the Ordovician.

Introduction

The numerical increase of marine orders, families, and genera in the Ordovician was so dramatic (Sepkoski et al. 1981; Miller 1997a,b, 2012; Sepkoski 1997; Harper 2006; Bassett et al. 2007; Alroy et al. 2008) that it has been termed the great Ordovician biodiversification event (GOBE; Webby et al. 2004; Servais and Harper 2018). This increase is thought to be especially rapid around the Dapingian/Darriwilian boundary (Servais et al. 2009; Hints et al. 2010; Rasmussen et al. 2016; Trubovitz and Stigall 2016), but factors influencing this increase remain obscure.

While the hypotheses for what these factors might be are many and varied, it is possible to group them into two classes. The first class of hypotheses is climate related. Specifically, a mid-Ordovician cooling is thought to have contributed to a more favorable global climate regime that in turn promoted diversification (Trotter et al. 2008; Rasmussen et al. 2016). The second class of hypotheses involves nutrient availability. For instance, periods of increased tectonics (Miller and Mao 1995) and/or tectonically induced volcanic activity (Botting 2002) are hypothesized to have led to more sedimentary and nutritional input into the oceans and habitat fractioning, both of which may have enhanced the diversification of taxa (Miller and Mao 1995). The establishment of nutrient-rich upwelling zones has been documented from Laurentia during the Middle Ordovician (Pope and Steffen 2003). These upwelling zones may have served as new ecospace (Rasmussen et al. 2016), and their subsequent occupation by both migrants and taxa that evolved in situ could have contributed to an increase in taxon diversity. Some of the mechanisms that were suggested to have an influence on the GOBE, such as tectonic activity and volcanism, likely acted at a regional rather than the global scale, and such regional changes are unlikely to be simultaneous (Miller 1997a,b, 2004; Zhan and Harper 2006). Therefore the influence of the different mechanisms should be examined separately for the different paleocontinents (Miller 1997a; Miller and Mao 1998) and paleoenvironments (e.g., Miller and Mao 1995; Miller and Connolly 2001; Novack-Gottshall and Miller 2003).

Continents drifted and experienced different spatial environments over the Ordovician. Whereas the center of the paleocontinent Laurentia was situated close to the equator throughout the Ordovician, Baltica rotated northward, starting in the southern mid-latitudes and ending south of the equator by the Late Ordovician (Cocks and Torsvik 2006; Torsvik and Cocks 2016a). These two continents can act as model systems for a relatively stable continent (Laurentia) versus a continent on which the physical environmental is relatively unstable due to continental movement (Baltica) during the time in question. In this paper, we explore whether greater environmental change on Baltica is associated with greater taxon turnover but lower genus richness compared with Laurentia.

Different paleoenvironments may also influence diversification rates. For instance, morphological novelties, as represented by ordinal-level originations, arose preferentially in onshore rather than offshore areas (Jablonski et al. 1983; Jablonski 2005). Based on his work on ordinal-level originations, Jablonski (Jablonski et al. 1983; Jablonski 2005) predicted that genus originations would preferentially occur where their higher-level taxa were already established. If this is true, we expect that relative origination dynamics might be similar for onshore and offshore areas, all things being otherwise equal. However, offshore areas are more shielded from environmental changes than onshore areas, where changes in sea levels lead to changes in the size of habitable areas (Holland and Christie 2013). Here, we explore whether greater environmental variability in onshore areas is associated with higher genus origination and extinction rates compared with offshore areas during the Ordovician.

We present Ordovician genus origination and extinction rates based on capture–recapture models for the revised global Ordovician stages (Bergström et al. 2009; Harper and Servais 2013; Lindskog et al. 2017) both on a global scale, for the paleocontinents Baltica and Laurentia, and for onshore and offshore areas.

Methods and Data

Data

We downloaded data from the Paleobiology Database (PBDB) on 26 January 2018 (https://paleobiodb.org/data1.2/occs/list.csv?taxon_reso=genus&interval=Cambrian,Aeronian&show=class,acconly,ecospace,coll,coords,loc,paleoloc,lithext,geo,refattr). These data consist of 130,367 occurrences of taxa identified to at least the genus level that span the Cambrian to the Aeronian (second stage in the Silurian) and their metadata. Only accepted genus names were included in our data analyses (73,735 occurrences). Although our analyses are focused on the Ordovician, we included Cambrian and Silurian data to ameliorate edge effects (see subsection on diversification rates).

The temporal resolution associated with the data we retained from the PBDB download range from 0.4 to 66.2 Myr. Considering only Ordovician time intervals, the resolution is from 1.7 to 50.8 Myr (median: 7.8 Myr; mean: 8.6 Myr). To increase the temporal resolution of our analyses, we only included PBDB data that were assigned to time bins smaller than 12 Myr (see Supplementary Material text and Supplementary Fig. S1). For occurrences reported in PBDB with regional Ordovician stages, we randomly assigned point estimates using a uniform distribution between the reported minimum and maximum ages and then assigned those estimates to global Ordovician stages (median: 7.7 Myr; mean: 5.9 Myr). This was to ensure comparability to standard, global treatments of Ordovician data. We present only one iteration of parameter estimates in our main figures. But because random assignments to time intervals vary, we also present averaged parameter estimates from 100 iterations of such random assignments to the global data set or subsets for each of our analyses in the Supplementary Material.

We created a non-observation/observation matrix from the downloaded PBDB data. In this matrix, each row of data is a string of observations (1) or non-observations (0) of a genus within sequential global stages. We estimated extinction, origination, and sampling rates by applying a capture–recapture modeling framework to these data. Although there are many different models within the capture–recapture modeling framework (Nichols and Pollock 1983; Liow and Nichols 2010), we follow previous paleontological applications (Connolly and Miller 2001b; Liow et al. 2015; Kröger 2017) in using the Pradel seniority model (Pradel 1996) for estimating diversification rates. In addition, we use the POPAN model (Schwarz and Arnason 1996) for estimating genus richness.

We note that up to 65% of our global genus data set is made up of arthropods, brachiopods, and mollusks. The remaining 35% of the phyla in the data set are bryozoans, chordates, cnidarians, echinoderms, hemichordates, and poriferans.

Diversification Rates

The Pradel seniority model (Pradel 1996) was developed for estimating population parameters, namely survival probability φ, seniority probability γ, sampling probability p, and population growth λ. Because our data represent genera rather than individuals, these parameters reflect, respectively, genus survival probability (the complement of which is genus extinction probability), genus “seniority” probability (the complement of which is genus origination probability), sampling probability, and diversification rates. Very briefly, we assume that genera are sampled within time intervals, i. We number these from older (i) to younger geological intervals (i + 1). In this framework, survival (φ) is the probability that a genus that is extant in a given time interval continues to be extant in the next time interval. Extinction probability (1 − φ) is hence the probability of a genus going extinct in the next time interval, given that it was extant in the time interval in question. Similarly, seniority probability γ describes the probability that if a genus was extant in time i + 1, it was already extant in time i. In this case, 1 − γ is the origination probability between the two time intervals. Diversification λ describes how much a group of taxa is growing or declining from one time interval to the next, where the net diversification rate is λ_i − 1. If the number of genera increases from one time interval (i) to the next (i + 1), the net diversification rate will be positive, while if the number of genera is decreasing, it will be negative. A value of zero indicates no net increase or decrease in the number of genera. Given that the Pradel seniority model has been previously detailed for paleontological data sets (Connolly and Miller 2001a,b, 2002; Liow and Nichols 2010; Liow et al. 2015), we refer readers to those publications for further details. It is sufficient to note that a logit link is used and that the parameters are estimated in a maximum-likelihood framework.

As already mentioned, the basic parameters of the Pradel seniority model are survival, seniority, and sampling probabilities. These basic parameters can be specified in different ways. For example, we can assume that they are constant through time (time-constant model). Using origination probability as an example, this means that time interval i has the same origination probability as time interval i + 1, and so on, such that there is only one estimate for the whole of the Ordovician. Because event probabilities increase as the sampling intervals increase, we specify φ^t or γ^t, with t being the duration between two sampling occasions such that origination and extinction probabilities can be interpreted in this time-constant model. Another model we could use is a fully time-varying one, meaning that estimates (e.g., for origination) are allowed to be different for every time interval. Here, to allow comparisons of our estimates across these stages, we transform the estimated probabilities into rates by using a Poisson model. Specifically, Φ_i = −log[1−(1−φ)_i]/T _i for the extinction rate Φ_i, Γ_i = −log[1−(1−γ)_i]/T _i for the origination rate Γ_i, and P _i = −log[1 − (1 − p)_i]/T _i for the sampling rate P _i (Liow et al. 2015), with T _i being the duration of one time interval. Additionally, we can add covariates to the model, such as affiliation with Laurentia or Baltica, in order to estimate separate parameters for the two groups simultaneously. We specified and examined many models (see “Results”) but present only a selection of these in our main text (see Supplementary Material for additional results).

In a fully time-varying Pradel seniority model, survival and sampling in the last time interval (given by p _iφ_i−1, with i in this case being the last time interval) and seniority and sampling in the first time interval are not separately identifiable (Connolly and Miller 2001b). To obtain estimates for survival, seniority, and sampling probabilities for all of the Ordovician, we include two Cambrian time bins before and two Silurian after the first and last Ordovician stages, respectively. In that way we avoid the described boundary effects. Note that for the time-constant model, we only use occurrence data from the Ordovician, because there are no boundary effects and because we are only interested in comparing Ordovician probabilities.

Genus Richness

We estimated genus richness using the POPAN model (Schwarz and Arnason 1996). In contrast to the closed-population approach of the Pradel seniority model, the POPAN model is an open-population approach (Schwarz and Arnason 1996). This means that individuals (genera in our case) can enter the observed population (global or regional community in our case) from a superpopulation. In our context, the superpopulation (N) is the number of genera available to be sampled that may or may not be sampled throughout the whole study period, where $N = \sum\nolimits_0^{i-1} {B_i} $ (Schwarz and Arnason 1996). Here, B _i is the number of new genera originating between two time intervals (number of births in the original description of POPAN). It is estimated by B _i = pent _i−1N, where pent is the probability of entering the pool of genera. The genus richness in each time interval is then N _i = N _i−1φ_i−1 + B _i (Schwarz and Arnason 1996), where φ_i−1 is the survival probability into the time interval in focus. Similar to the edge effect described above, φ_i−1 cannot be estimated independently in a fully time-varying POPAN model. This is also true for the parameters N ₁, N _max, pent ₁, and pent _max−1. Thus, we include additional time bins before and after the Ordovician, as done for the Pradel seniority model.

Laurentia and Baltica

To assign the observations of genera documented in the PBDB to paleocontinents, we looked up their geoplate code from Müller et al. (2008), which is provided in the PBDB, and then assigned these to either Baltica or Laurentia using Cocks and Torsvik (2004) and Torsvik and Cocks (2016a,b) for guidance. Note that in some cases, geoplate codes were not assigned in the PBDB. In those cases, we used country (cc) or state, encoding modern positions, for assignments to paleocontinents in PBDB. We provide our assignments in the Supplementary Material. Some genera were observed as “endemic” to Baltica (194 genera) or Laurentia (731 genera), while others were described from both paleocontinents (480 genera). We restricted our analyses to those genera that were observed only on either one of the paleocontinents in question. “Endemic” is in quotation marks, because we assume that genera are found only where they are observed, and we do not explicitly model migration.

We note that 33% of the Baltic genera are arthropods, while 18% are brachiopods and 16% are mollusks. The remaining 33% of the Baltic genera are bryozoans, echinoderms, cnidarians, chordates, annelids, hemichordates, and others. Of the Laurentian genera, 25% are arthropods, 21% are mollusks, and 19% are echinoderms. The remaining 35% of Laurentian genera are bryozoans, brachiopods, poriferans, cnidarians, chordates, and others. More detailed information on the representation of different phyla in the data subsets can be found in Supplementary Tables S6–S8 in the Supplementary Material.

Onshore–Offshore

For our comparison of onshore and offshore areas, we assign genera using information on depositional environments given in the PBDB (see Supplementary Material for details). Roughly speaking, we consider onshore environments as those above fair-weather wave base in high-energy environments and offshore environments as those below storm wave base representing low-energy environments. We only consider genera that have their first occurrence record in either onshore (532 genera) or offshore (396 genera) areas, regardless of the locations of their subsequent observations. This means that extinction rates that are estimated with these data are not representative for the actual extinction that happens in onshore and offshore areas, respectively, because taxa may have migrated and became extinct somewhere else. To estimate extinction rates for onshore and offshore areas, we considered genera that have their last occurrence records in either onshore (537 genera) or offshore (411 genera) areas, regardless of the locations of their prior observations. In doing so, we assume that these first and last paleoenvironments of observations are the same as those where they actually originated or went extinct.

The analyses were done with the program MARK (White 2016) via the RMark interface for R (Laake 2013).

Results

Global origination rates are similar to extinction rates between the Tremadocian and Dapingian (Fig. 1A). At the Dapingian/Darriwilian boundary, the global origination rate is much higher than the extinction rate, which also dropped perceptibly compared with the earlier Ordovician. Extinction and origination rates become similar again before the Hirnantian mass extinction (Fig. 1A).

Figure 1. Global genus diversification and genus diversity dynamics. A, Origination (solid line and solid circles, O) and extinction (dotted line and open circles, E) rates as events per million years (Myr). B, Net diversification rate (solid line and triangles) and genus richness (dotted line and rectangles). Horizontal gray line represents zero net diversification. C, Sampling events per million years (Myr). Vertical lines are 95% confidence intervals. Stage abbreviations: Tr, Tremadocian; Fl, Floian; Dp, Dapingian; Dw, Darriwilian; Sa, Sandbian; Ka, Katian; Hi, Hirnantian.

The resulting net diversification is negative between the Tremadocian and Dapingian, corresponding to the extinction rates being slightly greater than origination rates during these time intervals (Fig. 1B). Between the Dapingian and Darriwilian, net diversification dramatically increases, with the number of genera doubled. Net diversification is lower again thereafter (i.e., close to zero), before it becomes negative, driven by the high extinction rate during the Hirnantian mass extinction. Estimates of genus richness largely correspond to the changes in net diversification changes (Fig. 1B), despite different assumptions of “population closure” used to estimate diversification and richness (see “Methods”). For example, net diversification rate between the Dapingian and Darriwilian is about 1; hence it follows that the point estimate of the number of genera about doubled from 900 in the Dapingian to 1900 in the Darriwilian. After the increase in the Darriwilian, genus richness stays relatively high until the Hirnantian extinction. There is an average increasing trend in sampling rates throughout the Ordovician (Fig. 1C).

Origination probabilities for Baltica are 0.190 (95% CI 0.155, 0.230) per million years for the whole Ordovician, while Laurentian origination probabilities are 0.088 (95% CI 0.079, 0.097), given a Pradel seniority model that imposes constant origination, extinction, and sampling probabilities through time. Extinction probabilities, estimated with the same time-constant model, are 0.169 (95% CI 0.134, 0.211) for Baltica and 0.069 (95% CI 0.061, 0.080) for Laurentia. Net diversification probabilities are 0.026 (95% CI 0.014, 0.038) for Baltica and 0.020 (95% CI 0.014, 0.025) for Laurentia.

Baltica shows greater interval-to-interval origination rates than Laurentia throughout the Ordovician (Fig. 2A). At the Tremadocian/Floian and Dapingian/Darriwilian boundaries, origination rates are at their highest on Baltica. On Laurentia, origination rates are at their highest at the Cambrian/Tremadocian and Dapingian/Darriwilian boundaries. Note that variation in estimates and their 95% confidence intervals is much greater until the Dapingian/Darriwilian boundary (see also Supplementary Fig. S7). However, after the Dapingian/Darriwilian boundary, Baltic origination rates are consistently higher than Laurentian ones.

Figure 2. Genus diversification dynamics on Baltica and Laurentia. A, Origination rates as events per million years (Myr). B, Extinction rates as events per million years (Myr). C, Net diversification rates where the horizontal line represents zero net diversification. D, Sampling rates as events per million years (Myr). Vertical lines represent 95% confidence intervals in A to D. Baltic rates (B) are indicated by solid lines and circles, while Laurentian rates (L) are indicated by dashed lines and triangles. Stage abbreviations as in Fig. 1. Estimated probabilities or confidence intervals thereof that are one or zero are indications of poorly constrained estimates and are hence removed from the plots. Net diversification rates with confidence intervals spanning more than 10 are also removed for the same reason.

From the Tremadocian/Floian boundary and thereafter, extinction rates are higher on Baltica than on Laurentia (Fig. 2B). On both paleocontinents, extinction rates decrease on average from the Tremadocian/Floian boundary until the Darriwilian/Sandbian boundary. Thereafter they show an increasing trend. Note that the 95% confidence intervals for extinction rates are largely overlapping for the two paleocontinents until the Dapingian/Darriwilian boundary (cf. Supplementary Fig. S7). Thereafter, Baltic extinction rates are consistently higher than Laurentian extinction rates.

As a result of the extinction and origination rates (Fig. 2A,B), net diversification rates are greater on Baltica than on Laurentia at the Tremadocian/Floian boundary (Fig. 2C). At the Floian/Dapingian boundary, net diversification rates are slightly greater on Laurentia than on Baltica. Thereafter, net diversification rates are slightly greater on Baltica than on Laurentia, until the Katian/Hirnantian boundary, although 95% confidence intervals are overlapping between the two paleocontinents throughout.

Sampling rates increase on average from the Floian to the Hirnantian for Laurentia (Fig. 2D). Sampling rates for Baltica can only be estimated from the Darriwilian onward, during which sampling rates are lower on Baltica compared with Laurentia.

Baltica had relatively low genus richness in the Early and early Middle Ordovician (Tremadocian to Dapingian) (Fig. 3). Laurentia had greater genus richness than Baltica during the same time, although the trajectories of both paleocontinents are very similar. The peak in net diversification (cf. Fig. 2) is expressed as a fivefold increase of genus richness on Baltica and an almost threefold increase on Laurentia between the Dapingian and Darriwilian. This results in an increase from about 13 to 64 genera on Baltica, and from about 99 to 286 genera on Laurentia for the higher taxa we examined. Although the first substantial increase in genus richness happened at the Dapingian/Darriwilian boundary, net diversification rates close to or slightly above zero until the Katian/Hirnantian boundary (Fig. 2C) led to an Ordovician peak in genus richness during the Katian (Fig. 3).

Figure 3. Genus diversity dynamics on Baltica (B) and Laurentia (L). Estimated genus richness on Laurentia (dashed line, triangles) and Baltica (solid line, circles) based on the POPAN model. Stage abbreviations as in Fig. 1. Vertical lines represent 95% confidence intervals.

Origination and extinction rates in onshore and offshore areas are barely distinguishable throughout the Ordovician. Onshore origination probabilities are 0.069 (95% CI 0.060, 0.079) per million years for the whole Ordovician, while offshore origination probabilities are 0.060 (95% CI 0.052, 0.070), given a Pradel seniority model that imposes constant origination, extinction, and sampling probabilities through time. Extinction probabilities for the whole Ordovician, using the reversed data set but the same model, are 0.070 (95% CI 0.058, 0.085) for onshore areas and 0.085 (95% CI 0.069, 0.104) for offshore areas.

Interval-to-interval origination rates in onshore and offshore areas are very similar throughout most of the Ordovician, given a time-varying Pradel seniority model with separate estimates for onshore and offshore taxa (Fig. 4). After the Cambrian/Tremadocian boundary, onshore rates, which start higher than offshore rates, converge. At the Dapingian/Darriwilian boundary, origination rates are high in both paleoenvironments, reflecting the global dynamics at that time.

Figure 4. Genus origination and extinction rates in onshore and offshore areas. A, Origination rates as events per million years (Myr), based on first observation in either onshore or offshore areas. B, Extinction rates as events per million years (Myr) based on last observation in either onshore or offshore areas. Vertical lines represent 95% confidence intervals in A and B. Onshore rates (on) are indicated by solid lines and circles while offshore rates (off) are indicated by dashed lines and triangles. Stage abbreviations as in Fig. 1. Estimated probabilities or confidence intervals thereof that are one or zero are indications of poorly constrained estimates and are hence removed from the plots. Net diversification rates with confidence intervals spanning more than 10 are also removed for the same reason.

Until the Floian/Dapingian boundary, extinction rates are higher in onshore areas than in offshore areas (Fig. 4B). At the Dapingian/Darriwilian boundary, extinction rates from onshore and offshore areas are similar, with greatly overlapping confidence intervals. Thereafter, until the Katian/Hirnantian boundary, extinction rates are higher in offshore than in onshore areas.

For diversification and sampling rates, based on first and last observations in either onshore or offshore areas, which are not central for our hypothesis, see Supplementary Material (Supplementary Figs. S9 and S10).

In this study we are specifically interested in changes of diversification rates on different paleocontinents or onshore/offshore areas through time. Therefore, we choose to present results from the fully time-varying model for all parameters and covariates. For completeness, we tested the fit of different models with different parameter specifications through time by comparing their Akaike information criterion values. Results of these model comparisons can be found in the Supplementary Material (Supplementary Tables S1–S4).

Discussion

Global Ordovician Dynamics

We confirmed, by explicitly modeling incomplete sampling, that the increase of genus richness during the Ordovician was very rapid around the Dapingian/Darriwilian boundary, as suggested by many authors, including Servais et al. (2009), Hints et al. (2010), Harper et al. (2013), Rasmussen et al. (2016), and Trubovitz and Stigall (2016). However, this increase could have been driven by a lowered extinction rate and/or a higher origination rate. Here, we show for the first time, that global origination rates for the Ordovician peaked at the Dapingian/Darriwilian boundary, right before the spike in genus richness during the Darriwilian (Fig. 1). In addition, the extinction rate at the same boundary is slightly lower than in the previous boundary (Floian/Dapingian boundary). Net diversification rates were in fact close to or lower than zero before the Dapingian/Darriwilian boundary, after which they rocketed up with a doubling of genera (from c. 800 to c. 1900 genera; Fig. 1) during the Darriwilian. This doubling of genera is greater than what was inferred recently by Kröger (2017) using a similar approach. Kröger estimated a 1.4 times increase in genus richness over the same time period (see Kröger 2017), although his genus richness estimate for the Darriwilian (ca. 2000) is similar to ours (ca. 1900). In addition, Kröger's maximum genus richness for Ordovician stages stands at ca. 2200 genera during the Katian, in contrast to ours, which occurs at the Darriwilian. These differences, despite our using the same model for genus richness inference, might be due to differential temporal assignments of genus records from the PBDB via RNames in the Kröger study. RNames (Kröger and Lintulaakso 2017) is a dynamic interface that correlates selected stratigraphic opinions with the purpose of increasing stratigraphic resolution. Dividing stages into stage slices reduces the number of observations per time bin, possibly depleting some stage slices of data, so confidence intervals of the estimates would subsequently increase, especially for our paleocontinental and paleoenvironmental analyses. The increases in confidence intervals would be even more severe if observations are restricted to specific paleocontinents or bathymetric groups. In addition, the rules by which stages are subdivided into stage slices in RNames are currently not well-documented enough for our use. For these reasons, we estimated diversification rates on a coarser time resolution, but with better understood data for the benefit of more confident estimates. There is some concern that bias is introduced by assigning fossil occurrences to global stages when the occurrences were originally associated with regional stages. However, our sensitivity analyses show that reducing the data set to only fossil occurrences assigned to global stages does not change our qualitative inferences (Supplementary Fig. S6 in the Supplementary Material).

In their pioneering global study, Connolly and Miller (2001b) simultaneously estimated origination, extinction, and sampling probabilities for the Ordovician. They found that both origination and extinction probabilities have higher values approximately around the Darriwilian, although they found no exceptional diversification peak (Connolly and Miller 2001b), in contrast to what we found. While Connolly and Miller (2001b) used the same model we did, they estimated origination and extinction probabilities between time intervals of equal length, because probabilities of events over unequal time intervals cannot be directly compared. Rather than forcing these data to conform to geologically unnatural but equal time intervals, we used unequal but geologically natural Ordovician stages, but we converted probabilities into rates (see “Methods”). The general diversification trends estimated in Connolly and Miller may also differ from ours due the broader set of marine invertebrates used in our estimation. For completeness, we divided our data set into subsets of those taxa Connolly and Miller (2001a) studied (see Supplementary Figs. S12, S15–S17). We find that the dynamics of these taxa broadly conform to the global patterns we have estimated. However, others have shown that graptolites (Cooper et al. 2014) and echinoderms (Sprinkle and Guensburg 1995, 2004), unlike the taxa we investigated (arthropods, brachiopods, and mollusks; see Supplementary Material), show an increase in genus richness during the Early Ordovician. Given that the model used in Connolly and Miller (2001b) and this study explicitly accounts for sampling probabilities, we do not expect the addition of more observations (of correctly identified and dated samples) to change the estimates of mean extinction and origination considerably, but to increase the sampling probabilities while tightening the confidence intervals for extinction and origination.

As already mentioned, many hypotheses have been raised for why diversification rates were so high during the GOBE, some of which likely acted more globally and others more regionally. While our global estimates integrate over the processes that occur at local to global levels, details of regional dynamics give us clues as to what types of processes might be operating. In other words, if the GOBE in different regions is driven by different factors, we would expect diversification dynamic patterns to be different. For instance, if increased sedimentary input influenced diversification rates, the rates of taxa closest to volcanically and tectonically active areas may experience more changes than those in the oceans farther away, all other things being equal. To explore the contrast between global and regional controls of Ordovician diversification, we examine the dynamics on Laurentia and Baltica and in onshore and offshore areas and discuss these in relation to our global estimates.

Contrasting Laurentia and Baltica

It is striking that the peaks in origination rates happened on both Baltica and Laurentia between the Dapingian and the Darriwilian (Fig. 2), a process that led to the largest increase in genus richness over the Ordovician at the Darriwilian for both continents. Prior to this peak, origination and extinction rates appear to be uncoordinated on the two paleocontinents (see also Supplementary Fig. S7). After the Dapingian/Darriwilian boundary, changes in diversification dynamics appear to be more coordinated. Based on this observation, we suggest that we can subdivide Ordovician diversification into three different phases. In the first phase (the Early Ordovician: Cambrian/Tremadocian and Tremadocian/Floian boundaries), the variation in diversification rates seems relatively uncoordinated on Baltica and Laurentia. In the second phase (the Middle Ordovician: Floian/Dapingian and Dapingian/Darriwilian boundaries), origination rates peak dramatically on both paleocontinents. In the last phase (the Late Ordovician: Darriwilian/Sandbian to Katian/Hirnantian boundaries), diversification rates seem to change in concert on the two paleocontinents. The three phases hint at regional effects on diversification being stronger in the Early Ordovician, and global effects being stronger in the Middle to Late Ordovician. We note in passing that our onshore and offshore data subsets appear also to be more uncoordinated before the Dapingian/Darriwilian boundary than after that, just like the dynamics we estimated for the two paleocontinents.

In the Early Ordovician, Laurentia was stably situated close to the equator, while Baltica was rotating from the southern midlatitudes toward lower latitudes (Torsvik and Cocks 2016a), moving through different climatic zones. Active volcanism occurred during the whole Ordovician (Barnes 2004), accompanied by great siliciclastic input in adjacent areas (Servais et al. 2009), especially on the eastern side of Laurentia, where the Taconic orogeny had started during the Middle Ordovician (Holland and Patzkowsky 1996; Novack-Gottshall and Miller 2003). This greater siliciclastic input purportedly led to a greater genus richness on Laurentia (Miller 1997a). The increase in genus richness from the Dapingian to the Darriwilian common to Laurentia and Baltica corroborates previous observations restricted to specific organismal groups. For instance, Trubovitz and Stigall (2016) found a peak in brachiopod species richness on both Baltica and Laurentia in the middle Darriwilian by using range-through data compiled by Rasmussen et al. (2007). Similarly, Paris et al. (2004) recorded the onset of the GOBE in chitinozoans at the Dapingian/Darriwilian boundary on both Baltica and Laurentia, based on range-through data. It is noteworthy that their actual peak of Laurentian chitinozoans was found in the early Late Ordovician, while the Baltic peak was detected spanning from the Darriwilian to the early Sandbian (Paris et al. 2004). These findings are similar to our finding that peak genus richness occurred in the Katian Laurentia and Baltica (Fig. 3).

While the diversification dynamics on Laurentia and Baltica are broadly similar after the Dapingian/Darriwilian boundary, there are also differences between the two paleocontinents. At the Tremadocian/Floian boundary, there appear to be peaks in both origination and extinction rates on Baltica that are not expressed on Laurentia. Given the size of the 95% confidence intervals (see also Supplementary Fig. S7), we refrain from overinterpreting these observations. However, these peaks may indicate changes on Baltica before the Dapingian/Darriwilian peak that enhanced increased turnover, although this did not result in an increase in genus richness.

The abovementioned regional effects could have been dominating factors in the early Ordovician, while globally acting changes could have contributed to a coordinated increase of origination rates on Baltica and Laurentia during the GOBE (i.e., starting in the Dapingian and continuing into the Darriwilian). One such global change suggested as a driver of the GOBE is the proposed cooling that began during the Middle Ordovician (Trotter et al. 2008; Rasmussen et al. 2016) and was sustained until the Late Ordovician (e.g., Saltzman and Young 2005; Trotter et al. 2008). The suggested timing of the onset of cooling fits broadly with the increase in origination rates that we have found in our study. Frequent asteroid impacts during the Early and Middle Ordovician represent another worldwide condition that would have been experienced by the Ordovician seas (Schmitz et al. 2001), and these extraterrestrial inputs have also been linked to GOBE (Schmitz et al. 2008). This asteroid-driven GOBE hypothesis has recently been refuted on the basis of refined timings of these impacts using new isotopic data (Lindskog et al. 2017), but we argue that we cannot completely reject this hypothesis before we have temporal estimates of diversification that are also more highly resolved.

Ecological theory and contemporary empirical studies have shown that larger terrestrial areas can harbor more species than smaller areas (e.g., MacArthur and Wilson 1967; Rosenzweig 1995). However, this relationship is less clear for the marine realm (Roy et al. 1998; Valentine 2009), in part due to fluid boundaries. Despite uncertainties, we were still curious as to whether we could detect a greater total genus richness on Laurentia than Baltica (as found by Miller 1997a), given that circumstantial evidence indicates that Laurentia was much larger than Baltica (Torsvik and Cocks 2016a). Here, we assume, as in many paleontological studies (e.g., Sepkoski 1991; Finnegan and Droser 2005), that we can use genus richness as a proxy for species richness. The estimated genus richness for the Ordovician differs quite substantially between the two paleocontinents (Laurentia, 646 genera [95% CI 627, 683]; Baltica, 206 genera [95% CI 185, 294]), supporting our naive expectation. In addition, there is also temporal variation that is not coordinated through time (Fig. 3).

We restricted our paleocontinental analyses to genera that were observed on either of the paleocontinents but not both. In doing so in our analysis framework, we assumed that the genera in question actually originated on their assigned paleocontinents and that there are no unobserved occurrences on other paleocontinents prior to our first observations. We are aware that this assumption might be violated for some genera, but the Pradel seniority model does not consider observed or unobserved migration. In addition, we want to emphasize that while the Pradel seniority model assumes a closed population, the POPAN model assumes an open population. This may explain why the POPAN-estimated genus richness on Laurentia continuously increases substantially between Darriwilian and Katian, despite the Pradel-estimated net diversification rates on Laurentia only being slightly above zero.

Contrasting Onshore–Offshore Environments

In addition to examining paleocontinental differences, we also explored origination dynamics in two contrasting environments (onshore and offshore) in order to understand what may influence diversification rates. In addition to differing energy levels, sea-level changes may also more greatly affect onshore areas than offshore areas (Holland and Christie 2013).

Jablonski suggested that origination at the genus level should occur mostly where the higher-level taxa in question were already established (Jablonski et al. 1983; Jablonski 2005), while Tomašových et al. (2014) found evidence, using data from the Eocene and Plio-Pleistocene, that genus origination and extinction rates in onshore areas are greater than those in offshore areas. We note that Ordovician origination and extinction rates are only minutely different between on- and offshore environments (Fig. 4), results that are more consistent with Jablonski's hypotheses than with Tomašových's. Analyses based on model selection giving extra support for the similarity in origination rates between onshore and offshore areas are presented in the Supplementary Material.

In comparing onshore and offshore genus origination rates, we had to assume that any unobserved occurrences preceding the first record of the genus in question in a given environment are in the same environment as the first record. Given the incomplete temporal sampling of any genus, we realize this may not be a robust assumption, but we lack the tools to deal with this uncertainty at this time. However, given that we discarded genera with multiple paleoenvironmental associations in their first time interval of observation in the PBDB, we believe that “rogue” genera that occur in different environments before their first observation are rare. An additional caveat is that environmental classifications in the PBDB are not clearly defined, such that anyone entering data into the PBDB may have to make a subjective decision as to which environmental classification to assign to the fossil in question. As a result of uncertainty, many environmental assignments end up in the “indet.” (indeterminate) category. The onshore and offshore data set we used is quite a lot smaller than what was potentially available (see Supplementary Material). However, we have no reason to suspect that data points for which environment is marked “indet.” in the PBDB are biased toward any of the two environments in which we are interested. We are also aware that our division of onshore and offshore environments is more conservative than others published (e.g., Bottjer and Jablonski 1988) and that a different subset of data might suggest different results.

Conclusions

The global increase in genus richness at the Dapingian/Darriwilian boundary was the result of a dramatic increase in origination rates, while extinction rates stayed at relatively low levels. This global increase in origination rates is mirrored on Laurentia and Baltica, as well as in onshore and offshore areas, pointing toward a global mechanism behind the peak in origination rates at the Dapingian/Darriwilian boundary, which is likely the greatest contributor to the GOBE. Based on the patterns of paleocontinental and paleoenvironmental diversification dynamics, it appears that regionally varying factors were prominent during the Early Ordovician, while globally acting controlling factors appear to have been dominant during the Middle Ordovician and onward. All else being equal, our results suggest that origination probabilities are slightly greater in areas exposed to greater environmental change (i.e., greater disturbance through sea-level changes or continental movement).