Diversity dynamics of microfossils from the Cretaceous to the Neogene show mixed responses to events.

Microfossils have a ubiquitous and well-studied fossil record with temporally and spatially fluctuating diversity, but how this arises and how major events affect speciation and extinction is uncertain. We present one of the first applications of PyRate to a micropalaeontological global occurrence dataset, reconstructing diversification rates within a Bayesian framework from the Mesozoic to the Neogene in four microfossil groups: planktic foraminiferans, calcareous nannofossils, radiolarians and diatoms. Calcareous and siliceous groups demonstrate opposed but inconsistent responses in diversification. Radiolarian origination increases from c. 104 Ma, maintaining high rates into the Cenozoic. Calcareous microfossil diversification rates significantly declines across the Cretaceous-Palaeogene boundary, while rates in siliceous microfossil groups remain stable until the Paleocene-Eocene transition. Diversification rates in the Cenozoic are largely stable in calcareous groups, whereas the Palaeogene is a turbulent time for diatoms. Diversification fluctuations are driven by climate change and fluctuations in sea surface temperatures, leading to different responses in the groups generating calcareous or siliceous microfossils. Extinctions are apparently induced by changes in anoxia, acidification and stratification; speciation tends to be associated with upwelling, productivity and ocean circulation. These results invite further micropalaeontological quantitative analysis and study of the effects of major transitions in the fossil record. Despite extensive occurrence data, regional diversification events were not recovered; neither were some global events. These unexpected results show the need to consider multiple spatiotemporal levels of diversity and diversification analyses and imply that occurrence datasets of different clades may be more appropriate for testing some hypotheses than others.

Preserved remains of marine plankton and algae are highly abundant, cosmopolitan, and well represented in the geological record (Bé & Hutson 1977; Norris 2000). These marine micro‐organisms are inextricably linked to their environment and their species and genera richness fluctuates through time (e.g. O’Dogherty & Guex 2002; Erba 2004; Lazarus et al. 2014; Fraass et al. 2015). The exceptional preservation of microfossils provides a unique opportunity to investigate diversity dynamics, selectivity of extinction events, and sensitivity of certain ecological groups to different global‐scale environmental forcing events (Bown et al. 2004; Ezard et al. 2011; Lowery et al. 2020), including fluctuating atmospheric CO2 concentrations, sea surface temperatures (SSTs), and ocean circulation patterns. Additionally, biological interactions, such as symbiotic relationships, are catalysts for the evolutionary changes in major microfossil groups (e.g. photosymbiotic species (Norris 1996; Birch et al. 2012)). However, these must be assessed against impacts upon microfossil preservation, abundance, diversity and assemblages (Ezard et al. 2011; Pälike et al. 2012; Campbell et al. 2018). The improved temporal resolution of all major marine microfossil records has allowed investigation of the impacts of macroevolutionary events upon species‐level diversity (Lowery et al. 2020).

Numerous extinction events and major perturbations in climate have punctuated Earth’s geological history, varying spatially and temporally in intensity, tempo and selectivity, evidenced in the sedimentological record (Raup & Sepkoski 1982; Erwin 2001; Foote 2010). From the Cretaceous and throughout the Cenozoic, a myriad of different events occurred, each with different drivers and unique impacts upon and within different microfossil groups. The Cretaceous–Palaeogene boundary (K/Pg) is an example of one of the most significant and instantaneous extinctions triggered by a bolide impact on the Yucatán carbonate platform at c. 66 Ma (Alvarez et al. 1980; Schulte et al. 2010; Dinarès‐Turell et al. 2014). This global‐scale event instigated prolonged darkness, global cooling, and severe disruption to the ocean’s biological and chemical processes (Coxall et al. 2006) and subsequently impacted significantly on diversity and survivability at all trophic levels. Examples of large‐scale climate changes have also impacted significantly on microfossil diversity and assemblage composition, such as the global temperature increase of 5–8°C at the Paleocene–Eocene transition (c. 56 Ma) (McInerney & Wing 2011). The origination of the event is debated, yet common hypotheses suggest volcanism or extensive dissociation of methane hydrates triggering rapid short‐term warming, giving rise to the Paleocene–Eocene Thermal Maximum (PETM; Dickens et al. 1995; Petrizzo 2007; Andrade et al. 2018). This event subsequently affected microfossil diversity, dissolution and preservation, with many groups experiencing large shifts in geographical range, rapid evolution and changes in trophic ecology (Scheibner et al. 2005; Mutterlose et al. 2007; Jiang et al. 2018). Today, the current trajectories of contemporary climate change show that anthropogenic influence will continue to alter our complex marine ecosystems, posing an imminent threat of transitioning into the next global extinction (Ceballos et al. 2015; IPCC In Press). Rapidly increasing carbon dioxide levels and subsequent atmospheric temperature increase is just one of the ways that humanity is putting marine ecosystems under persistent stress (IPCC In Press). Despite current climate change having some unique characters relative to the geological record, the resulting successive biodiversity losses necessitate an investigation of the impact of past extinction events on microfossil diversity and assemblages to allow for predictions of how different marine micro‐organism groups may respond and recover today (Erwin 2001; McInerney & Wing 2011).

Microfossil diversity studies frequently restrict analyses to the first and last occurrences, and the ranges established by those first and last occurrences of taxa (Ezard et al. 2011; Fraass et al. 2015). These methods ignore key components that are integral to the cause and maintenance of diversity: the complex interactions of speciation and extinction (Silvestro et al. 2014, 2019). Relying solely on the raw fossil occurrence data poses problems with incompleteness and underrepresentation of many lineages, inconsistent sampling between clades (Foote et al. 2008), and it is unlikely to represent the true lifespan of a species (Alroy 2010). Several microfossil studies have tried to improve upon these methods by using subsampling and classic rarefaction curves (Rabosky & Sorhannus 2009; Lloyd et al. 2012, ; Renaudie & Lazarus 2013; Kocsis et al. 2014; Lazarus et al. 2014) to standardize samples and account for heterogeneity in sampling (Close et al. 2018). This, however, promotes misleading interpretations of occurrence data, compresses species richness ratios, and flattens diversity trajectories (Cárcer et al. 2011; Chao & Jost 2012). Recent methodological advances in phylogenetic methods provide detailed assessments of evolutionary dynamics through time in many organisms. Planktic Foraminifera have had their evolutionary relationships comprehensively examined throughout the Cenozoic (Aze et al. 2011; Ezard et al. 2011; Fordham et al. 2018), providing insights into diversification and ancestry. This method is often limited in power when estimating diversity as it is dependent on accurate phylogenetic tree estimation of evolutionary rates and times (Warnock et al. 2015; Stadler et al. 2018; Mitchell et al. 2019; Magee & Höhna 2021), which can be exacerbated when extant species are rare and extinction rates are high (Silvestro et al. 2014, 2018, 2019). Extinct and extant data have been incorporated into stochastic birth–death models to improve estimations of diversity and diversification rates (Liow et al. 2010), however, few integrate sampling and fossilization processes within diversity estimates, promoting significant biases (Etienne & Apol 2009).

Microfossil groups have different preservation potentials, largely related to test composition. Preservation of planktic foraminiferans and calcareous nannofossils varies through time with the positioning of the carbonate compensation depth (CCD; Van Andel 1975; Pälike et al. 2012). The CCD is the depth at which carbonate delivered to the bottom of the ocean is matched by dissolution, as stability of calcium carbonate declines with increasing depth, resulting in sediments below this threshold devoid of carbonate (Ridgwell & Zeebe 2005; Greene et al. 2019). In contrast, radiolarians and diatoms build with silica, forming ornate skeletons and frustrules (Kling & Boltovskoy 1999). Biogenic silica formation depends on the balance between organic productivity and other controls such as dissolution, upwelling and accumulation rates (Egge & Aksnes 1992; Kidder & Tomescu 2016). Undersaturation of silica can initiate the dissolution of biogenic silica in the water column and diagenesis of siliceous marine sediments (Van Cappellen & Qiu 1997). Thus, radiolarians and diatom oozes have a reduced preservation potential compared to calcareous oozes and are often under sampled (D’Hondt 2005; Dutkiewicz et al. 2016).

