Cascading trend of Early Paleozoic marine radiations paused by Late Ordovician extinctions.

The greatest relative changes in marine biodiversity accumulation occurred during the Early Paleozoic. The precision of temporal constraints on these changes is crude, hampering our understanding of their timing, duration, and links to causal mechanisms. We match fossil occurrence data to their lithostratigraphical ranges in the Paleobiology Database and correlate this inferred taxon range to a constructed set of biostratigraphically defined high-resolution time slices. In addition, we apply capture-recapture modeling approaches to calculate a biodiversity curve that also considers taphonomy and sampling biases with four times better resolution of previous estimates. Our method reveals a stepwise biodiversity increase with distinct Cambrian and Ordovician radiation events that are clearly separated by a 50-million-year-long period of slow biodiversity accumulation. The Ordovician Radiation is confined to a 15-million-year phase after which the Late Ordovician extinctions lowered generic richness and further delayed a biodiversity rebound by at least 35 million years. Based on a first-differences approach on potential abiotic drivers controlling richness, we find an overall correlation with oxygen levels, with temperature also exhibiting a coordinated trend once equatorial sea surface temperatures fell to present-day levels during the Middle Ordovician Darriwilian Age. Contrary to the traditional view of the Late Ordovician extinctions, our study suggests a protracted crisis interval linked to intense volcanism during the middle Late Ordovician Katian Age. As richness levels did not return to prior levels during the Silurian-a time of continental amalgamation-we further argue that plate tectonics exerted an overarching control on biodiversity accumulation.

Critical transitions in the state of the Earth that forced global ecosystems to adjust to changes in the physical environment fundamentally impacted Phanerozoic biodiversity levels (1). The Phanerozoic record of metazoan life denotes multiple events of abrupt change that affected the planet’s ability to sustain life. However, temporally highly resolved estimates of biodiversity change through geological time are lacking, making it difficult to assess when and how biotic or abiotic changes affected biodiversity in deep time.

The Early Paleozoic, and notably the Ordovician Period with its dramatic fluctuations in taxonomic richness, is a particularly relevant interval of the fossil record to test whether thresholds in the physical environment affect biodiversity accumulation. Previous richness curves through the Early Paleozoic show a modest increase during the earliest Cambrian, initiated by the Cambrian “Explosion,” followed by a more sustained radiation known as the Great Ordovician Biodiversification Event (GOBE) (2–6). However, the precise timing and duration of these events, or indeed whether or not they represent the same extended diversification pulse, are not sufficiently resolved as estimates of their timing differ by tens of millions of years (2, 6–8). In addition, there are diverging opinions regarding what happened after the GOBE: Did it mark the start of a Paleozoic biodiversity plateau only punctuated by one major mass extinction event during the Early Paleozoic (5, 6)? Or, did the diversification continue through the Silurian with only a minor disruption during the latest Ordovician (2, 3)? These inconsistencies in richness estimates originate due to different approaches to taxon counting. Whereas the classic work of Sepkoski (6, 9, 10) focused strictly on generic first and last appearances, Alroy et al. (2, 3) used genus-level occurrences to address sampling and preservation bias in a temporal binning framework of equal duration (Fig. 1). However, this latter approach came at the cost of the temporal resolution, which was then lowered from Sepkoski’s stage level binning (5 My) to an 11-My binning scheme. We complement these earlier works by providing a synthesis on biodiversity change during the first 120 My of the Phanerozoic partitioned into time slices with an average duration of 2.3 My.

Fig. 1.

Comparison of previous estimates on Early Paleozoic richness. The red shading [Sepkoski (9)] shows range interpolated presence/absence data partitioned into global stages; the green shading after Alroy (2) shows the sample standardized trend separated into 11-My time bins. Note that the y axis is arbitrary due to different estimates, with Ordovician peak in GOBE used to scale the two curves. Hence, this figure only shows the trends and relative timing of events.

We base our study on fossil occurrences compiled within the Paleobiology Database (PaleoDB), the data source that also served as the framework for the studies of Alroy et al. (2, 3). Entries in the PaleoDB consist of published taxon occurrences in specific strata and localities. Depending on the stratigraphic resolution of the published sources, it is possible to bin the occurrences in time intervals. Because time binning is not trivial and high-resolution chronostratigraphic bins are not available for the Early Paleozoic, these problems have been addressed previously by creating long ranging time slices in which generic ranges could be binned (2–4). To overcome these limitations, we established a set of 53 time slices through the entire Early Paleozoic based on biozones that can be correlated on a global scale by applying published chronostratigraphic schemes (11–15). We binned the PaleoDB genus occurrences within our time slices by using the lithostratigraphic information associated with each PaleoDB occurrence entry () (16).