Development of methods to quantify diversification rates has identified the need to reconstruct preservation rates as a key input. The recently developed software PyRate jointly determines lineage duration and diversification dynamics by teasing out speciation and extinction rates independently, accounting for fossilization and sampling processes from occurrence data, while integrating sources of uncertainty within the analysis to limit bias in results (Silvestro et al. 2014, 2019). Using all possible data often neglected in other methods, including extant species (Silvestro et al. 2014), PyRate implements a robust hierarchical Bayesian framework to determine the most appropriate birth–death model for the data to estimate diversification rate changes through time (Silvestro et al. 2018). Using comprehensive microfossil occurrence data from the Neptune database (Renaudie et al. 2020), we reconstruct changes in diversification rates through 140 myr of major transitions to provide insights into the primary drivers of ecological and evolutionary change within microfossils from the Mesozoic to the Pliocene.

Material and method

Data source

Foraminiferan data were retrieved on 20 February 2020, while other fossil group data were downloaded on 26 May 2020 from Neptune Sandbox Berlin (NSB), the latest development of the Neptune database (Renaudie et al. 2020). NSB is openly accessible but requires interested scientists to request a username and password. The dataset holds information of over one million species occurrences of several microfossil groups pulled from scientific ocean drilling expeditions (Spencer‐Cervato 1999): the Deep‐Sea Drilling Project, Ocean Drilling Program, the Integrated Ocean Drilling Program, and International Ocean Discovery Program. Four fossil groups were downloaded: planktic foraminiferans, calcareous nannofossils, diatoms and radiolarians using the package ‘NSB Companion’ v2.1 (Renaudie et al. 2020) within R v4.0.2 (R Core Team 2013). Downloaded species occurrence data spanned c. 150–0 Ma, containing corresponding sampling details and palaeogeographical information (see Data S1). Microfossil datasets used had different time ranges; occurrence records began at 150 Ma for calcareous nannofossils and radiolarians, planktic foraminiferan records commence from 120 Ma, and diatom records start from 60 Ma.

Data manipulation

Using commands within ‘NSB Companion’, synonymous species were corrected in the dataset, removing non‐valid taxa. A dating uncertainty of ±0.25 myr was applied to occurrence sample ages, taken from standard deviations of NSB data (Spencer‐Cervato et al. 1994) and subsequent improved stratigraphy. The packages devtools v2.3.1 (Wickham et al. 2020), RCurl v1.98‐1.2 (Lang 2020), and packages within the tidyverse v1.3.0 (Wickham et al. 2019) were used to format the data into a table of taxon‐occurrence range entries. Utility functions provided alongside PyRate were used to extract single point ages for each occurrence sampled from a uniform distribution, bounded by the occurrence ranges that were input into the main PyRate program v3.0 (Silvestro et al. 2014). This extraction procedure was replicated ten times to account for the range of uncertainty associated with each occurrence (see Data S2). Despite downloading data from 150–0 Ma, only the data between 140 and 5 Ma were analysed due to large uncertainty margins in both the earliest occurrences and in extant data. This uncertainty is associated with truncated maximum and minimum ages in the dataset, which produced erroneous results.

Diversification modelling

PyRate estimates preservation rates using occurrences to infer speciation and extinction ages of individual taxa within a hierarchical Bayesian framework, from which rates of speciation and extinction through time were reconstructed (Silvestro et al. 2014). Incorporating all fossil occurrences removes the artificial modifications associated with subsampling methods. PyRate implements three different preservation models (the default non‐homogeneous, homogeneous, and time‐dependent Poisson processes) that vary the allowed heterogeneity of the preservation rates through time (Table 1). A gamma‐distributed parameter was also incorporated into these models to account for heterogeneous preservation between species. Each of the four major groups was run in a separate analysis, while compositional effects on preservation rates were assessed by grouping calcareous microfossils (foraminiferans and nannofossils) and siliceous microfossils (radiolarians and diatoms) into two further analyses. PyRate analyses for each taxonomic group were run as array jobs of the ten replicates generated on the Blue Pebble cluster at the University of Bristol (Advanced Computing Research Centre; see Table 1) using the recommended reversible jump Markov Chain Monte Carlo algorithm, generating 107 iterations, and sampling every 5000 iterations. Convergence was assessed in Tracer v1.7 (Rambaut et al. 2018) using the effective sample size (ESS) to ensure that the posterior distribution was suitably represented. Although not all model runs completed the iterations within the 72‐hour wall time, all model runs converged (ESS > 200) providing representative values for preservation rates and species longevity of each species.

Table 1

PyRate preservation models used to model speciation and extinction rates in all fossil groups, with the specifications and modifications made to find the best model for the data.

Model	Modifications
NHPP: Non‐homogeneous Poisson process of preservation Preservation rates change following a normal (bell‐shaped) distribution, for the duration of the lineage lifespan.	‐mG: to test the Gamma model of rate heterogeneity on all models. Accounts for preservation heterogeneity. The shape of the distribution is estimated from the data. Runs were completed with and without Gamma. After model testing, the model deemed the best would always be run with Gamma included.
	Sampling age errors to test for model sensitivity; all of the following ages (my) were run: ±0.25, ±0.20, ±0.15, ±0.10, ±0.05 myr.
	‐A 4: Reversible Jump Markov Chain Monte Carlo (RJMCMC): used on all other models as it is deemed the most accurate (Silvestro et al. 2014).
	‐edgeshift: when calculating rate shifts, the maximum and minimum boundaries of occurrences can cause apparent rate shifts in the data that reflect a potential sampling bias. We used ‐edgeShift to calculate the times of rate shift within fixed boundaries of 140–5 Ma to avoid erroneous results.
HPP: Homogeneous Poisson process Preservation rate is static throughout time and across lineages.	‐mG: see NHPP ‐edgeShift: see NHPP.
TPP: Time‐variable Poisson process A model that assumes preservation rates will change across certain time intervals but remain constant within timeframes; it tests whether preservation alters as a function of time rather than among lineages.	‐mG: see NHPP ‐edgeShift: see NHPP ‐pP: the preservation rate distribution shape and rate were changed from the default of ‐pP 1.5 1.5 to ‐pP 1.5 0; setting the rate to 0 allowed a rate parameter to be estimated from the data to reduce any subjectivity; the results showed no significant change to default and the default settings were used in subsequent TPP models. ‐qShift: the time intervals tested were changed; both epochs and ages of the geological record were run in the model. Using ages as pre‐defined timeframes produced the higher quality results with a reduced error margin and was used in all TPP tests on all microfossil groups.

Results of model selection are shown in Table 2.

Table 2

Taxon‐specific selection of the preservation model for each PyRate analysis using maximum likelihood.

PyRate Run	HPP		NHPP		TPP
	AICc	δAIC	AICc	δAIC	AICc	δAIC
Calcareous microfossil groups	1266659.0	43557.3	1223101.2 *	–	1235929.0	12827.4
Siliceous microfossil groups	610008.2	18433.2	591574.9 *	–	594251.4	2676.5
Foraminifera	349132.1	9439.3	339692.8 *	–	345945.1	6252.3
Nannofossils	908162.4	34773.7	873388.7 *	–	887865.2	14476.5
Radiolarians	425443.9	17520.6	407923.3 *	–	410398.7	2475.4
Diatoms	178146.1	11608.4	177212.4	10674.7	166537.7 *	–

Corrected Akaike Information Criterion (AICc) represents the AIC adjusted for small sample size. δAIC represents the difference in AIC from the statistically best model (highlighted in bold). Significance of the best model is indicated.

Significance at p < 0.01.

Model selection

Model selection for each microfossil group used a maximum‐likelihood framework to allocate the best fitting preservation. The fit of each model was calculated in the Akaike Information Criterion (AIC), providing a statistical method for assessing the fit of qualitatively different preservation models. Our main analyses used a standard precision of ±0.25 myr for each occurrence age, however, the precision of the Neptune data varies and is more precise than this. The sensitivity of the analyses to occurrence dating precision was investigated by running four additional PyRate analyses with decreasing occurrence ranges (±0.20, ±0.15, ±0.10, ±0.05 myr) using default preservation model settings (NHPP; Table 1). The results show no significant difference, suggesting dating uncertainty at this level does not impact upon the results (see Fig. S1).