Addressing sampling and preservation biases is fundamental. Therefore, we incorporated three different approaches to calculate richness estimates that address potential bias. Our primary approach is based on capture–recapture (CR) modeling (17–19). This method was conceptually transferred from ecology to fossil data following Nichols and Pollock (18) and Liow and Nichols (20), and allows for an explicit estimation of sampling probability (). We compared the results of the CR-modeling (DCR) approach with results from shareholders quorum subsampling, DSQS (2, 3), and from estimates of the Shannon entropy Hill number, DHill (17, 21). The results of the three approaches show a robust general trend, but the DSQS and the DHill estimates differ in having a larger overall volatility (). We interpret this volatility as an effect of uneven sampling (). Additionally, we compared our results partially against richness curves resulting from the alternative, automated time binning approach of Ordovician PaleoDB occurrences by Kröger and Lintulaakso (22). The two approaches yielded almost identical results ().

Results and Discussion

Our CR-richness curve (Fig. 2; ) shows a rapid rise in the early Cambrian followed by two falling phases corresponding to the Botomian and Marjuman extinctions. The Early Ordovician marks a slightly increased base level followed by a distinct Mid Ordovician rise in richness—the GOBE—peaking in the early Late Ordovician. This is followed by a prolonged fall in richness levels, corresponding to the end Ordovician mass extinctions. The early Silurian is characterized by a weak rebound, whereafter richness levels fall yet again later in the period.

Fig. 2.

Generic biodiversity history for the Early Paleozoic based on an average estimate after 50 model runs using the CR method. The temporal resolution of the dataset is shown by the alternating gray/white bars corresponding to the 53 time slices partitioned on the basis of biostratigraphically defined boundaries and matched with lithostratigraphical information extracted from the PaleoDB. Confidence intervals envelope the blue richness line. Major biotic events are annotated. In addition, Early Paleozoic biodiversity accumulation is set in the context of potential abiotic determinants. The sea level and strontium trends overall mirror the richness curve, reflecting that plate tectonics is a first-order determinant on speciation. The main radiation during the GOBE correlates strongly to cooling climate and a rise in oxygen levels. However, these determinants do not seem to accelerate biodiversity accumulation during the Silurian. Last, note the apparent association between Large Igneous Provinces (LIPs), abrupt, positive 13C values, and extinctions phases, suggesting a volcanic component in these intervals. Strontium isotope data are based on the LOWESS fitted line (77); LIPs are based on refs. 78–82, with darker red bands corresponding to large volcanic eruptions (69, 83–85); climate proxies are based on refs. 33, 35, and 63; the 13C record is compiled from ref. 11, and sea level is based on refs. 30, 32, 72, and 86.

The Cambrian–Middle Ordovician Interval.

The temporal resolution applied herein reveals a first-order stepwise rise in richness through the Cambro–Ordovician periods with the Cambrian Explosion and the GOBE being two separate diversification events. As such, the Cambrian Explosion not only marks a phylum-level expansion in biodiversity, but also a substantial increase in generic richness during a short burst in the earliest Cambrian Epoch 2. This follows the global stage-level data on invertebrate fossils (4), but our estimate shows a more significant rise in generic richness that is further sustained for ∼50 My despite the intervention of significant mass extinction events. This Cambrian–Early Ordovician biodiversity plateau, however, is low in comparison with the GOBE, which according to the current study took place during a narrow phase spanning just 15 My. This is in line with Sepkoski’s original definition of the Ordovician Radiation, as well as Miller and Foote (23), who used sample standardized, regional series-level data of Laurentia to calibrate the timing of this event. Since those early studies, opinions on the onset and duration of the events during the Cambro–Ordovician periods became more blurred, with some studies suggesting a late Cambrian onset of the GOBE (2, 3, 8). In our view, this opinion of a late Cambrian start of the GOBE more possibly reflects the longer time binning of previous richness studies (2, 3) and a lack of consensus as to what “biodiversity” refers to within the concept of the GOBE. There is undoubtedly an increase of abundance (as in number of specimens) from the Early Ordovician onward (24), but that does not necessarily reflect an increase in richness. Although a few higher clades originate during this interval, such as certain groups within the primary producers (8), in our view this does not justify expanding the timing of the GOBE into the Cambrian.

A number of studies have found a first-order correlation of richness to outcrop area and number of formations (25, 26). More strata deposited and preserved during times of high sea level stands would yield higher biodiversity levels because more habitats would be available in the first place and at the same time more rock volume would be available for the fossil record [the so-called “common-cause hypothesis” (27, 28)]. With the Ordovician probably representing the greatest transgressed rock-volume preserved of the Phanerozoic (29) and the early Katian part of the Late Ordovician possibly representing a Phanerozoic sea level maximum (30, 31), we expected higher biodiversity accumulation at this time. We tested for a correlation of the first differences between our DCR and a modified compilation of published sea level curves (30, 32) and the number of formations extracted from the PaleoDB (). The analysis supports this expectation. The first differences in richness correlate with first differences in number of formations (r = 0.836; P < 0.001) and a smoothed five-bin averaged sea level curve (r = 0.616; P < 0.001) (see for details). This is in accordance with previous analyses and with the common-cause hypothesis (25–28). The richness increase is associated with an increase in rock volume, but at the same time, the most drastic increase in richness, observed during the Middle Ordovician, occurs against a background of major sea level fall (Fig. 2). Cooling climate, as expressed in temperature estimates (33, 34) and rising oxygen levels (35, 36) may have been key factors for the diversification during this time (see further discussion on abiotic drivers below).

The Late Ordovician–Silurian Interval.

We find the highest biodiversity levels of the entire Early Paleozoic to be reached during the Sandbian–early Katian interval, but hereafter a three-phased fall in richness levels starts that ended during the earliest Silurian Rhuddanian Age (time slice Rhud2). This 10- to 12-My-long extinction interval does not follow the traditional view of a swiftly operating, two-phased end Ordovician mass extinction event confined to the Hirnantian Age (37).

Based on our estimates, the main drop in richness occurs during the earliest to mid-Katian. This is in contrast to previous estimates (38) and challenges the traditional link that is made to greenhouse–icehouse–greenhouse shifts (37). Evidence for icehouse conditions during the Hirnantian is well documented and undisputed (38, 39), but the timing of the onset of the biodiversity decrease seems to be diachronous based on the best resolved datasets, which arguably come from brachiopods (40) and graptolites (41). These analyses show peak diversities at different points in the late Katian. Our overall net loss in genus richness suggests an earlier onset of the extinctions and, thus, a prolonged survival interval compared with the classic perception of an approximately 1-My-long crisis phase during the Hirnantian Age (37). This implies that the global late Katian climatic optimum known as the Boda Event (42), as well as the regionally well-developed Richmondian Invasion (43), should be included in the crisis interval of the end Ordovician mass extinctions. These events were characterized by faunal dispersal and immigration, notably toward Laurentia, but also with enhanced migration at much higher latitudes where invading taxa regionally increased richness (44, 45). This broadening of generic geographic ranges enabled immigrants, not new species, to fill niches that had become vacant regionally. Thus, overall global biodiversity accumulation started to fall during the Katian (46). Dispersal (and migration) continued extensively during the terminal Ordovician Hirnantian Age, notably with the shelly benthos being characterized by the coldwater Hirnantia and Dalmanitina faunas reaching lower latitudes (47).

As Silurian richness levels did not surpass the Late Ordovician levels, our estimate is consistent with the “Paleozoic Plateau” suggested by Sepkoski (6, 48). However, it does not support the marked richness increase in the middle–late part of the Silurian shown by Alroy (2). This discrepancy likely reflects the low temporal resolution of the Alroy curve where just two time slices constitute the bulk of the Silurian Period. Thus, the latest Silurian is binned together with the early Devonian that likely pulled late Silurian richness levels up in that study (). More surprising in our results, is the lack of a biodiversity rebound during the Silurian. Given the robustness and high confidence level of our approach, we regard this signal as being well supported. This suggests that the end Ordovician extinctions had a far longer lasting impact on global biodiversity than previously recognized. This is concurrent with a comparably low-level ecological reorganization during the early Silurian Llandovery Epoch (41, 49, 50). Previous studies have shown that the ecological disruption succeeding the end Ordovician extinction events was of lower order than, for instance, after the end-Permian and end-Cretaceous mass extinctions (51, 52). Recovery intervals have been conceptualized as short intervals of a few million years with a geologically rapid reorganization of faunas where many clades adapted to the changing environments, conquered newly vacated ecospace and, as a consequence, radiated into new niches (48). Here, we point out that the Silurian interval differs from these more drastic recovery intervals in being comparatively protracted. This is supported by clade-level studies on rhynchonelliform brachiopods that show a diachronous recovery phase with a strong paleogeographical and bathymetrical component (53, 54). Brachiopod faunas recovered much faster during the earliest Silurian on Laurentia (55) than seen globally (56). However, recovery in generic richness was, to a large extent, the result of taxa that moved into vacant niches in shallow-water, tropical settings, and then eventually evolved into new families (50, 57, 58).