The parameters estimated using PyRate that are presented in the Results section below are speciation rate, extinction rate and subsequent net diversification rates, as well as species longevity. These diversification rates are the number of speciation/extinction events per lineage per unit time (here per 1 myr, shown as ‘lineage−1Ma−1’; Silvestro et al. 2014). For net diversification rates, positive values represent speciation whereas negative values represent extinction. Where appropriate we also include the 95% credibility interval (95% CI). Note that in the birth–death process, mean species longevity equals the inverse of the extinction rate. These diversification rates are calculated for the four individual microfossil groups, and calcareous and siliceous groups. The significance of rate changes is calculated using standard log‐Bayes factors (BF) thresholds that determine statistical support (Kass & Raftery 1995). A BF of 2 represents a positive signal while a BF of 6 represents a strong signal.

Results

Sampling and preservation

Preservation rates tend to increase towards the present day (Fig. 1), however, there are notable deviations in calcareous groups. Foraminifera have relatively high preservation rates consistently throughout the Cretaceous with the lowest preservation rates in the Cenozoic at 60–20 Ma, while nannofossils have intermittent peaks of preservation from 150–60 Ma before steadily increasing towards the Recent. Sampling frequency of all four fossil groups improves towards the present day (Fig. 2). Radiolarians have extremely low sampling rates in the Mesozoic, prior to the occurrence of diatoms in the NSB. Nannofossils are the most sampled group, notably over 80 000 occurrences recorded in the Pliocene and Pleistocene. Nannofossil observations from the Cretaceous to the Recent number 280 934, followed by foraminiferans (157 958), radiolarians (116 527) and diatoms (114 845).

Fig. 1

Preservation rates of each of the four microfossil groups from 150 Ma to the Neogene (occurrences lineage−1Ma−1). Preservation in all groups increases towards the Recent. Calcareous groups generally have higher preservation rates than siliceous groups. Solid black line is the preservation rate and grey shaded regions represent the associated uncertainty. Vertical green bars represent the main events discussed in the text: VWAE, Valanginian–Weissert Anoxia Event; OAE‐1a, 1b, 1d, 2, Oceanic Anoxia Event 1a, 1b, 1d, 2; K/Pg, Cretaceous–Palaeogene Extinction Event; PETM, Paleocene–Eocene Thermal Maximum; EECO, Early Eocene Climatic Optimum; MECO, Middle Eocene Climatic Optimum; EOT, Eocene–Oligocene Transition; OMT, Oligocene–Miocene Transition.

Fig. 2

The number of occurrences of each microfossil group from 150 Ma to the Neogene in the NSB database, indicating sampling frequency through time.

Non‐homogeneous Poisson process (NHPP) preservation models were preferred for all microfossil groups, except diatoms (for which a time‐variable Poisson process (TPP) model was preferred; Table 2). For each microfossil group, the results presented below are from the respective best fitting model (incorporating taxic preservation rate heterogeneity).

The Cretaceous and the K/Pg

Calcareous microfossils

Following an initial decline in nannofossils at c. 135 Ma, their speciation and extinction rates remain stable through much of the Early Cretaceous (Figs 3A, B, 4A). Foraminifera experience a substantial but abrupt pulse (c. 2 myr) of rapid speciation around 105 Ma (Fig. 4B), shortly after the earliest records in this dataset. A more prolonged period of significantly higher rates of extinction (0.043 lineage−1Ma−1; 95% CI: 0.017 to 0.055 lineage−1Ma−1) in calcareous microfossils (BF > 6; Fig. 3B) between 95–90 Ma is associated with an overall net negative diversification rate (−0.002 lineage−1Ma−1; 95% CI: −0.010 to 0.005 lineage−1Ma−1) and a drastic decrease in species longevity (Fig. 3C, D). While net diversification rates return to prior levels, species longevity is reduced for the remainder of the Late Cretaceous until a brief spike in diversification before the K/Pg (Fig. 3), predominantly in the nannofossils (Fig. 4A).

Fig. 3

Calcareous microfossil diversification dynamics from 140–5 Ma. Plot showing the combined diversification rates (speciation/extinction events lineage−1Ma−1) of planktic foraminiferans and calcareous nannofossils. The solid rate lines represent the mean diversification rate, while the shaded area represent the 95% credibility interval (CI). The vertical green bars represent the events discussed in this study from the Cretaceous to the Neogene (see Fig. 1). A, speciation and extinction rates in calcareous microfossils. B, frequency of rate shift plot showing the significance of extinction and speciation rate change; PyRate is designed to generate a frequency of rate shift plot that corresponds to speciation and extinction plots; frequency of rate shifts are calculated using prior probabilities of rate shifts through time to compute posterior sampling frequencies; dashed lines indicate statistical support and are determined using standard log‐Bayes Factor (BF) thresholds; the lower line represents a BF = 2 while the top line represents a BF = 6, suggesting a positive or strong signal respectively (Kass & Raftery 1995). C, net diversification of the two fossil groups from 140–5 Ma; PyRate generates net diversification plots by subtracting extinction rates from speciation rates through time. D, mean species longevity of the two fossil groups.

Fig. 4

Extinction and speciation rates (left) of the four individual fossil groups with their corresponding frequency of rate shift plot (right) (speciation/extinction events lineage−1Ma−1). Solid rate lines represent the mean diversification rate; shaded area represents the 95% credibility interval (CI). Vertical green bars represent the events discussed (see Fig. 1). Dashed lines on the frequency of rate shift plots represent log‐Bayes Factor thresholds (see Fig. 3). A, calcareous nannofossil diversification dynamics. B, planktic foraminiferan diversification dynamics. C, radiolarian diversification dynamics. D, diatom diversification dynamics.

The largest decline in calcareous diversification in the whole sequence occurs at the K/Pg: extinction rates peak at 66 Ma (0.07 lineage−1Ma−1; 95% CI: 0.052–0.09 lineage−1Ma−1; Fig. 3A, B) and mean species longevity declines to its lowest overall level (c. 15 myr; Fig. 3D). Both foraminiferans and nannofossils are substantially affected (approximately 14‐fold and 5‐fold increases in mean extinction rates respectively) with very strong support for a change in rates (BF > 6; Figs 3A, B, 4A, B). Speciation rate increases across the boundary, simultaneously with increases in speciation rate within nannofossils (Fig. 4A).

Siliceous microfossils

There are no Mesozoic diatom data in NSB and consequently the siliceous groups are represented only by radiolarian data during the Cretaceous and K/Pg. Speciation rate is larger than extinction rate at the beginning of the record (140 Ma) at 0.12 lineage−1Ma−1 (95% CI: 0.06–0.23 lineage−1Ma−1; Fig. 5A, B) but gradually declines until a substantial decrease in mean speciation rate (0.10 lineage−1Ma−1) shifts at c. 125 Ma. Net diversification rate falls below 0.0 lineage−1Ma−1 (95% CI: −0.08 to +0.08 lineage−1Ma−1; Fig. 5C) and continues to decline, with a corresponding reduction in mean species longevity (c. 13 myr; Fig. 5D) by c. 120 Ma. Increasing extinction rates further reduce net diversification to −0.08 lineage−1Ma−1 (95% CI: −0.17 to −0.01 lineage−1Ma−1) by 115 Ma (Fig. 5A, B). In the early Late Cretaceous, speciation rates increase rapidly while, coupled with a corresponding gradual decrease in extinction rates, net diversification continues to increase throughout the Late Cretaceous to reach a peak at c. 76 Ma (Fig. 5A, B). Species longevity increases substantially in the latest Late Cretaceous to 40 Ma (Fig. 5D). Extinction rates remain constant across the K/Pg, however, speciation rates in radiolarians crash to about one‐third of prior levels (0.22 lineage−1Ma−1 to 0.07 lineage−1Ma−1; Fig. 4C).

Fig. 5

Siliceous microfossils diversification dynamics from 150–5 Ma, showing the combined diversification rates of radiolarians and diatoms (speciation/extinction events lineage−1Ma−1). Solid rate lines represent the mean diversification rate; shaded area represents the 95% credibility interval (CI). Vertical green bars represent the events discussed (see Fig. 1). A, speciation and extinction rates in siliceous microfossils. B, frequency of rate shift plot showing extinction and speciation rate change; dashed lines represent log‐Bayes Factor (BF) thresholds (see Fig. 3). C, net diversification of siliceous microfossils. D, mean species longevity of siliceous microfossils.

The Cenozoic

After the K/Pg, net diversification continues to decline until c. 63 Ma (Fig. 3C). This is concurrent with the origination of planktic foraminifer Praemurica inconstans yet remains negative until the onset of the Eocene. Extinction rates in nannofossils remain elevated for c. 16 myr at 0.045 lineage−1Ma−1 (95% CI: 0.030–0.055 lineage−1Ma−1) and remain higher than speciation rates until c. 40 Ma (Fig. 4A), whereas foraminiferal extinction rates are elevated for only c. 6 myr but at c. 0.12 lineage−1Ma−1 (95% CI: 0.09–0.16 lineage−1Ma−1) following the K/Pg (Fig. 4B). Mean species longevity of calcareous microfossils does not recover to its levels in the Late Cretaceous, increasing by c. 12 myr in the Palaeogene (Fig. 3D). Net diversification of calcareous groups is stable until c. 33 Ma when it increases above 0 due to a decline in extinction rate (Fig. 3C). Nannofossils and foraminiferans do not display significant extinction rate shifts in the Oligocene and Miocene (BF < 2; Fig. 4A, B).

At c. 64 Ma, mean speciation rates drop abruptly (Figs 5A, B) by 0.15 lineage−1Ma−1 and mean net diversification rate declines by 0.19 lineage−1Ma−1 (Fig. 5C). Species longevity erratically varies by over 30 myr between the later Paleocene and middle Eocene (Fig. 5D) associated with significant increases in the extinction rate in two pulses in diatoms (c. 46 Ma and c. 42 Ma, BF > 6; Figs 4D, 5A, B) that causes negative net speciation at 56 Ma (Fig. 5C). Extinction shifts are not acknowledged in the radiolarian record (Fig. 4C); however, a significant extinction rate shift occurs in the diatom record at 38 Ma (Fig. 4D). Simultaneously siliceous group species longevity declines to its lowest value of c. 6 myr (Fig. 5D). Prior to the Eocene–Oligocene transition (EOT), net diversification rates stabilize at around 0.0 (95% CI: −0.01 to +0.01 lineage−1Ma−1; Fig. 5C), which lasts until the Palaeogene–Neogene boundary when an increase in mean speciation rates of diatoms drives an increase in net diversification (0.08 lineage−1Ma−1 and 0.03 lineage−1Ma−1 increases in speciation and net diversification respectively; Figs 4D, 5C). An increase in extinction rate coincides with the Miocene–Pliocene transition in the diatom record (Fig. 4D).

Discussion

Cretaceous anoxia events

Global perturbations of the carbon cycle extensively affected the world’s oceans and their ecosystems. Diverse biotic and geological data provides evidence that Cretaceous oceans were extremely warm and high, and frequently experienced widespread oceanic anoxia events (OAEs) with important evolutionary impacts (Haq 1973; Leckie 1989; Wilson & Norris 2001; Leckie et al. 2002). Periodic and transient horizons of distinct carbon‐rich deposits in Cretaceous sediments are clear indicators of OAEs (Gambacorta et al. 2016), which indicate widespread changes in ocean geochemistry (Jenkyns 2010), sea levels, ocean circulation and vertical stratification (Leckie et al. 2002). These events caused widespread shifts in microfossil preservation and diversification dynamics (Lowery et al. 2020).

Valanginian–Weissert Anoxia Event (VWAE)

The calcareous nannofossil speciation decline at c. 135 Ma (Figs 3A–B, 4A) coincides with a positive carbon isotope excursion (CIE), detected during the late Valanginian (Weissert 1989; Lini et al. 1992; Wortmann & Weissert 2000). Widespread eutrophication (Bottini et al. 2018) associated with continental break‐up, basalt flooding and increased crustal production (Soua 2016) resulted in increased transport and nutrification triggering a ‘nannoconid decline’ at the onset of the OAE (Mattioli et al. 2014). Large error margins associated with the rate change in nannofossils at 135 Ma (Fig. 4A), may instead represent a false speciation rate shift due to truncated maximum and minimum sample ages in the dataset (Silvestro et al. 2014). However, nannofossil preservation reaches its highest value of the Cretaceous during this timeframe (Fig. 1) and therefore probably presents a representative shift in nannofossil diversification. Despite an overall decrease in calcareous diversification rates due to a ‘biocalcification crisis’ (Williams & Bralower 1995; Erba et al. 2004), nannofossil genera that indicate ocean fertilization such as Diazomatolithus display a marked increase at the VWAE (Barbu & Melinte‐Dobrinescu 2008). This possibly suggests that excess CO2 and warm surface waters resulted in higher fertility that promoted selectivity for genera that thrive in eutrophic conditions (Bersezio et al. 2002). Planktic Foraminifera are not present in the NSB data until c. 120 Ma since the records from the Berriasian to Barremian are scattered stratigraphically and geographically, making interpretations of their evolution difficult (Petrizzo et al. 2020). However, Silva & Sliter (1999) and Coccioni et al. (2007) noted the appearances of the first planispiral Globigerinelloides in the late Valanginian, suggesting that this OAE was not a driving force in foraminiferal evolution. Radiolarians experienced a biotic turnover during this event (Erba et al. 2004; Bottini et al. 2018), which is not reflected in this study (Fig. 4C), probably due to low preservation rates (Fig. 1) and low sampling frequency (<500 occurrences; Fig. 2) in the middle–late Cretaceous.

OAE‐1a

During the early Aptian (c. 120 Ma), OAE‐1a, a negative CIE event, caused an interval of organic‐rich black shale deposition in marine sediments from all major ocean basins (Robinson et al. 2004; Midtkandal et al. 2016; Bauer et al. 2019), leading to a decline in siliceous microfossil speciation and diversification rates (Figs 4C, 5A–C). The direct cause of OAE‐1a is still contested, however, increased carbon burial is probably the result of increased volcanism in the south‐west Pacific Ocean and an augmented greenhouse effect (Li et al. 2008). Subsequent changes in oceanic circulation patterns (Larson & Erba 1999), rising sea levels (Haq et al. 1988; Leckie et al. 2002), and enhanced weathering is likely to have caused alterations in biogeochemistry and primary production (Gomes et al. 2016). Under these conditions, deeper‐dwelling radiolarians experience species extinctions (Erbacher & Thurow 1997), probably as a result of reduced oxygen saturation or interruptions to vertical ocean stratification creating an unsuitable habitat. It is likely that this led to the last occurrences of deep‐water radiolarian genera Pantanellium and Podobursa at c. 120 Ma (Erbacher & Thurow 1997). Furthermore, a decline in sulphur concentrations and development of ferruginous conditions in the oceans are likely to have promoted anoxia, interrupting planktic micro‐organisms that rely on photosymbiosis (Bauer et al. 2019), suggesting that anoxia is significant in controlling siliceous diversity.

Nannofossil diversification does not fluctuate during this event, opposing an expected decline in calcification due to high pCO2 as seen in modern taxa (Erba & Tremolada 2004). However, it has been demonstrated that OAE‐1a only impacted certain groups (Lowery et al. 2020) and coupled with low sampling and preservation during this period (Figs 1, 2), it is probable that calcareous nannofossil decline would not show up in this study. In contrast, other foraminiferan studies demonstrate an increase in extinction rate (Leckie et al. 2002), which is not seen here (Fig. 4B). Foraminiferal occurrence data in NSB commences at 120 Ma with exceedingly low preservation rates (Figs 1, 2), hence OAE‐1a is unlikely to be captured in this analysis.