The Role of Abiotic Drivers in Early Paleozoic Biodiversity Accumulation

The current study enables a better comparison with environmental covariates considered to control biodiversity accumulation through time (Fig. 2). We selected a number of potentially relevant covariates and conducted simple pairwise tests for correlation of first-differences estimates using Pearson’s r as correlation coefficient ().

On the longest term, plate tectonics are believed to control biodiversity accumulation through the changing cycles in breakup and amalgamation of continents (59). At times with many geographic entities, provinciality is high, leading to increased isolation of faunas, which then become more endemic. Such plate tectonic cycles have been demonstrated on the Phanerozoic scale to yield higher global biodiversity (59, 60), and this has also been invoked as causal mechanisms behind the GOBE (61, 62). We have tried to test this indirectly, by comparing our richness estimate to the secular trend in 87Sr/86Sr, but no dependencies can be detected. There are two explanations for this. First, the strontium trend operates over too long timescales to be reflected by the current highly resolved binning framework. A second possible explanation is that more than one determinant is acting in concert during certain intervals of the studied interval and thus blur dependencies at our level of observation. The temperature trend, for instance, is coordinated without being directly correlated to richness. In Fig. 2, a first interval shows a decreasing temperature curve associated with roughly constant richness levels until the earliest Darriwilian. The succeeding early Darriwilian interval shows a temperature decrease associated with a richness increase. Hereafter, in a third interval that lasted until the late Katian, temperature and richness run largely parallel. During a fourth interval with highly volatile temperature estimates during the Hirnantian–Rhuddanian, the richness decreased until the Aeronian, but from the mid-Telychian both curves are again roughly parallel. These coordinated changes suggest complex links between the two trends. This is further illustrated in Fig. 3, which compares first differences of atmospheric oxygen and temperature to richness in a time series analysis, where time slices with discrepancies in change points are shaded in red. We found no dependencies between temperature and richness (), but, markedly, there are coincidences between major change points and trends (33, 63). These are best shown when richness and temperature estimates are smoothed as in Fig. 3. During the Middle Ordovician Darriwilian Age, richness reaches a maximum increase in the same time slice as data on conodont clumped-isotope thermometry show equatorial sea surface temperatures to fall within present-day conditions (33). There is now considerable evidence for the onset of Middle Ordovician glaciations during the earliest Darriwilian (34, 45, 64–66), and data from Baltica further suggest the onset of a Quaternary-scale ice age with a corresponding sea level fall in the order of 150 m at this time (34). Thus, the causal link between the onset of the GOBE and cooling climate is further corroborated by the current CR estimate.

Fig. 3.

Dependencies of global marine genus richness on pO2 and temperature. (A and B) Correlation between atmospheric oxygen (35) and richness. (A) Moving averages of pO2 and richness; (B) smoothed first differences of moving averages (Pearson’s r = 0.647; P < 0.005). (C and D) Coordinated change between richness and tropical sea surface temperature estimates (33, 63). (C) Moving averages; (D) smoothed first differences of moving averages (Pearson’s r = 0.306; P = 0.146). Averages are calculated over five time bins. First differences are calculated with a time lag of 1. The red overlays show time slices with discrepancies in change points. See for details.