OAE‐1b

Characterized by cooling and sea‐level fall, OAE‐1b (late Aptian, c. 116 Ma) goes undetected in calcareous microfossils, despite it being recognized as the largest extinction within the Cretaceous for planktic Foraminifera (Fig. 3; Premoli Silva & Sliter 1999; Fraass et al. 2015; Lowery et al. 2020), and preservation is high during this period (Fig. 1). This may be the result of geographically and stratigraphically limited data during the Albian. Surviving planktic foraminiferans resided in the mixed layer and were typically microperforate species (Petrizzo et al. 2020), whereas deeper dwellers were predominantly impacted by reorganization of North Atlantic deep‐water circulation and break down of water column stratification (Huber & Leckie 2011). Consequently, limited Foraminifera sampling of deeper oceanic environments during this period (see Fig. S2) may have led to this extinction event going undetected in this analysis. OAE‐1b presents a dramatic reduction in foraminiferal test size during this time, with average sizes dropping by c. 100 µm (Huber & Leckie 2011). Associated with fundamental changes of foraminiferal ultrastructure, large and ornate taxa are shown to disappear in the late Aptian, replaced with simpler and smaller forms (Premoli Silva & Sliter 1999; Leckie et al. 2002). Assemblage composition change is associated with prolonged cooling and possible acidification from Kerguelen volcanism (Leckie et al. 2002; Petrizzo et al. 2020).

OAE‐1b is not reflected by a pronounced change in siliceous diversification since it is unlikely to have recovered from OAE‐1a due to the greenhouse nature of the ocean–climate system (Leckie et al. 2002). Net diversification of siliceous microfossils kept stable at c. –0.08 lineage−1Ma−1 (95% CI: −0.15 to −0.01 lineage−1Ma−1) for the remainder of the Aptian, with elevated extinction rates (Fig. 5A–C) spanning the duration of the mid‐Cretaceous. The poor siliceous preservation rate (Fig. 1) during the Cretaceous may have prevented key biotic events from being detected by PyRate. Additionally, OAE‐1b had a generally regional impact constrained to certain regions such as the North Atlantic and Mediterranean ocean basins, contrasting the global distribution observed at OAE‐1a (Premoli Silva & Sliter 1999).

OAE‐1d

The late Albian to early Cenomanian oceanic anoxic event (c. 104–99.5 Ma) represents a global interval of black shale deposition, distinguished by a moderate (0.5–2‰) positive carbon isotope (δ13C) excursion (Bralower et al. 1993; Leckie et al. 2002; Watkins et al. 2005; Meyers et al. 2006; Giorgioni et al. 2012). OAE‐1d has been proposed as an interval of global warming, sea‐level rise, and changes in water column stratification due to widespread oxygen deficiency (Watkins et al. 2005; Rodríguez‐Cuicas et al. 2019; Petrizzo et al. 2020), leading to an evolutionary turnover of planktic Foraminifera (Fraass et al. 2015). A significant radiation of planktic Foraminifera is observed at c. 105 Ma (Fig. 4B), associated with the origination of keeled, deep‐dwelling, and morphologically complex taxa due to increased ecospace (Hart 1990). Extinction associated with a turnover is not reflected in Foraminifera results (Fig. 4B), although occurrences of both Ticinella and Bitticinella, which developed during OAE‐1b, cease at c. 103 Ma (Fraass et al. 2015). New mixed layer genera originated, including Paracostellagerina and Planomalina (Petrizzo et al. 2008), suggesting they were speciating to occupy the ecological niche space, vacated by mixed‐layer genera which disappeared during OAE‐1b (Lowery et al. 2020). OAE‐1d is enigmatic and intermittent stratification interruptions strongly influenced anoxia in bottom waters, initiating a shift from calcareous to siliceous deposition (Gambacorta et al. 2016). Formation of a pycnocline (Tiraboschi et al. 2009) would have destabilized calcium carbonate and established an oxygen minimum zone, initiating a loss of deep‐sea habitats (Erbacher & Thurow 1997). Consequently, it is probable that limited occurrence data during this period, especially from pelagic ocean basins, resulted in an extinction event going undetected in Foraminifera.

Eutrophic conditions, rising sea levels, and the collapse of density gradients that developed during OAE‐1d probably led to a turnover in Radiolaria (Wilson & Norris 2001; Leckie et al. 2002). Radiolarians are at a net negative diversification rate of –0.02 lineage−1Ma−1 (95% CI: −0.18 to +0.2 lineage−1Ma−1) at 105 Ma, before steeply increasing towards the Early–Late Cretaceous boundary (Fig. 4C). Erbacher & Thurow (1997) described the disappearance of 15 radiolarian taxa and origination of 12 new taxa. The abiotic events at the Albian–Cenomanian boundary include an eustatic maximum, which may have allowed radiolarians to diversify precipitously due to a rapid expansion of new niches and increased temperatures (O’Dogherty & Guex 2002; Wang et al. 2019). This long‐term sea‐level rise, in addition to oxygen deficiency in bottom waters, is likely to have caused increased marine productivity, upwelling of nutrient‐rich waters, and controlled the shift from calcareous to siliceous deposition, making trends in sea level and eutrophic conditions key factors in radiolarian diversification.

OAE‐2

The Cenomanian–Turonian boundary (c. 94 Ma) marks one of the most widespread carbon cycle perturbations of the Phanerozoic Eon (Turgeon & Creaser 2008; Clarkson et al. 2018). Oceanic anoxia increased globally by at least a factor of three compared to the Recent (Montoya‐Pino et al. 2010); this, coupled with high atmospheric CO2 and rapid warming due to the eruption of the Caribbean large igneous province (LIP), affected the biota (Leckie et al. 2002; Lowery et al. 2020). Previous studies have generally shown a stepwise extinction in Radiolaria (Erbacher & Thurow 1997; Lowery et al. 2020); our results, however, show speciation rates in siliceous microfossils continuing to increase following OAE‐1d (Fig. 5A, B), and net diversification exceeding 0 lineage−1Ma−1 by c. 95 Ma (Fig. 5C). Warmer temperatures during OAE‐2 triggered rapid pulses of terrigenous nutrients weathered from LIPs (Kuroda et al. 2007), which may in turn have promoted an increase in oceanic silica saturation (Jenkyns 2010). During OAE‐2, climate changes were influential on deep circulation dynamics (Gambacorta et al. 2016). Warm and saline bottom waters triggering the vertical advection of nutrients and expansion of the oxygen minimum zone are likely to have created favourable conditions for the origination of new radiolarian taxa (Huber et al. 2002; Musavu‐Moussavou et al. 2007; Wang et al. 2019). The vertical distribution and ecological preferences of radiolarians are vitally important to the reported net diversification (Fig. 5C), since deeper dwellers tended to go extinct whereas those in surface waters and shallower depths speciated and increased in diversity (Wang et al. 2019). Anoxic bottom waters and warm conditions typically promote high preservation potential in siliceous sediments due to a lack of bioturbation to regenerate dissolved silica back into the water column (e.g. Kidder & Tomescu 2016). However, the number of occurrences (Fig. 2) in NSB within siliceous groups during this period may be responsible for an extinction event going undetected.

Calcareous microfossils depict an opposing trend with a significant increase in extinction rate commencing at c. 95 Ma (Fig. 3A–C). Breakdown of vertical density gradients during the later Cenomanian may have contributed to extinctions within deeper‐dwelling Foraminifera, such as Globigerinelloides bentonensis (Leckie et al. 2002). Evidence of declining calcareous nannofossils diversity across the Cenomanian–Turonian boundary (Lowery et al. 2020) also contributed to the decreased calcareous diversification rates. Calcareous species longevity simultaneously plummets by over 40 myr (Fig. 3D), probably influenced by the dominance of short‐lived, smaller, and opportunistic taxa able to tolerate low oxygen conditions (Coccioni & Galeotti 2003). High‐nutrient waters promoted the proliferation of the foraminiferal genus Schackoina; it is likely that the elongated and ornamented chambers were important to maximize oxygen and nutrient uptake (Georgescu 2012). At c. 90 Ma, diversity and mean species longevity begins to recover but they do not reach levels prior to the elevated extinction. This would suggest that the anoxic event had a long‐lasting impact on calcareous populations, with net diversification in decline for c. 5 myr. Gradual yet not complete recovery of foraminiferal assemblages is also noted by Coccioni & Galeotti (2003), resulting from a well‐developed oxygen minimum zone. Despite the anoxia event lasting only c. 440 000 years, habitat loss and mass extinction resulted in severe assemblage disruption (Pogge Von Strandmann et al. 2013) and deterioration of ecological space (Parente et al. 2008).

Cretaceous–Palaeogene extinction event

The largest and most highly significant extinction rate increase of the calcareous record occurs at c. 66 Ma with a precipitous decline in net diversification rate (Fig. 3A–C). The extinction in both terrestrial and marine realms immediately before and across the Cretaceous–Palaeogene boundary appears to have been driven primarily by the Chicxulub impact (Alvarez et al. 1980; Schulte et al. 2010; Vellekoop et al. 2016; Hull et al. 2020). An increase in extinction rate is seen to occur in both calcareous nannofossils and planktic foraminiferans (Fig. 4A, B) with the largest impact observed in planktic foraminiferans (c. 0.14 lineage−1Ma−1 increase), as c. 95% of taxa within planktic foraminiferans go extinct across the K/Pg caused by incredibly volatile and hostile conditions (Arenillas et al. 2016). Aerosols, soot from wildfires, and ejecta sparked a short‐impact winter (Wolbach et al. 1990; Tabor et al. 2020), reducing climate forcing and ceasing photosynthetic activity (Vellekoop et al. 2014). In turn, pelagic stratification broke down causing the mixed layer to plummet from its normal depth of c. 100–150 m to >2500 m (Brugger et al. 2017). The combined impact of these stressors promoted high extinction and low speciation rates in planktic foraminiferans. Nannofossils display both high extinction and high speciation rates (Fig. 4A), reflected in the calcareous record (Fig. 3A, B). Following the short‐term cooling, aerosols and vaporized carbonates resulted in severe warming (Kring 2007). Enhanced ocean mixing and additional nutrients from settling out ejecta, caused increased nutrient transport from the deep ocean to the surface, and promoted primary production (Brugger et al. 2017). Consequently, it is likely that aberrant, opportunistic ‘disaster taxa’ originated in regional blooms shortly after the impact in calcareous groups (Arenillas et al. 2018). The pathway the biota took to recover from the mass extinction, and precisely how to define ‘recovery’, is an active subject of current research (e.g. Hull et al. 2011; Schueth et al. 2015; Lowery et al. 2018). Corresponding with species longevity plummeting, these taxa would have been short‐lived, colonizing new morphospace, before diversity began to increase again, resulting in high turnover rates (Lowery & Fraass 2019). Net diversification in calcareous groups takes c. 3–4 myr to begin to recover (Fig. 3C). This is concurrent to the re‐establishment of planktic foraminifer photosymbiosis at c. 63.5 Ma (Norris 1996; Birch et al. 2012). The recorded recovery, however, does not appear to be globally synchronous (Hull & Norris 2011). Danian taxa tended to be unusually small (<150 µm; Gallala et al. 2009) since smaller, more angular forms provided a selective advantage for harbouring symbionts in the first photosymbiotic taxa (Birch et al. 2012). Deep‐sea transfer of organic matter via the biological pump was severely interrupted (Coxall et al. 2006) inhibiting deep‐sea carbonate sedimentation (Peters et al. 2013), but evidence of new originations suggests that it was less severely impacted and for less time than previously thought (Birch et al. 2012). Extinction remains high in nannofossils well into the Cenozoic, taking 10 myr to reach pre‐extinction diversity (Alvarez et al. 2019). This crisis interval is dominated by short‐lived ‘bloom’ taxa such as Futuyania and Prinsius that thrive in eutrophic conditions. These are replaced by oligotrophic taxa as the climate warms and global thermal stratification increases into the Paleocene (Hilting et al. 2008). Shallow, oligotrophic taxa were increasingly uncompetitive and the most successful in the Paleocene, despite experiencing some of the highest extinction rates at the boundary (Jiang et al. 2010), perhaps suggesting the simultaneous significant speciation and extinction rates (Fig. 4A). High extinction and speciation rates across the K/Pg in calcareous fossil groups demonstrate that time is required to rebuild ecological niche space to promote novel radiations and increasingly complex morphologies to promote diversification (Lowery & Fraass 2019), though the durations of these radiations are different in this analysis.

Radiolarians, on the other hand, demonstrate continually high speciation and diversification rates throughout the Late Cretaceous and initially across the K/Pg, while extinction rates remain low (Fig. 5A–C) with species longevity reaching its highest values across the boundary (c. 42 myr; Fig. 5D). Previously analysed marine sediments suggest moderately rich radiolarian assemblages, with both common and rare species discovered (Jianbing & Aitchison 2002; Hollis & Strong 2003). Data for radiolarians and other siliceous groups appears poorly represented across the K/Pg (Lowery et al. 2020) relative to other parts of their record, although preservation increases at c. 65 Ma (Fig. 1). At several sites noted by D’Hondt (2005), concentrations of silica increased at the K/Pg and remained high for at least 1–2 myr following the mass extinction, interpreted as an interval of high biosiliceous productivity (Hollis & Strong 2003). Nutrient‐rich waters and cool upwelling zones in the earliest Paleocene are likely to have advanced siliceous plankton speciation (Strong et al. 1995). This suggests that SSTs, nutrient input, and ocean circulation patterns are integral in radiations of radiolarians.

Paleocene and Eocene

Delayed recovery in calcareous taxa following the K/Pg parallels the prolonged extinction in nannofossils, which remain elevated for c. 16 myr (Fig. 4A). Net diversification of calcareous taxa gradually increases from c. 63 Ma (Fig. 3C) and reaches equilibrium just after the Paleocene–Eocene boundary (56 Ma). The Paleocene–Eocene transition is indicated by an abrupt onset of a negative CIE (Kennett & Stott 1991; Zachos et al. 2001) associated with the most significant climatic event of the Cenozoic (Paleocene–Eocene Thermal Maximum (PETM); Zachos et al. 2008; Jiang et al. 2018) and widely used as an analogue for anthropogenic warming (Hönisch et al. 2012; Zeebe & Zachos 2013). During the PETM, calcareous microfossils indicate significant rate changes in extinction and minor rate changes in speciation (Fig. 3B), probably a consequence of cold‐water taxon extinctions and speciation among warm‐dwelling taxa (Mutterlose et al. 2007). Nannofossil taxa considered proxies of increased SST (e.g. Discoaster araneus, Rhomboaster cuspis and Tribrachiatus bramlettei) originated between 56 and 55 Ma (Erba & Tremolada 2004; Mutterlose et al. 2007). Conversely, low‐latitude taxa were dramatically impacted by the PETM (e.g. losses within Fasciculithus). Dissolution linked to calcification depletion from increased CO2 and subsequent acidification caused cold water taxa to suffer (Raffi et al. 2005). Despite the ecological response at the PETM, Gibbs et al. (2006) noted that the rate of turnover was relatively modest and both cool‐water and warm‐water taxa appeared and disappeared making significant extinctions quite rare (Fig. 4A). As a result, the main impact of the PETM was the rapid change in evolutionary turnover as opposed to an individual environmental factor (Gibbs et al. 2006), promoting the high latitude extinctions of taxa living at their ecological limits. Planktic foraminiferan diversity increased overall (Fig. 4B; Ezard et al. 2011; Fraass et al. 2015; Lowery et al. 2020). Acidification at the PETM is of the same magnitude as the K/Pg extinction (Penman et al. 2014), yet diversification increases overall; indicating that acidification is an unlikely main control on calcareous extinctions (Ridgwell & Schmidt 2010; Hönisch et al. 2012), while SST fluctuations thus potentially impact calcareous diversity more. Species longevity decreases across the PETM (Fig. 3D), owing to the origination and subsequent extinctions of short‐lived excursion taxa such as Morozovella allisonensis (Kelly et al. 1998; Bralower & Self‐Trail 2016).