The Middle Ordovician also saw changes in atmospheric oxygen levels that started to rise slowly during the early Darriwilian (35). This modeled estimate on pO2 correlates well with our richness estimates (0.647; P < 0.005), which appear to be strongly associated with changes in richness through major parts of the latest Furongian–Llandovery interval (Fig. 2). The time series analysis of first differences between richness and atmospheric oxygen show an overall good correlation, with time slices containing discrepancies mostly being concentrated during the Early–Middle Ordovician (Fig. 3). This suggests that richness generally is dependent on pO2, but also that particularly the Middle Ordovician interval is better explained by temperature changes. The slow rise in pO2 seen in Fig. 2 may reflect a lag time for the oxygenation of the atmosphere to build up compared with the rapid radiation of marine metazoan life suggested by our CR estimate. Molybdenum isotopes have been applied as a proxy for marine oxygenation levels (36). They too suggest a rise in the oxygenation levels of the deep ocean during the Middle Ordovician, corroborating the trend seen in atmospheric oxygen. Likely the increase in oxygen concentrations benefitted metazoan life and continued to sustain increased biodiversity accumulation as sea levels rose (32), expanding the habitable ecospace on the drowned epicontinental platforms. This is characteristic for the later parts of the Ordovician Period (29, 67). During the Silurian, pO2 increased rapidly to concentrations higher than Ordovician levels without a corresponding rise in richness (Fig. 3). This may suggest that oxygen levels, as with temperature, have the greatest impact on richness within a certain interval of boundary conditions.

Our data also reveal perspectives on the end Ordovician mass extinctions as well. Traditionally, this event has been linked to icehouse conditions during the Hirnantian Age, but more and more evidence for pre-Hirnantian glaciations has emerged within the past decade (34, 64, 65, 68), suggesting a less steep temperature gradient toward the end of the Ordovician than previously thought (33). We do not observe a distinct and sudden richness drop but rather a long, three-phased interval characterized by successively reduced generic loss. This is a markedly different trend than previous estimates suggest (6), but further studies specifically focusing on an estimation of the turnover rates are needed to resolve exactly when the main taxonomic loss takes place. Our observation invokes other causal drivers for the event than climate. Reviewing the published literature on Large Igneous Provinces (LIPs), and increased volcanism, shows a temporal correlation both to the early Cambrian Botomian Extinctions, as well as the early Katian drop in richness. Early Katian volcanic events have previously been suggested as a driver for this extinction event (69), but the temporal correlation was lacking. With the protracted crisis interval of our CR estimate, a volcanic driver no longer can be overlooked as a potential candidate for a causal determinant. This is further supported by increased weathering as implied by the Sr record suggesting less provinciality and the rapidly fluctuating 13C curve for the Late Ordovician. In fact, our review indicates a tendency for all major extinction intervals up through the Early Paleozoic to be slightly predated by large positive excursions in the carbon record. Positive excursions in 13C are generally believed to reflect increased primary production, as for example has been argued for in the case of the positive Middle Ordovician “Middle Darriwilian Isotopic Carbon Excursion” (MDICE) (34). However, as the MDICE is characterized by a slow increase toward heavier values, this contrasts with the repeated, abrupt positive pulses seen associated with the extinction intervals (Fig. 2). These, instead, could reflect rapidly increased atmospheric carbon events, indicating a hitherto-overlooked volcanic component in some of these Early Paleozoic extinction events, if not all. Therefore, a Late Ordovician scenario could be that increased volcanism triggered the extinctions with extreme warming and cooling phases, stressing faunas, and forcing them to migrate—in some cases successfully—but ultimately favoring those best adapted to either climate or facies. The stressed faunas during the prolonged survival phase were further challenged by the changing paleogeographic setting, most notably the configuration of continents during the latest Ordovician–Silurian interval. Decreasing provincialism and habitat destruction caused by amalgamating continents (70) further hindered a fast rebound of biodiversity accumulation during the earliest Silurian (56).

Conclusions

Our primary focus has been to produce a thoroughly tested Early Paleozoic richness estimate in high temporal resolution. In our view, this is needed to leverage the understanding of the effects of potential abiotic causal drivers in biodiversity accumulation through time. The associated intrinsic aspects of this accumulation, such as the turnover rates and the evolution of certain morphological features, are beyond the scope of the current study. The above discussion must, therefore, be carefully vetted and tested against such potential factors by future studies.

While controls on biodiversity change have been long discussed (71), better temporal precision on the chain of eco-evolutionary events and dating of the richness fluctuations themselves may resolve this debate. The two stepwise pulses of biodiversity increase observed in our curve echoes the classic ideas of the kinetic model of taxonomic richness suggesting that an equilibrium was reached twice during the Early Paleozoic (10). With our approach, the results of Sepkoski’s analyses are revived. This paves the way for future discussions as it suggests that relatively sudden shifts in environmental conditions, such as temperature and pO2, are critical in facilitating increased biodiversity accumulation. The Middle Ordovician Darriwilian Age may represent one such example of a fundamental change in the state of the Earth as global temperatures dropped while oxygenation levels increased in the marine and atmospheric realms. This appears to have had a dramatic impact on the accumulation of biodiversity seen during the GOBE, but the underlying evolutionary rates of the changes in richness still need to be resolved. The possible link between continental amalgamation and the end Ordovician biodiversity loss deviates from the traditional view of climate-induced end Ordovician mass extinctions. This alone indicates that favorable environmental conditions are not the only prime factor controlling biodiversity accumulation through time.