Following the initial decline after the K/Pg, siliceous extinction rates increase significantly (Fig. 5A–C), and species longevity declines (Fig. 5D) simultaneous to the onset of the PETM. Siliceous diversity throughout the Palaeogene is dynamic, probably resulting from several turnover intervals (Oreshkina 2012) corresponding to significant changes in diatoms (Fig. 4D). Previous studies have noted a turnover in siliceous groups with warm‐water taxa migrating poleward (Liu et al. 2011; Oreshkina & Radionova 2014), an extinction of cold‐water taxa, and a rapid proliferation of tropical siliceous taxa across the PETM (Mitlehner 1996). It is likely that these previously seen results are not detected here due to the weakness of the silicate record in this interval, reflected in low preservation and sampling rates (Figs 1, 2). Poor preservation is a known feature of silicate records studies (Oreshkina 2012), and a potential avenue for work moving forward.

After a c. 3 myr negative perturbation in siliceous taxa diversification, speciation begins to increase. The rise in biodiversity at c. 52–50 Ma has been connected to the Early Eocene Climatic Optimum (EECO; Oreshkina 2012). The EECO has been causally linked to high atmospheric CO2 levels thought to correspond to intensified volcanic emissions (Zachos et al. 2008), though the extent of these high concentrations is poorly constrained (e.g. Yapp 2004), and the impacts of EECO appear regionally disparate. Extreme changes in regional temperatures and changes in chemical weathering (Zachos et al. 2008), are likely to have promoted biosilica depositions in the Arctic (Barron et al. 2015). An alternative suggestion is that the Tasmanian Gateway started opening at c. 50 Ma, terminating the EECO and leading to an abrupt biosilica deposition (Bijl et al. 2013); supported by tectonic studies of seafloor spreading (Livermore et al. 2007). Following the EECO, increases in speciation rate occur in two stages in the diatom record with spikes at c. 46 Ma and c. 42 Ma. These speciation events coincide with the middle Eocene cooling trend and Middle Eocene Climatic Optimum (MECO) respectively (Barron et al. 2015). Middle Eocene cooling driven by obliquity cycles initiated extensive diatom deposition in the Atlantic Ocean (Vahlenkamp et al. 2018), reflected here (Figs 4D, 5C). Evidence suggests that Synedropsis, a genus present at c. 47–46 Ma, are sea‐ice dependent and regular melting episodes promoted silica‐rich water, hence increased speciation rates (Stickley & Koc 2009; Stickley et al. 2009). Transient warming at c. 40 Ma interrupting the long‐term Eocene cooling trend characterizes the MECO (Barron et al. 2015), coinciding with a spike in siliceous diversification and diatom speciation rate. Supporting the results of Witkowski et al. (2014), high diatom and radiolarian accumulation rates are observed, resulting from increased nutrient supply and subsequent eutrophication from elevated SSTs.

During the Eocene, global carbonate production was in decline, controlled by eccentricity‐scale changes resulting in sudden dissolution episodes (Rivero‐Cuesta et al. 2019). However, this was only recorded regionally (Arimoto et al. 2020) and consequently not picked up by these analyses. Bleaching was evident in several symbiont‐bearing planktic foraminiferans due to the warming (Edgar et al. 2013), yet net diversification rates remain constant within both calcareous fossil groups, despite a fluctuating Eocene climate and highly variable CCD (Fig. 4A, B; Lyle et al. 2005). Siliceous diversity peaked during the MECO warming which abruptly ended at c. 38–39 Ma as oxygenation recovered and SST decreased (Bohaty et al. 2009). Rapid reoxygenation is likely to have initiated a breakdown in water column stratification resulting from decreased fresh‐water input at the cessation of the MECO (Boscolo Galazzo et al. 2013), concurrent with a significant change in siliceous diversity. At 38 Ma, siliceous groups experience an abrupt yet short‐lived decline in net diversification (Fig. 5C) and species longevity (Fig. 5D). A decreased sea level and a widespread deep‐sea hiatus occurring at c. 38.5 Ma has been put forward as the explanation for a decline in diatom abundance (Wade & Kroon 2002). Discourse as to the trigger of this event is still ongoing. It has been proposed as the outcome of the opening of Drake Passage and Tasman Gateway (Egan et al. 2013) or the onset of the Northern Component Water (Borrelli et al. 2014), both initiating changes in ocean circulation patterns and influencing the diversity of siliceous taxa. This perturbation is temporally short and siliceous net diversification increases sharply after, paralleled with recorded increases in opal deposition due to increased upwelling (Diekmann et al. 2004).

Eocene–Oligocene Transition

The descent into an icehouse regime continued through the Eocene and a continental ice sheet developed at the Eocene–Oligocene Transition (EOT; Pearson et al. 2008). Ice sheet growth was triggered by orbital forcing controlling CO2 concentrations (Coxall & Pearson 2007), leading to SSTs dropping by 3–5°C (Galeotti et al. 2016). Hypothesized as a turnover or extinction event in all fossil groups, it was expected that EOT diversity would reflect a move to cold‐water taxa dominance while earlier Paleocene–Eocene taxa became extinct. The EOT is typically described as a large extinction in calcareous taxa, due to decreased pCO2 and subsequent cooling (Wade & Pearson 2008; Lloyd et al. 2012; Lowery et al. 2020), while siliceous taxa undergo a turnover, with diatoms gaining c. 30 new species (Lazarus et al. 2014; Barron et al. 2015). In contrast to these previous findings, we show a significant shift in extinction rate (Fig. 3B) resulting in a slight increase in diversification in calcareous taxa (Fig. 3C) and a decrease in diversification in siliceous groups (Fig. 5C), though the latter appears to start prior to the EOT. The CCD during the early Oligocene deepened (>1 km) within the tropical Pacific, roughly concurrent with ice sheet growth and sea‐level decline, as deep‐water deposition went from opal‐rich to carbonate‐rich sediments (Coxall et al. 2005), thus increasing carbonate preservation potential. These changes, however, are not reflected here.

The Eocene–Oligocene boundary is indicated by the rapid extinction of all remaining species within the Hantkeninidae (Wade & Pearson 2008; Houben et al. 2012), with the last occurrence in NSB of Hantkenina australis at 34.6 Ma. The Hantkeninidae extinctions were due to changes in marine stratification and pronounced productivity changes (Pearson et al. 2008). Despite evidence of small rate shifts within nannofossils (Fig. 4A), these do not exceed the BF thresholds, resulting in calcareous diversification remaining stable from the EOT until the end of the record at 5 Ma. Since radiolarians are well‐preserved and well‐sampled (Figs 1, 2) during the EOT, the decline in siliceous diversification rate may be the result of local radiolarian extinctions with 13 radiolarian taxa disappearing within a 100 000 year period in middle–low latitudes (Kamikuri & Wade 2012). Poor sampling in high latitudes has possibly resulted in speciation going undetected and net diversification subsequently declining.

Oligocene–Miocene transition

Palaeoclimate records suggest that congruence of eccentricity and obliquity cycles initiated a major transient glaciation at the Oligocene–Miocene Transition (OMT; c. 23 Ma; Zachos et al. 2001). However, recent research proposes that it may be the result of steadily decreasing CO2 levels like the glaciation at the EOT (Greenop et al. 2019). Despite its controversial onset, the glaciation during this period is likely to have led to sea‐level decline (Zou et al. 2016), and fluctuations of the CCD into the Miocene (Pälike et al. 2012). However, these changes do not correlate with increased carbonate species diversity here (Fig. 3C). Increased vertical mixing and intensified ocean circulation led to enhanced latitudinal temperature gradients, promoting an increase in nutrient supply and silicate saturation (Naish et al. 2001). Amplified ocean productivity and cool SSTs subsequently drove an increase in diatom species richness (Suto 2006), reflected in Fig. 4D. The important genus Chaeotoceros diversifies during this transition, a key development for modern oceans with more than 400 extant species (Suto 2006). Atlantic and Indian Ocean basins in the Miocene have been suggested to have a deep thermocline periodically interrupted by vigorous upwelling, promoting increased productivity at the surface (Smart & Thomas 2006). This upwelling, along with intermittent glacial blanketing of the Antarctic continent, interrupted chemical weathering feedbacks and caused a rise in pCO2. Warming of up to 2°C is likely to impact foraminiferal assemblages and could have posed a selective advantage for larger taxa. Furthermore, shoaling of the CCD (600 m) at c. 18.5 Ma is likely to have caused a ‘carbonate famine’ (Lyle 2003; Pälike et al. 2012) and may have promoted a transition to surface‐dwelling, smaller taxa.