Methods

Our analysis is based on a download of all Cambrian–Devonian (Lochkovian) fossil occurrences registered in the PaleoDB (downloaded September 14, 2017) resulting in a total number of 199,796 occurrences. We divided the Cambrian–Silurian into 53 time bins plus 1 for the Lochkovian to account for range interpolation during the latest Silurian Pridoli Epoch (). The absolute (numeric) geochronological time spans of the individual bins are adopted from Ogg et al. (11), which also formed the primary basis for the global biozonal correlations conducted herein. Exceptions from this are during the lowermost Cambrian where we follow the biozonation used by Peng et al. (72), and for the Ordovician, where we used a modified time slice partitioning adopted from Bergström et al. (13) for the Lower–Middle Ordovician, and Ferretti et al. (14), in conjunction with Ogg et al. (11), for the Upper Ordovician. For the Silurian Period, we followed Cramer et al. (12) and Ogg et al. (11).

All lower Paleozoic formational names were extracted together with the fossil occurrences from the PaleoDB (September 14, 2017). Initially, all redundancies in formational names were removed. Hereafter, we partitioned the formational data into Cambrian (640 entries), Ordovician (1,650 entries), and Silurian (1,533 entries) datasets (Datasets S1–S4). Formations were assigned to our time slice partitions using the two latest renditions of The Geologic Time Scale (11, 15, 72, 73), as well as relevant literature covering regional stage names, and chronostratigraphic partitions no longer in use.

We manually time binned all PaleoDB lithological formation data (Datasets S1–S4) and filtered the occurrences for a resolution of two time bins or higher. Additionally, we used a PaleoDB download of all Ordovician fossil occurrences and a download of all Ordovician stratigraphic units from the RNames database (rnames.luomus.fi/) for an alternative automated time binning approach using the Ordovician time bins of and filtered the resulting list of occurrences for a resolution of two time bins or higher [see Kröger and Lintulaakso (22) for details on methods]. Both downloads are from September 14, 2017.

Subsequently, we filtered the two resulting datasets for taxonomic determinations at genus level and higher. This accrued a total of 185,348 occurrences for the manually binned dataset and 188,850 for the combined automated and manually binned datasets for the sample coverage methods and 25,298 and 25,597 occurrences, respectively, for the CR approach (see below). Based on these two different datasets we calculated richness trends using three different methods ():

The shareholder quorum subsampling approach (2), herein DSQS, using John Alroy’s R function, version 3.3 (bio.mq.edu.au/∼jalroy/SQS-3-3.R), was used.

The Shannon entropy Hill number (17, 21), herein DHill, was calculated using the R Package iNext, version 2.0.12 (74).

Additionally, the CR approach, herein DCR, was used for richness estimation. The method was transferred from ecology data to fossil data following the approach of Liow and Nichols (20), assuming that each genus is equivalent to a captured and recaptured organism, and that the total genus number is equivalent to the size of the population. A presence–absence matrix was constructed based on the filtered and binned PaleoDB genus occurrences for the time bins. This matrix served for the fitting of explicit models for richness estimation with time-varying probabilities of survival, sampling/preservation, and origination. We fitted the Jolly-Seber model following the POPAN formulation, also known as the “superpopulation approach” (75) (herein DCR). The CR estimates have been calculated using the program MARK (www.phidot.org/software/mark/) and the R Package RMark, version 2.2 (76).

The first two approaches are very similar, because they both compare samples of equal completeness and not equal sample size, such as, for example, classical rarefaction and they are based on the concept of sample coverage (see ref. 21 for a review). These methods differ in that DHill uses a unified sampling framework that seamlessly links rarefaction and extrapolation models depending on sample coverage, whereas DSQS uses rarefaction, exclusively.

The CR approach differs radically from the other methods in fitting explicit models for each dataset and extrapolating time specific probability-based diversities.

All data produced and analyzed by this study are available in , and the code, as well as associated R files, are available for download at https://doi.org/10.5281/zenodo.2586976.