From eccentricity triggered Antarctic ice sheet expansion in the middle Miocene (Shevenell et al. 2004), to a modest warming interval at the Miocene–Pliocene transition, the Neogene is an important interval in diatom evolution (Lazarus et al. 2014). Climate in the Late Miocene and Early Pliocene is likely to have been driven by a combination of pCO2 fluctuations and tectonically driven partial or complete closure of several oceanic gateways (Lazarus et al. 2014), probably leading to an increased extinction rate within diatoms (Fig. 4D). The gradual closure of the Central American gateway at roughly the Miocene–Pliocene transition redirected tropical Atlantic waters northwards, possibly interrupting surface to deep‐sea circulation and diatom productivity (Berger 2007; Butzin et al. 2011). Rapid evolutionary turnover in the North Pacific Ocean is also noted due to changes in gyral intensity causing extinction rates to increase to 3 to 4 events per 0.5 myr (Barron & Baldauf 1990; Barron 2003). Following this extinction event, almost 50% of contemporary extant species originated (Lazarus et al. 2014), suggesting that diversity–environment relationships during the Miocene and Pliocene were integral to the siliceous diversity of the modern ocean.

Implications of the microfossil occurrence record

Our results selected the NHPP model for estimating preservation rates in each of the groups analysed (Table 2). This reflects a bell‐shaped distribution of preservation rates through the range of a species, corroborating previous research using a different approach that estimated species’ geographical extent (Liow & Stenseth 2007). While we do not separate geographical effects in our results, it is likely that the drivers of the two models are the same: establishment of a species leads to increased geographical range that supports a higher population. This increased range and population then facilitates preservation and sampling at more sites until species decline reverses this trend (Enquist et al. 1995).

The Neptune Database is extensive, comprehensive, and heavily curated, representing one of the highest‐quality palaeontological occurrence datasets. Included fossils are derived from ocean drilling cores, with consistent sampling efforts between samples and groups. The variability in the number of occurrences throughout the interval studied thus represents an accurate picture in fossil preservation bias, both within and between microfossil groups (Figs 1, 2). Despite the extent of this dataset, the results of this study differ from previous studies on the diversity of microfossil groups, missing the fluctuation of nannofossil diversity at OAE‐1a (Erba & Tremolada 2004). The largest change of foraminiferal diversity in the Cretaceous at OAE‐1b (Fraass et al. 2015; Lowery et al. 2020), and the turnover of siliceous taxa at the PETM (Mitlehner 1996; Liu et al. 2011; Oreshkina & Radionova 2014) are three prime examples. The diversification rates for many of the groups reconstructed here are largely static, barring the most major events (Fig. 4). We suggest two reasons behind this inconsistency.

First, preservation biases will change the number of occurrences and the strength of the diversification signal through time, and microfossil preservation will be particularly affected by changes in the marine chemistry (e.g. CCD). While PyRate accounts for variability in preservation (Silvestro et al. 2014), periods of reduced occurrences are accompanied by increased uncertainty in the reconstruction of diversification rates (Figs 3, 4, 5). Previous applications of PyRate to the fossil record of vertebrate groups Rhinoceratoidea and marine mammals show increased uncertainty in their reconstruction of diversification rate where occurrences are less frequent (Silvestro et al. 2014, 2019). We argue that results produced by PyRate, and its underlying Bayesian methodology, benefit greatly from incorporating and presenting this uncertainty. The Neptune Database has highly precisely dated occurrence data (±0.25 myr, as applied in this study), yet the potential smearing effects described above apparently reduced the otherwise lively diversification record of the microfossils included (particularly calcareous forms; Fig. 4). It is expected that similar effects will apply to other groups, and so these effects should be considered when interpreting the results.

Second, many of the previously identified changes that are not picked up in our analyses are regionally limited or temporally sweep through different areas, for example OAE‐1b and the EECO (Premoli Silva & Sliter 1999; Zachos et al. 2008). As we use a global microfossil record, our analyses may thus smooth over all but the most abrupt or widespread events. This highlights the importance of considering the spatiotemporal resolution of occurrence data in analyses (Close et al. 2020). We recommend that future research should consider regional variation in diversity signal as well as global trends. In the study of groups that have a relatively poorer fossil record, less frequent sampling, and/or occurrences, applying the relevant resolution is also important.

Conclusion and future study

Diversity within all microfossil groups fluctuated through time. PyRate tests existing hypotheses associated with major transition periods by ascertaining speciation and extinction rate variations through time. All microfossil groups were uniquely impacted by the different environmental perturbations, but their responses have not been consistent from one extinction to the next. Climatic and tectonic influences appear to be the primary drivers of changes in speciation and extinction (although diversity dependence (Ezard et al. 2011) was not considered); temperature, silicate weathering and nutrient input, ocean circulation patterns altered by oceanic gateways, and upwelling zones strongly influence diversification within both calcareous and siliceous groups. Perturbations associated with warm waters, with a coincident deepening of the CCD, often promoted diversification in planktic foraminiferans and calcareous nannofossils, while cooler SSTs, high latitudes and elevated nutrient inputs are associated with radiolarian and diatom diversification. However, opposing behaviours in calcareous and siliceous taxa were occasionally inconsistent; for example, during Cretaceous anoxia events. The EOT results contrast against previous studies of the glaciation event, with little evidence of a biotic turnover; instead net diversification increased unexpectedly in calcareous forms.

PyRate provides a leap forward in quantitative testing of diversification rate changes throughout geological time. In combination with raw fossil data, phylogenetics and palaeoclimate data, PyRate could provide an increasingly vigorous method to quantitatively determine diversification dynamics. However, it is necessary to acknowledge that PyRate does improve its reliability when using a robust dataset. This study therefore highlights the need for higher resolution records of siliceous biota since they are underrepresented in microfossil databases, particularly at high latitudes, whenever possible, as they are integral to analysing glacial intervals. Additionally, extending the record further back in time could resolve discrepancies associated with OAE events. Nonetheless, PyRate provides a novel method within microfossil diversity to statistically tease out speciation and extinction rates, with the potential for further exciting developments. Since major environmental issues affecting the ocean today include rapid onset warming and precipitous increases in carbon dioxide and the subsequent acidification, it is fundamental to study events such as OAEs and the PETM on both regional and global scales. Investigating mass extinctions such as the K/Pg and major transitions in climate such as the EOT provide insights into ecosystem recovery and evolutionary response to major global change. Determining what drives these changes in diversity, assemblage composition, productivity, and selectivity across events, is imperative to improve our understanding of ancient oceans and how our ocean ecosystems might respond today.

Author contributions

BCM and AJF conceptualized and supervized the project; KMJ carried out the formal analyses, investigation, visualization, and wrote the first draft; all authors contributed to the data curation and methodology. All authors validated and agreed on the final version of the manuscript.

Fig. S1. Sensitivity of PyRate analyses to occurrence age.

Click here for additional data file.

Fig. S2. Number of species occurrences at different sampling depths.

Click here for additional data file.

Fig. S3. Latitudinal differences in sampling frequency of microfossil occurrences.

Click here for additional data file.

Data S1. Example data extraction from Neptune Sandbox Berlin.

Click here for additional data file.

Data S2. Example PyRate Script used for diversification analysis.

Click here for additional data file.