Late Cretaceous ammonoids show that drivers of diversification are regionally heterogeneous

Domain: earth_science
Published: 2024-06-27
Editorial outcome: major_revision

Abstract
Palaeontologists have long sought to explain the diversification of individual clades to whole biotas at global scales. Advances in our understanding of the spatial distribution of the fossil record through geological time, however, has demonstrated that global trends in biodiversity were a mosaic of regionally heterogeneous diversification processes. Drivers of diversification must presumably have also displayed regional variation to produce the spatial disparities observed in past taxonomic richness. Here, we analyse the fossil record of ammonoids, pelagic shelled cephalopods, through the Late Cretaceous, characterised by some palaeontologists as an interval of biotic decline prior to their total extinction at the Cretaceous-Paleogene boundary. We regionally subdivide this record to eliminate the impacts of spatial sampling biases and infer regional origination and extinction rates corrected for temporal sampling biases using Bayesian methods. We then model these rates using biotic and abiotic drivers commonly inferred to influence diversification. Ammonoid diversification dynamics and responses to this common set of diversity drivers were regionally heterogeneous, do not support ecological decline, and demonstrate that their global diversification signal is influenced by spatial disparities in sampling effort. These results call into question the feasibility of seeking drivers of diversity at global scales in the fossil record.

Introduction
Inference of patterns and processes of diversification in the fossil record has a long history of palaeobiological interest1–4. Ever larger compilations of fossil occurrence data and improved chronostratigraphic constraints have enabled precise quantification of the magnitude and timing of interspersed mass extinctions and evolutionary radiations5–7. In turn, this has permitted investigation of the drivers of diversification, broadly divided into intrinsic regulation of taxonomic richness by diversity-dependent mechanisms (the Red Queen Hypothesis)8–12 versus extrinsic, abiotic forcings imposed by the Earth-Solar system (the Court Jester Hypothesis)4,13,14. The influence of the Red Queen remains contentious due to the difficulty of establishing whether biodiversity has ever truly entered a curtailed, diversity-dependent regime15–20, with the additional observation that interactions between life and the Earth-Solar system through geological time have dynamically altered Earth’s carrying capacity21,22. By contrast, a series of Court Jester mechanisms are repeatedly hypothesised to have driven patterns of diversification in individual clades to entire biotas, for example, the linked effects of atmospheric CO2 concentration and temperature, the configuration of the continents, or eustatic sea-level variations19,23–29.

It is also widely recognised that diversity patterns in the fossil record are skewed by geological and anthropogenic biases1,6,30–33, fuelling development of increasingly sophisticated methods for quantifying diversification dynamics from incomplete, biased fossil occurrence data. In the last decade, Bayesian approaches, which couple birth-death and preservation processes have enabled estimation of sampling-corrected origination and extinction rates from fossil occurrences34–37, avoiding the problems of inferring these fundamental rates from extant phylogenies38–42. In turn, lineage birth and death rates can be modelled as functions of their potential drivers43–45, permitting separate consideration of the factors that promoted origination or drove extinction. These models have enjoyed widespread uptake by the palaeontological community46–61 but the paradigm of global biodiversity analysis, whether of individual clades or entire biotas, is challenged by the observation that the spatial distribution and extent of the fossil record varies through geological time31,32,62. This spatial variation induces geographic sampling biases that distort our view of taxonomic richness and diversification rates even after correction for geological sampling biases31. Not only is it inadvisable to treat the fossil record as a representative sample of varying global diversity, with some workers questioning whether global patterns are biologically informative in the first place63, but this also overlooks its well-established biogeographic nuances, necessitating spatially sensitive approaches31,32,64–67.

Given that diversity shows spatial heterogeneity, it is reasonable to also expect this of its drivers and such variation is evident in the present-day biosphere, for example latitudinally structured covariation in irradiance, climate and species richness68–70. Some palaeontologists, however, have continued to compare ‘global’ fossil records to potential drivers without considering regional heterogeneity in diversification processes64,71,72. Here we investigate diversification dynamics of Late Cretaceous ammonoids in a regionalised framework (Fig. 1) to determine whether the drivers of diversity show this expected variation. Ammonoids were an iconic clade of pelagic shelled cephalopods, which originated in the Devonian and survived multiple mass extinction events up until their demise at the Cretaceous-Paleogene (K-Pg) boundary73,74. Late Cretaceous ammonoids provide an ideal study system for our investigation as their fossil record is taxonomically mature and extremely well-sampled compared to many other fossil clades75. Ammonoids additionally present an intriguing controversy regarding whether they were in ecological decline prior to their definitive extinction at the end of the Cretaceous73,76,77, with previous authors highlighting the importance of combining global and regional analyses to tackle this problem78–80. Palaeontologists are also rapidly recognising the wealth of ‘dark data’, present in museum collections but unrecorded in published literature, that has the potential to transform the scale and completeness of palaeobiological analyses66,81–86. We, therefore, take the opportunity to unite publicly available fossil occurrence data with ‘dark’ datasets to address this emerging gap in both knowledge and practice.

Here we show that a regionally heterogeneous interplay of Red Queen versus Court Jester processes underpinned spatial variation in ammonoid diversification in the Late Cretaceous. These results highlight the challenges that face confident inference of the drivers of diversity in the geographically biased fossil record and call for more nuanced consideration of spatiotemporal complexity of the processes that jointly underlie past biodiversity and its present-day geological remnants.

Results
Global ammonoid diversification dynamics
We captured and spatially standardised regional samples of fossil occurrence data for the Antarctic, Tethys, Western Interior Seaway (WIS), Atlantic and Gulf, East Pacific, West Pacific, West Africa, and South Africa (Fig. 1), then estimated genus and species-level origination and extinction dynamics from these data in addition to global dynamics (Fig. 2), using a birth-death-sampling model implemented in a Bayesian framework (PyRate). In all analyses, species-level origination (λs) and extinction (μs) rates were more volatile compared to the genus-level rates (λg, μg), leading to sharper fluctuations in taxonomic diversity (Supplementary Figs. 1–9). The differential quality of the genus and species identifications for the same set of fossil occurrence data led to drops in species richness below genus richness during some intervals in the regional analyses (Fig. 3), even though the former should minimally equal the latter (i.e., at least one species present per genus), highlighting a limitation of our data, but ultimately patterns of diversity at both taxonomic levels were qualitatively similar in all analyses (Fig. 3).

At the global level, μg was broadly stable through the Cenomanian to Santonian stages, punctuated by a spike at end of the Cenomanian (Fig. 2A). A further spike occurred in the early Campanian, followed by a reduced rate for the remainder of the stage. μg then intensified gradually through the Maastrichtian before spiking at the end of the Cretaceous (Fig. 2A). The end-Cenomanian spike, Maastrichtian intensification and end-Cretaceous spike were mirrored by μs but were accompanied by additional pulses at the ends of the Turonian, Coniacian, and Santonian, with stepwise decreases in the intervening background rate (Fig. 2A). μs then increased during the early Campanian to an elevated background rate, although the uncertainties surrounding this rate increase indicate the potential for a much more prominent spike during this interval. λg remained stable through the Cenomanian to the middle Coniacian until a sudden drop, then showed continued stability for the remainder of the Cretaceous (Fig. 2A). λs displayed pulses during the late and terminal Coniacian, preceding reduction to a slightly lower, but stable rate from the early Turonian to the end of the Coniacian (Fig. 2B). λs then followed a pattern of volatile spikes and crashes through the early Cenomanian, a final rapid intensification in the latest Cenomanian, and gradual decline through the Maastrichtian (Fig. 2B).

The outcome of these rates was a broadly stable pattern of genus richness through the Cenomanian to Santonian punctuated by losses at the beginning of the Turonian (coincident with Ocean Anoxic Event 2–OAE2) and Coniacian, a phase of decline through the early Campanian, then a second interval of stability until the end-Cretaceous mass extinction (Fig. 2D). At the species level from the Cenomanian to Santonian, diversity dropped suddenly at the beginning of the stage followed by a more gradual recovery, leading to a pattern of stepwise decline through this interval (Fig. 2D). Species richness crashed successively at the start and middle of the Campanian, with the intervening time displaying high levels of uncertainty in global diversity estimates, recovered suddenly in the early Maastrichtian, then declined in the late Maastrichtian prior to extinction at the K-Pg boundary (Fig. 2D).

Regional ammonoid diversity trends
Spatially varying diversification rates and diversity trajectories demonstrate that ‘global’ ammonoid dynamics in the Late Cretaceous were heterogeneous between regions and taxonomic levels (Figs. 2 and 3). There is a general and expected trend of increasing diversity with region size under the species-area effect, although the WIS is the most diverse region despite being smaller than Tethys (Fig. 2E, F). Qualitatively, the global trend is primarily a synthesis of patterns from the larger, historically well-sampled regions within the Global North (Tethys, WIS, Atlantic and Gulf, East Pacific), although these regional patterns still display their own nuances (Fig. 2E, F). Notably this correspondence between the global curve and specific regional trends is still apparent despite the former including data that was not present in any of the other regions, for example, more sporadically sampled occurrences from South America, Australia and Madagascar (Fig. 1).

Declines in genus and species richness occurred during the interval containing OAE2 in Tethys, Atlantic and Gulf, WIS, and East Pacific, but the onset of declining species richness was earlier in the latter (Fig. 3F). Trends outside of the Global North were more varied, with genus richness increasing substantially in West Africa through the Cenomanian (Fig. 3E), but only marginally in the West Pacific and Antarctic (Fig. 3A, G), while South Africa experienced no change (Fig. 3C). Genus richness then declined from the late Cenomanian into the early Turonian in the West Pacific, while West Africa and the Antarctic experienced minor diversity losses during OAE2, and South Africa instead experienced a rapid increase across the same time interval. Sharp species loss took place in all four regions in the Global North at the end of the Turonian, but corresponding loss of genera did not take place in the WIS (Fig. 3D). Declines in species and genus richness also took place in South Africa, West Africa and East Pacific, while the Antarctic displayed little to no change across this boundary. In South Africa, this led to continued decline into the Campanian while the other regions showed recovery during the intervening Coniacian and Santonian stages. Generic richness declined gradually in the WIS, Tethys and the Atlantic and Gulf at the beginning of the Campanian, while the East Pacific experienced a much sharper crash several million years later (Fig. 3F). Species-level dynamics at this time were even more varied, with a gradual decline in Tethys, a sharp spike then equally sudden crash in the East Pacific and the Atlantic and Gulf, and turbulent, highly uncertain fluctuation in the WIS (Fig. 3B, D, F, H). Recovery of diversity varied in its intensity across the Global North in the late Campanian and Maastrichtian, with Tethys displaying the greatest increase in species richness and most stable genus richness through the latter interval (Fig. 3B). Away from the Global North, genus richness also recovered towards the end of the Cretaceous (Fig. 3A, C, E, G), but this recovery began earlier in the Antarctic and West Pacific in the late Campanian, compared to the Maastrichtian in South and West Africa.

Despite their greater volatility, trends in species richness in these regions generally support the trends recovered at the genus-level (Fig. 3).

Biotic and abiotic drivers of diversification
We used a multivariate birth-death model (PyRateMBD) to identify time-continuous correlates of our regional origination and extinction rates (Fig. 4). Biotic factors considered were the effects of ammonoid diversity on their own origination and extinction rates (diversity dependence), the influence of food availability (based on nannoplankton and planktic foram diversity) and the pressure of pelagic predators/competitors (based on ichthyolith abundance), while abiotic factors were long-term changes in sea level, atmospheric CO2 concentration, and sea surface temperature. At the global level, μg and μs showed significant positive relationships with long-term sea level (+ve), significant negative relationships with ichthyolith abundance and sea surface temperature, μs displayed significant positive relationships with ammonoid and planktic foram diversity, and a significant negative relationship with atmospheric CO2 concentration.

At a regional scale, species-level rates generally displayed a greater number of significant relationships with our chosen drivers compared to genus-level rates (Fig. 4), with significant relationships also becoming more common in the smaller sampling regions. Extinction rates at both taxonomic levels most consistently showed significant negative and positive relationships with sea surface temperature and long-term sea level, respectively, aside from at the genus level in Antarctica, where the relationship was instead negative. Ammonite diversity also displayed recurring significant negative relationships with origination rates at the species level and the genus level to a lesser extent. Beyond these broad trends, however, drivers of diversification showed complex variation between regions and taxonomic levels, encompassing cases where relationships remained significant but changed their sign or instances where significance changed entirely (Fig. 4).

Taxonomic predictors of extinction risk
Alongside time-continuous predictors, we also used a multi-trait extinction model (PyRateMTE) to investigate the effects of discrete traits on regional ammonoid extinction risk (Fig. 5). Rather than define and measure a series of ecologically meaningful traits across all Late Cretaceous ammonoid taxa, a task that may not even be possible for some taxa depending on the quality of their fossils, we instead used their suborder and superfamily classifications as these categories are considered to represent distinct ecomorphological guilds which summarise a range of taxonomically diagnostic, functional ecological characteristics83–85. At the global scale, suborders displayed predictive power for both μg and μs, while superfamilies were only predictive for the latter (Fig. 5). Similarly, suborders displayed greater predictive power for regional μg and μs with significant effects present in four regions, compared to two and three regions for superfamilies at the genus and species levels respectively. Predictive power varied regionally depending on the taxonomic level of the extinction rates under consideration and on the predictor clades themselves, although the Perisphinctina showed the greatest vulnerability across regions and taxonomic levels (Fig. 5).

Methods
Ammonite occurrence data compilation
We compiled a database of Late Cretaceous ammonite occurrences from several sources: the Palaeobiology Database (PBDB), a published occurrence dataset from Yaccobucci80, and three unpublished datasets compiled by the authors. This includes vast new occurrence data made accessible through large-scale museum digitisation initiatives123, alongside occurrences compiled from literature sources not present in the PBDB. To avoid edge effects in downstream analyses, we additionally included all ammonoid occurrences in the PBDB for the preceding Albian stage of the Early Cretaceous. All dataset compilation and revisions were conducted in R (v4.2.2)124. All Cretaceous cephalopod occurrences were downloaded via the PBDB application programming interface on 31/05/22 to avoid potential exclusion of ammonite occurrences misclassified within other cephalopod orders or with poor stratigraphic constraints. Occurrence chronostratigraphy in each dataset was updated according to the Geologic Timescale 2020 (GTS2020)125 using the chrono_scale() function from the fossilbrush R package126, then palaeocoordinates calculated from their midpoint ages under the PALEOMAP plate rotation model126–128 using the palaeorotate() function with default settings from the palaeoverse R package129. A small number of occurrences could not be rotated and were discarded from the dataset. PBDB occurrences all had locality information present, while missing localities in the other datasets were assigned according to unique combinations of their present-day longitude, latitude, and maximum and minimum stratigraphic ages. All four datasets were then combined, incorporating their occurrence and locality IDs, classifications at the species, genus, family, and order levels, maximum and minimum chronostratigraphic ages, and present-day and palaeogeographic longitude-latitude coordinates.

Next, we revised the taxonomy of the Late Cretaceous occurrences in our composite database. Subgenus classifications recorded within the genus column (commonplace within the PBDB) were deleted, leaving only the genus-level assignment. Suspect and outdated genera were manually reassigned, and family- and order-level classifications of all genera exhaustively updated according to the current literature and the taxonomic expertise of the authors. Further inconsistencies arising from spelling variations between genus names and any remaining discrepancies in higher-level taxonomy were flagged using the check_taxonomy() function from fossilbrush and manually resolved. Uncertain species entries (e.g., those with cf. or “?” modifiers) were stripped back to genus-level. Occurrences without genus-level assignments, those named from ammonite jaw elements, and all remaining non-ammonoid occurrences were discarded. Superfamily and suborder classifications were assigned to each occurrence, using the scheme of Yacobucci80 for the latter. All taxonomic revisions are documented in the electronic supplement.

Finally, we inspected the stratigraphic ranges of the remaining genera to identify and eliminate incorrectly classified occurrences using several methods. Firstly, stratigraphically suspect occurrences were detected by comparing occurrence stratigraphic ages in the composite database to two reference databases of ammonoid stratigraphic ranges: the Sepkoski Compendium130 and an unpublished compilation constructed with taxonomic expertise by J.D.W. The Sepkoski Compendium with GTS2020 dating was obtained from fossilbrush, while the compilation from J.D.W was updated to the same chronostratigraphic standard as above. Database comparisons were performed using the flag_ranges() function from fossilbrush, which identifies occurrences falling outside their reference ranges. Secondly, misclassified or poorly dated occurrences contributing to anomalously long tails (the portions of a genus stratigraphic duration falling outside a given confidence interval in the distribution of its occurrences) were flagged using the pacmacro_ranges() function from fossilbrush. Flagging by this latter method was performed using the kernel density and histogram implementations with tail proportions of 30, 35 and 40% (i.e., the % of the stratigraphic duration occupied by the tails) with the default step size of 0.1 Ma for densifying occurrence durations and default tail confidence intervals of 5% (i.e., the lowermost and uppermost portions of a stratigraphic duration assigned to the tails). Flagged occurrences from both methods were then manually inspected to determine their validity, leading us to discard them from the dataset, revise their genus-level classifications, or update their chronostratigraphic ages. Lastly, all remaining occurrences falling outside the Albian to Maastrichtian stages were discarded from the composite database along with all occurrences with stratigraphic age uncertainties exceeding 10 million years. All rejections are documented in the electronic supplement.

To our knowledge, our final database is the largest and most spatiotemporally comprehensive compilation of Late Cretaceous ammonoids to-date, with fully revised and internally consistent taxonomy, containing 19,536 mostly substage-level occurrences (mean and median age uncertainties of 3.5 and 2.5 Ma, respectively) classified at the order, suborder, superfamily, and genus levels, of which 15,149 (77.5%) are additionally classified to species level, and 9962 (50.1%) are occurrences from literature not yet incorporated into publicly accessible biodiversity databases, or true ‘dark’ data from museum collections. Consequently, it is well-suited for analysis of their diversity dynamics in the interval leading up to their demise at the end-Cretaceous mass extinction.

Spatial subsampling
To determine how ammonoid diversity dynamics differed between different parts of the globe and to circumvent potential biases introduced by spatially heterogenous sampling in the fossil record, we subsampled our composite dataset into eight biogeographically relevant regions with ammonoid occurrence records covering the full extent of the Albian to Maastrichtian: Tethys (corresponding to Europe and North Africa), West Africa, South Africa, Antarctic (corresponding to the James Ross Basin), West Pacific (corresponding largely to Japan), East Pacific (corresponding to the North American Pacific seaboard), the WIS, and the Atlantic Gulf (corresponding to the North American Atlantic and Gulf Coast). These regions were largely identified by plotting the palaeogeographic distribution of our composite dataset at five million-year intervals through the Late Cretaceous onto the PALEOMAP highstand reconstructions of Scotese and Wright (2018)128, while the decision was made to separate the otherwise geographically continuous Atlantic and Gulf from the WIS as they constitute different water masses and show marked endemic differences in their ammonoid faunas131,132, supporting their designation as distinct bioregions.

We followed the methodology of Flannery-Sutherland et al.65 to eliminate spatial sampling bias in each region, using the corresponding R functions provided in their ESM. We constructed fixed-size spatial windows to subsample the occurrence data for each region using the spacetimewind() function, permitting each window to slide in a constant direction between geological stages to accommodate continental drift of their target biogeographic regions. To remove any remaining effects of spatial sampling heterogeneity, the data in each region was further subsampled within each geological stage to target spatial extent based on minimum spanning tree (MST) length. MST length is correlated with a variety of other spatial properties12 and so is a robust metric for the spatial extent of a geographically resolved dataset. Subsampling was implemented using the spacetimestand() function which spatially bins the occurrence data within a given region and geological interval using a hexagonal grid generated by the hexagrid() function from the icosaR package133, thus avoiding any latitudinal distortions in grid cell area, calculates the MST connecting all occupied grid cell centres, and sums the great circle lengths of each of its segments. The function then progressively discards grid cells forming the tips of the MST, prioritising those with the fewest occurrences, until the closest possible match to the target MST length is achieved. After spatial standardisation, apparent changes in taxonomic diversity within a region may still be driven by temporally heterogeneous sampling. Consequently, a temporally standardised diversity curve should not be driven by any residual fluctuations in spatial extent if spatial standardisation was effective. We therefore calculated sampling-corrected genus diversity at stage level through time in each region by shareholder quorum subsampling (SQS) with a sampling quorum of 0.5, implemented within spacetimestand() using the estimateD() function from the iNEXT R package134 (Supplementary Figs. 1–9). SQS is the most robust subsampling method for estimating diversity from fossil occurrence data and so is well-suited to this purpose135. We then conducted one-tailed Pearson and Spearman correlations between genus-level SQS diversity, and longitude range, latitude range and MST length for each spatially standardised region to determine if increased spatial extent was associated with increased diversity. To determine how global ammonoid diversity dynamics related to its regional components through the Late Cretaceous, we also spatially standardised the entire composite dataset and tested for correlations between genus-level SQS diversity and spatial extent, following the above protocol (Supplementary Data 1). Collectively, our regions encompass 69.3% of the occurrences in the total dataset, while the standardised global dataset encompasses 94.5%.

Diversity dynamics
Diversity dynamics were quantified for each spatially standardised regional dataset in a hierarchical Bayesian framework using PyRate (v 3.0)34–36. PyRate jointly estimates the preservation rate underlying the occurrence records of an incompletely stratigraphically sampled set of fossil lineages, sampling-corrected times of origination and extinction of those lineages, and the corresponding origination and extinction rates arising from their true stratigraphic durations35. PyRate models origination and extinction rates using a birth-death process describing the expected numbers of lineage origination and extinction events per unit of time, and preservation rates using a Poisson process describing the expected number of fossil occurrences per lineage per unit of time, with posterior distributions of the parameters of each process estimated using Markov Chain Monte Carlo (MCMC). PyRate can implement several different Poisson preservation processes, including the potential for their variation between lineages and through time, and uses reverse jump MCMC sampling to estimate the number and timing of rate shifts in the birth-death process, along with their statistical support36. By jointly modelling the processes responsible for a set of taxonomically resolved fossil occurrence data, along with their uncertainties, PyRate provides substantially more accurate estimates of origination and extinction rates through geological time compared to traditional methods36,136 and the method has consequently received widespread use by palaeontologists working on a wide variety of different clades (e.g., ferns45, dinosaurs47, marine megafauna137). Input files for PyRate were generated in R from the global and regional datasets at genus and species levels. Each input file contains 10 age-randomised replicates of the original dataset where the age of each occurrence is randomly fixed within the bounds of its stratigraphic uncertainty. This procedure was performed such that occurrences from the same locality were consistently assigned the same randomised age in each replicate, following previous recommendations138. The best-fitting preservation model for each dataset was identified by maximum likelihood, using the -PPmodeltest command of PyRate) to compare Akaike’s information criterion (AIC)139 scores for a homogenous Poisson process (HPP), an inhomogenous Poisson process (NHPP), and a time-variable Poisson process (TPP) consisting of user-defined rate shifts bounding constant-rate intervals (Supplementary Data 2). For the TPP model, rate shifts (PyRate option -qShift) were set as the boundaries between the geological stages of the Late Cretaceous. The TPP model was identified as best-fitting for each regional dataset, except for Antarctica at the genus-level where the HPP model was best-fitting. Initial analysis under the HPP model, however led to unrealistic extension of Antarctic ammonoid genus ranges into the Cenozoic by several million years, despite the empirical fossil record clearly documenting their rapid demise at the end of the Cretaceous. Consequently, we elected to use the TPP model in all cases, as the HPP model was unable to adequately capture this particular case of total loss of preservation rate in the geological past resulting from sudden and complete clade extinction well before the present day. PyRate uses a gamma distribution with shape=1.5 as the default prior on the vector of bin-wise preservation rates when using the TPP model. We assigned a vague exponential hyper-prior on the rate parameter of the gamma distribution (PyRate option -pP 1.5 0) and permitted PyRate to estimate the empirical rates during the analysis, reducing subjectivity in prior selection and allowing the analysis to better adapt to the input data. Each analysis was additionally constrained to only search for rate shifts within the Albian to Maastrichtian to further reduce edge effects induced by the temporal bounds of our data (PyRate option -edgeshift 113.266). We ran PyRate on each age-randomised dataset for 100 million generations, sampling every 50,000, using the reverse jump MCMC algorithm to efficiently search for shifts in origination and extinction rates rather than set any such shifts a priori (PyRate options -A 4). Log files were analysed in Tracer (v1.7.2)140 to determine that stationarity had been achieved. The first 10% of each log file was discarded as burn-in, then each log file combined into a final log using the -combLog command of PyRate, enabling the age uncertainties in the original fossil occurrence data to be incorporated into our results. Origination and extinction rates, along with Bayes Factor (BF) support for any rate shifts, and mean times of origination and extinction for each lineage were calculated for each region using the -plotRJ and -ginput commands of PyRate (Supplementary Figs. 10–18). To assess the statistical significance of rate shifts, we took log BF > 2 and log BF > 6 as positive and strong support respectively141. Range-through diversity was calculated from the preservation-corrected lineage stratigraphic durations at 0.5 Ma intervals in R.

Drivers of origination and extinction
We used a multivariate birth-death (MBD) model to identify potential drivers of ammonoid diversification dynamics, implemented in PyRateMBD45. Origination and extinction rates within the birth-death process are modelled as linear or exponential correlations of time-continuous variables. Rates are derived from the preservation-corrected lineage durations estimated by a previous PyRate analysis, circumventing the need to account for sampling heterogeneity within the MBD model itself. The strength and sign of the correlation functions are jointly estimated for each variable using MCMC, with a horseshoe prior on each variable controlling for over-parameterisation and the potential effects of multiple testing by shrinking each correlation parameter around zero. Significant correlations consequently show shrinkage weights which differ substantially from zero (W > 0.5)45.

For biotic influences, we considered the effects of ammonoid diversity on their own origination and extinction rates (diversity dependence), the influence of nutrient availability and the pressure of pelagic predators/competitors. Ammonoid diversity was calculated directly within PyRateMBD from the preservation-corrected lineage durations used by the analysis. The diversities of calcareous nannoplankton and planktic foraminifera142 were used as proxies for particulate nutrient availability in the water column. Ichthyolith accumulation rates143,144 was used as a proxy for the abundance of large pelagic fish which were the most likely predators or competitors of ammonites145,146, as opposed to larger apex predators like mosasaurs. For abiotic influences, we considered long-term sea level106, atmospheric CO2147, and sea surface temperature148. All variables were taken from the electronic supplementary materials of each source, linearly interpolated at 0.1 Ma intervals, then rescaled to the range [0,1] to avoid biases induced by differences in their absolute magnitudes.

We ran PyRateMBD on the 10 sets of lineage origination and extinction times estimated in each regional genus- and species-level PyRate analysis for 50 million generations, sampling every 50,000, using exponential then linear correlation functions. To eliminate the potential impacts of total and near instantaneous lineage extinction at the end of the Cretaceous from affecting our results, all MBD analyses were temporally restricted to 100.5–68 Ma. Log files were analysed in Tracer to check for stationarity, then combined log files produced as above. We then compared the fits of the exponential versus the linear MBD models for each dataset in R using AIC for MCMC samples149 and Bayes factors141 calculated from the harmonic means of each model likelihood. In almost all cases we found that the fits of these equally complex models could not be readily distinguished from one another (ΔAIC < 2, log BF < 2; Supplementary Data 3), aside from the genus-level analyses for the Antarctic where the exponential model was favoured (log BF = 5.31). Regardless, the statistical significance of virtually all the correlation parameters were unchanged by the choice of correlation function, aside for the genus-level analyses for the West Pacific where the exponential model was marginally better fitting (ΔAIC = 0.05). Consequently, we interpreted results from the exponential models in all cases (Supplementary Figs. 19–27, Supplementary Data 4–12), although the results from linear models are also available in the electronic supplement (Supplementary Figs. 28–36, Supplementary Data 4–12).

Ecological correlates of extinction
In addition to examining the influence of temporally varying abiotic and biotic drivers on ammonoid diversification, we additionally tested for discrete ecological correlates of their extinction risk using a multi-trait-dependent extinction (MTE) model, implemented in PyRateMTE150. Likelihood of lineage durations within the birth-death process is modelled as a function of discrete, trait-based multipliers and a mean extinction rate. As with the MBD model, lineage durations are derived from the preservation-corrected speciation and extinction times estimated in a previous PyRate analysis. The trait multipliers enable lineage-wise deviations from the mean extinction rate, with mean extinction rate and multiplier strengths estimated using MCMC. Model priors and hyperpriors control for over-parameterisation by shrinking the multipliers towards equal values (i.e., no differential effects on extinction rate), along with Bayesian variable selection, which places 95% of the prior probability on the absence of an effect for each trait then removes traits deemed to have no substantial effect on rate variation. Traits selected as having an effect, which additionally show posterior frequencies greater than 51.4% (equivalent to log BF > 6) are taken as statistically significant.

Suitable ecological traits may be derived directly from categorical aspects of an organism’s life mode (e.g., filter feeder, deposit feeder, active hunter etc.) or from discretisation of a continuous trait (e.g., occupied depth range). Many of traits are difficult to directly characterise for ammonoids, but their higher taxonomic classifications are considered to map onto ecologically distinct life modes. Consequently, we used the suborder and superfamily classifications (five and twelve categories respectively) as ecological proxies of extinction risk in the MTE model, under the assumption that these classifications subsume sets of clade diagnostic characters that are also responsible for their distinct ecological modes. We ran PyRateMTE on the 10 sets of lineage origination and extinction times estimated in each regional genus- and species-level PyRate analysis for 10 million generations, sampling every 10,000. To eliminate the potential impacts of total and near instantaneous lineage extinction at the end of the Cretaceous from affecting our results, all MBD analyses were temporally restricted from 100.5 Ma to 68 Ma. Log files were analysed in Tracer to check for stationarity, combined log files produced as above, then each model parameter summarised as their means (Supplementary Data 13–14).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "I was really enjoyed to read this manuscript!\nThe driver of diversification is one of the most attractive topics in both paleontologists and biologists.\nThe authors adopted the state-of-the-art Bayesian method, PyRate, to recovery “global” and regional\ndiversity dynamics of Late Cretaceous ammonoids, and found strong heterogeneity in diversity\ndynamics among different regions and highly variable interactions between diversity and different\nbiotic and abiotic factors in variable tempo-spatial scales. The manuscript is will written, and the\nfigures are excellent.\nI have only a few concerns/comments as follow:\nI still can’t accept that species richness can lower than genus richness in some intervals in some\nregions. This is not true based on our natural classification system. The authors could treat, for\ninstance “Haploceras sp.”, as a species to overcome this problem.\nLines. 144-169: It’s a bit hard to read here. Why don’t you describe the diversity dynamics in the\neight regions in one paragraph. For instance: “Declines in genus and species richness occurred during\nthe interval containing OAE2 in all four,”, At first glance, I was confused, what are the four regions? I\ncan only find the answer (East Pacific, Tethys, WIS and Atlantic and Gulf) when I read the whole\nparagraph.\nLine 353. “Middle” => “Early”.\nLine 577-578: Do you mean five superfamily and twelve suborders? Please check.", "reviewer_2": "Spatiotemporal variations of biodiversity impact our understanding of past evolution. This manuscript\nfocuses on the debate of what drives diversification and whether we can faithfully identify these\ndrivers from the fossil record.\nI do not have major concerns on the technical part of the paper. However, the author(s) soundly\ncriticized past studies and, in some sentences, all palaeontologists. I am not sure if this is the best\nsolution to promote this paper, but this is up to the authors. From an absolutely objective perspective,\nmost researchers need a simple solution for many applications of biodiversity changes e.g., in\npalaeoenvironment studies. The appoarch provided by the author(s), though impressive and better, is\nnot simple, nor a solution to the general problem.\nNevertheless, this is a well-written manuscript that will add new insights into our understanding of\ndrivers of past extinction, recovery, and diversification.\nLine 69-71: True, but the words are rather harsh.\nLine 137: WIS first occurred here. But the abbreviation was only explained in line 421.\nLines 144-155: When reading this paragraph, the question popped up: Would the bias of sampling\nefforts completely rule out?\nLines 242-243: the sentence “regression-driven loss of habitat space is more commonly associated\nwith elevated extinction rates” needs references.\nLines 248-250: difficult to follow, explain why prevalence in the smaller sampling regions and absence\nin the larger occur.\nLines 250-252: regional ammonoid origination and extinction rates showed marked heterogeneity.\nHow much difference can be inferred as marked heterogeneity? Please explain.\nLines 288-293: Not very clear here. Taxonomic assignments are not affected by measurement\nlocation. This seems to contradict the idea of spatial heterogeneity mentioned earlier. And the\nfollowing sentence “fossil occurrences associated with taxonomic assignments display geographically\nbiased sampling” doesn’t seem to make sense. Explain in more detail and give some references.\nLines 309-312: What then was the reason (driver) for ammonoid extinction?\nLine 311-312: the sentence “the hypothesis of long-term global diversity decline prior to their final\nextinction at the K-Pg boundary” needs references.", "reviewer_3": "I laude the authors for their focus on regional diversity trends as a way to understand “global”\ndiversification patterns of ammonoids during the Late Cretaceous. As Geerat Vermeij stated in 2003,\nglobal diversity does not really mean what it implies given the well-established biases in the known\nfossil record at the global scale. This study is an interesting and important contribution to the\nunderstanding of a major group of animals as they approached their ultimate demise. The results\nindicate that despite parsing out global data into distinct regions, ammonoid diversity was not waning\nduring the Late Cretaceous as many workers have suggested. Moreover, the results suggest that\ndrivers of regional diversity patterns are more likely driven by regional drivers, such as temperature\nand sea-level change.\nOverall, this study highlights an exciting approach for assessing the evolutionary history of\nLate Cretaceous ammonoids and potential application to other major invertebrate groups in the fossil\nrecord. The methods also seem sound to me, particularly in their approach for controlling sampling\nbiases. I only have a few minor comments to make regarding the discussion and methodology,\notherwise I believe this manuscript should be accepted once these issues have been addressed.\nMinor comments:\nSea-level change has long been attributed to be a driver of both extinction and origination in\ndynamics in the fossil record. However, evidence suggests that sea-level change itself does not\nmeaningfully effect macroevolutionary dynamics, regardless of spatial scope (Holland 2012, Paleobiology).\nGiven that ammonoids are pelagic animals and even less likely to be affected by sea-level change, I fail to\nsee the reason to include sea-level analysis in this study. Sea-level change will have a stronger\nrelationship with the types of facies preserved, with relatively deeper-water facies being deposited\nduring times of sea level rise, which are more likely to preserve ammonoid fossils than shallower\nwater facies. Which leads into my next comment.\nAgain, while I appreciate the effort of the authors to consider the regional heterogeneity of\ndiversity dynamics, there is no mention regarding the nature of the stratigraphic record at the\nregional scale. While ammonoids are preserved across the majority of marine facies, they are most\nlikely to be live and preserve in relatively deeper-water environments (80–400m; Ritterbush 2014, J. of\nZoology; Lukeneder et al. 2010, Earth Planet. Sci. Lett.). As such, the structure of the stratigraphic\nrecord regionally would have a greater effect on diversity rather than sampling biases. I wonder if\nstratigraphic packages regionally would potentially have a 1st order control on diversity. I do realize,\nhowever, the difficulty in actually testing this, given the poor state of facies/deposition environment\ndata in the PBDB and lack of global stratigraphic packages in Macrostrat. Just a thought.\nLastly, while sea-surface temperatures might have mixed effects on diversity in this study, it is\npossible that temperature in conjunction with oxygen might be a better environmental component to\nassess. While modern cephalopods are sensitive to temperature, they have behavioral adaptations to\ncompensate for increased temperatures, but not for temperature-dependent hypoxia.\nOtherwise, great and interesting study!"}, "author_response": "Reviewer comments are given in regular text, our responses are in red.\nReviewer #1:\nI was really enjoyed to read this manuscript! The driver of diversification is one of the most\nattractive topics in both paleontologists and biologists. The authors adopted the state-of-the-\nart Bayesian method, PyRate, to recovery “global” and regional diversity dynamics of Late\nCretaceous ammonoids and found strong heterogeneity in diversity dynamics among\ndifferent regions and highly variable interactions between diversity and different biotic and\nabiotic factors in variable tempo-spatial scales. The manuscript is well written, and the\nfigures are excellent. I have only a few concerns/comments as follow:\nWe are delighted by the reviewer’s positive response to our work and thank them for their\nendorsement.\nI still can’t accept that species richness can lower than genus richness in some intervals in\nsome regions. This is not true based on our natural classification system. The authors could\ntreat, for instance “Haploceras sp.”, as a species to overcome this problem.\nWe agree that genus richness exceeding species richness is problematic under the nature of\nsystematics, but we do not posit that it is a true pattern. Instead we recognise it as an\nartefact of empirical palaeontological occurrence datasets that results from a deficit between\noccurrences identifiable at species level versus genus level. That this artefact is present in a\nclade as well sampled as ammonoids suggests that it is likely present elsewhere in the fossil\nrecord. The reality of incomplete sampling effects between taxonomic levels is an ugly one,\nbut we assert that this is an important limitation of empirical palaeontological occurrence\ndata which should be presented honestly, given the valid methodological concerns it raises.\nWe do not believe that adding placeholder species to a formal PyRate analysis would be a\njustifiable approach for three reasons:\n1. PyRate models preservation rates to help estimate true lineage durations, so the\ntaxonomic resolution of an occurrence will affect how it contributes to genus- versus species-\nlevel preservation rates. Dragging genus level occurrences to the species level would imprint\ngenus-level rates on the modelled species-level rates, distorting estimated species durations\nand their origination-extinction dynamics.\n2. While the purpose of placeholder species would be to meet the required richness ratio of\nat least one species per genus, it runs the risk of lumping genuinely separate (albeit\nunidentified) species into one composite record without any way of detecting this distortion.\n3. The logic behind assigning all genus-level occurrences of a taxon to one placeholder\nspecies means that that all occurrences which contribute information on the genus-level\nduration of that taxon should be treated in the same way. This would necessitate inclusion of\nspecies-level occurrences within their genus-wise placeholder species, despite already\nrepresenting their own lineages, and resulting in an analysis which is then largely collapsed\nback to the genus level. Conversely, inclusion of solely genus-level occurrences in\nplaceholder taxa would be an arbitrary choice to make species-level dynamics reflect genus\nlevel dynamics, even though they are ultimately different signals.\nA simpler, post hoc solution is for the reader to interpret the genus level curve as a lower-\nbound estimate on true species richness in instances where species richness falls below\ngenus richness. However, the high correlation values between genus and species richness\nin each region suggest that relative changes in the later still faithfully reflect changes in the\nformer, lending support to admittedly imperfect signal of species diversification we recover\nhere. Ultimately, we believe that it is better to present both sets of curves as they stand,\ngiven there is still no robust linkage in the literature between genus and species-level\ndynamics, and acknowledge their imperfections clearly in the discussion. To aid in this, we\nhave made explicit reference to the issue in the caption for fig. 3 (L1056–1058) and direct\nthe reader to the longer discussion of this issue in the text (L340–347).\nLines. 144-169: It’s a bit hard to read here. Why don’t you describe the diversity dynamics in\nthe eight regions in one paragraph. For instance: “Declines in genus and species richness\noccurred during the interval containing OAE2 in all four,”, At first glance, I was confused,\nwhat are the four regions? I can only find the answer (East Pacific, Tethys, WIS and Atlantic\nand Gulf) when I read the whole paragraph.\nWe agree that our original phrasing did not aid with manuscript clarity here. We now\ndescribe the dynamics in one paragraph as the reviewer suggests (L144-169). Our original\nstructure was designed to highlight how regions away from the Global North depart from the\n‘global’ diversity signal, so we have tried to retain this in the revised structure as we view it\nas a key result in our study.\nLine 353. “Middle” => “Early”.\nChange made (L401)\nLine 577-578: Do you mean five superfamily and twelve suborders? Please check.\nQuite right, these were misordered, change made (L625)\nI tried several times and can't access the code via the link provided above.\nWe apologise for this issue. We have revised our ESM in light of this review and have made\nthis new version available via the manuscript submission portal using Nature’s FigShare\nservice. We hope that this has addressed the problem, but if not then please get the editor to\ncontact us and we will endeavour to make the full ESM available in another way\nReviewer #2:\nSpatiotemporal variations of biodiversity impact our understanding of past evolution. This\nmanuscript focuses on the debate of what drives diversification and whether we can faithfully\nidentify these drivers from the fossil record. I do not have major concerns on the technical\npart of the paper. However, the author(s) soundly criticized past studies and, in some\nsentences, all palaeontologists. I am not sure if this is the best solution to promote this\npaper, but this is up to the authors.\nWe are pleased that the reviewer views our paper as technically sound. We agree that our\n that our\noriginal phrasing could be interpreted as a little damning of all palaeontologists, which is\ncertainly not our intention. We have modified our phrasing to reduce this possibility. We\ninstead view the need for spatially sensitive analysis as the ‘promoting’ point for the work.\nFrom an absolutely objective perspective, most researchers need a simple solution for many\napplications of biodiversity changes e.g., in palaeoenvironment studies. The approach\nprovided by the author(s), though impressive and better, is not simple, nor a solution to the\ngeneral problem. Nevertheless, this is a well-written manuscript that will add new insights\ninto our understanding of drivers of past extinction, recovery, and diversification.\nWe thank the reviewer for their positive perspective on the contributions of our manuscript to\nthe field. We agree that our analytical strategy is not the most straightforward to implement,\nbut this is largely because we aimed to analyse rates rather than standing diversity.\nMeasuring the latter in a spatially sensitive framework is much easier to achieve using\nstandard analytical tools and R packages and is a more achievable solution for future\nworkers. We view our take home message as the need for spatially sensitive analysis to\neffectively tackle questions of extinction, recovery and diversification in the fossil record.\nLine 69-71: True, but the words are rather harsh.\nWe have revised our phrasing to ‘some palaeontologists’ (L69). This is a more accurate\nreflection of the state of the field and implicitly acknowledges the excellent efforts of other\npalaeontologists to move away from global analyses towards spatially sensitive frameworks.\nLine 137: WIS first occurred here. But the abbreviation was only explained in line 421.\nWe have ensured that the abbreviation is now provided at the first instance of the complete\nphase (L95).\nLines 144-155: When reading this paragraph, the question popped up: Would the bias of\nsampling efforts completely rule out?\nOur analytical strategy is designed to minimise various sources of sampling bias as much as\npossible, but we agree that it is probably impossible to eliminate its effects entirely. We do\nnot address this point within the lines referenced by the reviewer to keep the results free\nfrom any discussion but have referred to the issue of residual sampling bias as a new\nsection later in the manuscript (L324–340). We believe that this inclusion greatly improves\nthe quality of our work and thank the reviewer for suggesting that we consider the issue of\nsampling in more detail\nLines 242-243: the sentence “regression-driven loss of habitat space is more commonly\nassociated with elevated extinction rates” needs references.\nWe have added references to support the impact of regression on extinction, but we have\nalso caveated that this relationship is not necessarily consistent across space, which helps\nto clarify why rising sea levels might instead coincide with elevated extinction rates (L243–\n250).\nLines 248-250: difficult to follow, explain why prevalence in the smaller sampling regions and\nabsence in the larger occur.\nThis was a remark on the prevalence of significant correlations between sp/ex and\nammonoid diversity in the smaller sampling regions, so we have added a reference to the\nfigure and reworked the phrasing for clarity. We have also added a comment with references\non a possible explanation for this finding based on the literature – that RQ effects are\noverwhelmed by CJ effects at larger spatial scales (L260–264).\nLines 250-252: regional ammonoid origination and extinction rates showed marked\nheterogeneity. How much difference can be inferred as marked heterogeneity? Please\nexplain.\nWe have removed ‘marked’ as this phrasing was vague with no clear basis for determining if\ndifference is marked or otherwise (L265).\nLines 288-293: Not very clear here. Taxonomic assignments are not affected by\nmeasurement location ….\nWe were a little uncertain of whether the reviewer is giving their interpretation of what we\nsaid in the manuscript in relation to the next point raised, or simply stating the above point\nplainly. We aimed to make the same point in the manuscript (L304–305), so we hope that we\nhave understood the reviewer here and that our response below addresses their intended\nconcern.\nLines 288-293 (cont): …. This seems to contradict the idea of spatial heterogeneity\nmentioned earlier. And the following sentence “fossil occurrences associated with taxonomic\nassignments display geographically biased sampling” doesn’t seem to make sense. Explain\nin more detail and give some references.\nWe disagree that independence of taxonomic assignment from spatial location contradicts\nthe idea of spatial heterogeneity. Our point is that the sampling of fossils as data points is\nspatially biased, but that this spatial bias will not have any impact on the taxonomic identities\nassigned to those data points (L304–309). Baculites sp. identified in one location will belong\nto the same taxonomic unit as Baculites sp. from another location, even if one of those\nlocations is disproportionately better sampled than the other. We have revised our phrasing\nto help ensure that this point is clear.\nLines 309-312: What then was the reason (driver) for ammonoid extinction?\nSentence and references added for the current stance in the literature on ammonoid\nextinction at the K-Pg (L354–358). We did this at a different point in the manuscript to the\nlines where the reviewer raises their comment to keep our conclusions as free as possible\nfrom inferences made elsewhere in the literature.\nLine 311-312: the sentence “the hypothesis of long-term global diversity decline prior to their\nfinal extinction at the K-Pg boundary” needs references.\nReferences added (L351).\nReviewer #3:\nI laude the authors for their focus on regional diversity trends as a way to understand\n“global” diversification patterns of ammonoids during the Late Cretaceous. As Geerat\nVermeij stated in 2003, global diversity does not really mean what it implies given the well-\nestablished biases in the known fossil record at the global scale. This study is an interesting\nand important contribution to the understanding of a major group of animals as they\napproached their ultimate demise. The results indicate that despite parsing out global data\ninto distinct regions, ammonoid diversity was not waning during the Late Cretaceous as\nmany workers have suggested. Moreover, the results suggest that drivers of regional\ndiversity patterns are more likely driven by regional drivers, such as temperature and sea-\nlevel change.\nWe are pleased by reviewer’s positive response on our work and thank them for drawing our\nattention to Vermeij (2003), whose point of whether global trends are actually meaningful is\nvery pertinent to our work, so we have added this reference to the introduction (L65).\nOverall, this study highlights an exciting approach for assessing the evolutionary history of\nLate Cretaceous ammonoids and potential application to other major invertebrate groups in\nthe fossil record. The methods also seem sound to me, particularly in their approach for\ncontrolling sampling biases. I only have a few minor comments to make regarding the\ndiscussion and methodology, otherwise I believe this manuscript should be accepted once\nthese issues have been addressed.\nWe are pleased that the reviewer considers our approach to be robust and thank them for\ntheir additional comments, which we found to be very thought provoking.\nMinor comments:\nSea-level change has long been attributed to be a driver of both extinction and origination in\ndynamics in the fossil record. However, evidence suggests that sea-level change itself does\nnot meaningfully effect macroevolutionary dynamics, regardless of spatial scope (Holland\n2012, Paleobiology). Given that ammonoids are pelagic animals and even less likely to be\naffected by sea-level change, I fail to see the reason to include sea-level analysis in this\nstudy. Sea-level change will have a stronger relationship with the types of facies preserved,\nwith relatively deeper-water facies being deposited during times of sea level rise, which are\nmore likely to preserve ammonoid fossils than shallower water facies. Which leads into my\nnext comment.\nWe included sea level due to its effect on the extent of shallow marine shelf area. While\nammonites are pelagic, we nonetheless expect increased volume of pelagic habitat overlying\nproductive shallow marine shelf area to promote their diversification. Additionally, many\nheteromorph ammonitezs, which became increasingly diverse through the Late Cretaceous\nand are well represented in our dataset, were nektobenthic rather than pelagic and so\nshallow marine shelf area will be even more relevant to these taxa We chose to use sea\nlevel as a proxy rather than a direct measurement of shelf extent due to ongoing uncertainty\nbetween palaeogeographic reconstructions, which would only be compounded by the\nuncertainty in sea level change itself.\nHolland (2012) did find that that the relationship between sea level and shallow marine area\nwas complex and asserted that at local scales sea level will not be a simple predictor for\ndiversity under the species area effect. Holland (2012) also showed, however, that at\nregional scales that sea level and shallow marine area remain positively correlated (Fig. 4 in\npaper). Additionally, Holland (2012) did not actually conduct any formal analyses of empirical\ndiversity versus changing area under sea level and simply predicted diversity as a function of\nbenthic area using the Arrhenius equation. More recent studies of the empirical fossil record\ninstead provide support for long term sea level as having predictive power at larger spatial\nscales (e.g., Roberts and Mannion 2019, loc cit).\nUltimately, we assert that sea level remains a worthwhile predictor to test, but we accept that\nits relationship with diversity may be complex. We have added new material to the\ndiscussion to properly address this possibility (L241–250). At best, sea level had a genuine\neffect and is captured by our analyses. At worst, its effects could be ignored as excluding it\nfrom the PyRateMBD analyses would not impact the modelled effects of the other predictor\nvariables.\nAgain, while I appreciate the effort of the authors to consider the regional heterogeneity of\ndiversity dynamics, there is no mention regarding the nature of the stratigraphic record at the\nregional scale. While ammonoids are preserved across the majority of marine facies, they\nare most likely to be live and preserve in relatively deeper-water environments (80–400m;\nRitterbush 2014, J. of Zoology; Lukeneder et al. 2010, Earth Planet. Sci. Lett.). As such, the\nstructure of the stratigraphic record regionally would have a greater effect on diversity rather\nthan sampling biases. I wonder if stratigraphic packages regionally would potentially have a\n1st order control on diversity. I do realize, however, the difficulty in actually testing this, given\nthe poor state of facies/deposition environment data in the PBDB and lack of global\nstratigraphic packages in Macrostrat. Just a thought.\nIt has been previously determined that PyRate provides robust results even when confronted\nby strong preservation rate heterogeneity. As such, our regionalised results will account for\npreservation rate variation induced by regional stratigraphic architecture. We still agree,\nhowever, that there will be a trade off between biological signal versus stratigraphic noise in\nrelation to region size. As the reviewer points out, this would be challenging without proper\nstratigraphic data for each region, but it is still worth pointing this out in the manuscript. No\nwork has been done on how facies architecture might bias PyRate analyses and this will\ncertainly be important to tackle in future work, particularly given the reviewer’s pertinent point\nregarding sea level controlling facies types. We have added a paragraph to the discussion to\noutline these issues and believe that it significantly improves the presentation of our study\nlimitations (L324–340).\nLastly, while sea-surface temperatures might have mixed effects on diversity in this study, it\nis possible that temperature in conjunction with oxygen might be a better environmental\ncomponent to assess. While modern cephalopods are sensitive to temperature, they have\nbehavioural adaptations to compensate for increased temperatures, but not for temperature-\ndependent hypoxia.\nAs with the previous issue of the lack of stratigraphic databases to test for facies-driven\nbiases, our temperature and redox data are unfortunately insufficient to properly tackle this\nissue as they are one dimensional time series which summarise global environmental states,\nrather than providing the necessary information on how temperature and hypoxia spatially\nco-varied between regions. A couple of studies have convincingly demonstrated this effect in\nthe fossil record, but all have spatially resolved estimates of oxygen availability and\ntemperature from earth system models coupled with metabolic-limitation models of marine\ninvertebrates. Such analyses fall beyond the scope of this study, but they also provide\nexcellent examples of spatially explicit approaches to modelling diversification which can\nadditionally solve the issues we raise regarding spatially biased proxy records. We have\nadded reference to these model-based studies in our discussion to provide a clear pathway\nfor future work to pursue (L390–392), as well as using the empirical point of temperature-\ndependent hypoxia to present the potential limitations of our study more thoroughly (L251–\n258), although recent meta-analysis of cephalopod temperature tolerances found\nconsistently negative effects of warming across lineages and ontogenies (Borges et al.,\n2023, loc cit).\nOtherwise, great and interesting study!\nThanks again!"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I appreciated the authors' effort to address my comments.Looking forward to see the publish of this\nintriguing paper!", "reviewer_2": "Dear authors and editor,\nThanks for revising the manuscript.\nI have another careful read of this paper and believe it is ready for publication.\nBest regards,\nYadong Sun", "reviewer_3": "I originally reviewed this manuscript. The authors responded and addressed all my comments and\nhave substantially improved the manuscript. The authors responses were thoughtful to not only my\ncomments but also those made by the other reviewers as well. With the improvement and\nclarifications made in the new submission, I believe this manuscript warrants publication in Nature\nCommunications."}, "author_response": "Reviewer comments are given in regular text, our responses are in red.\nReviewer #1:\nI was really enjoyed to read this manuscript! The driver of diversification is one of the most\nattractive topics in both paleontologists and biologists. The authors adopted the state-of-the-\nart Bayesian method, PyRate, to recovery “global” and regional diversity dynamics of Late\nCretaceous ammonoids and found strong heterogeneity in diversity dynamics among\ndifferent regions and highly variable interactions between diversity and different biotic and\nabiotic factors in variable tempo-spatial scales. The manuscript is well written, and the\nfigures are excellent. I have only a few concerns/comments as follow:\nWe are delighted by the reviewer’s positive response to our work and thank them for their\nendorsement.\nI still can’t accept that species richness can lower than genus richness in some intervals in\nsome regions. This is not true based on our natural classification system. The authors could\ntreat, for instance “Haploceras sp.”, as a species to overcome this problem.\nWe agree that genus richness exceeding species richness is problematic under the nature of\nsystematics, but we do not posit that it is a true pattern. Instead we recognise it as an\nartefact of empirical palaeontological occurrence datasets that results from a deficit between\noccurrences identifiable at species level versus genus level. That this artefact is present in a\nclade as well sampled as ammonoids suggests that it is likely present elsewhere in the fossil\nrecord. The reality of incomplete sampling effects between taxonomic levels is an ugly one,\nbut we assert that this is an important limitation of empirical palaeontological occurrence\ndata which should be presented honestly, given the valid methodological concerns it raises.\nWe do not believe that adding placeholder species to a formal PyRate analysis would be a\njustifiable approach for three reasons:\n1. PyRate models preservation rates to help estimate true lineage durations, so the\ntaxonomic resolution of an occurrence will affect how it contributes to genus- versus species-\nlevel preservation rates. Dragging genus level occurrences to the species level would imprint\ngenus-level rates on the modelled species-level rates, distorting estimated species durations\nand their origination-extinction dynamics.\n2. While the purpose of placeholder species would be to meet the required richness ratio of\nat least one species per genus, it runs the risk of lumping genuinely separate (albeit\nunidentified) species into one composite record without any way of detecting this distortion.\n3. The logic behind assigning all genus-level occurrences of a taxon to one placeholder\nspecies means that that all occurrences which contribute information on the genus-level\nduration of that taxon should be treated in the same way. This would necessitate inclusion of\nspecies-level occurrences within their genus-wise placeholder species, despite already\nrepresenting their own lineages, and resulting in an analysis which is then largely collapsed\nback to the genus level. Conversely, inclusion of solely genus-level occurrences in\nplaceholder taxa would be an arbitrary choice to make species-level dynamics reflect genus\nlevel dynamics, even though they are ultimately different signals.\nA simpler, post hoc solution is for the reader to interpret the genus level curve as a lower-\nbound estimate on true species richness in instances where species richness falls below\ngenus richness. However, the high correlation values between genus and species richness\nin each region suggest that relative changes in the later still faithfully reflect changes in the\nformer, lending support to admittedly imperfect signal of species diversification we recover\nhere. Ultimately, we believe that it is better to present both sets of curves as they stand,\ngiven there is still no robust linkage in the literature between genus and species-level\ndynamics, and acknowledge their imperfections clearly in the discussion. To aid in this, we\nhave made explicit reference to the issue in the caption for fig. 3 (L1056–1058) and direct\nthe reader to the longer discussion of this issue in the text (L340–347).\nLines. 144-169: It’s a bit hard to read here. Why don’t you describe the diversity dynamics in\nthe eight regions in one paragraph. For instance: “Declines in genus and species richness\noccurred during the interval containing OAE2 in all four,”, At first glance, I was confused,\nwhat are the four regions? I can only find the answer (East Pacific, Tethys, WIS and Atlantic\nand Gulf) when I read the whole paragraph.\nWe agree that our original phrasing did not aid with manuscript clarity here. We now\ndescribe the dynamics in one paragraph as the reviewer suggests (L144-169). Our original\nstructure was designed to highlight how regions away from the Global North depart from the\n‘global’ diversity signal, so we have tried to retain this in the revised structure as we view it\nas a key result in our study.\nLine 353. “Middle” => “Early”.\nChange made (L401)\nLine 577-578: Do you mean five superfamily and twelve suborders? Please check.\nQuite right, these were misordered, change made (L625)\nI tried several times and can't access the code via the link provided above.\nWe apologise for this issue. We have revised our ESM in light of this review and have made\nthis new version available via the manuscript submission portal using Nature’s FigShare\nservice. We hope that this has addressed the problem, but if not then please get the editor to\ncontact us and we will endeavour to make the full ESM available in another way\nReviewer #2:\nSpatiotemporal variations of biodiversity impact our understanding of past evolution. This\nmanuscript focuses on the debate of what drives diversification and whether we can faithfully\nidentify these drivers from the fossil record. I do not have major concerns on the technical\npart of the paper. However, the author(s) soundly criticized past studies and, in some\nsentences, all palaeontologists. I am not sure if this is the best solution to promote this\npaper, but this is up to the authors.\nWe are pleased that the reviewer views our paper as technically sound. We agree that our\noriginal phrasing could be interpreted as a little damning of all palaeontologists, which is\ncertainly not our intention. We have modified our phrasing to reduce this possibility. We\ninstead view the need for spatially sensitive analysis as the ‘promoting’ point for the work.\nFrom an absolutely objective perspective, most researchers need a simple solution for many\napplications of biodiversity changes e.g., in palaeoenvironment studies. The approach\nprovided by the author(s), though impressive and better, is not simple, nor a solution to the\ngeneral problem. Nevertheless, this is a well-written manuscript that will add new insights\ninto our understanding of drivers of past extinction, recovery, and diversification.\nWe thank the reviewer for their positive perspective on the contributions of our manuscript to\nthe field. We agree that our analytical strategy is not the most straightforward to implement,\nbut this is largely because we aimed to analyse rates rather than standing diversity.\nMeasuring the latter in a spatially sensitive framework is much easier to achieve using\nstandard analytical tools and R packages and is a more achievable solution for future\nworkers. We view our take home message as the need for spatially sensitive analysis to\neffectively tackle questions of extinction, recovery and diversification in the fossil record.\nLine 69-71: True, but the words are rather harsh.\nWe have revised our phrasing to ‘some palaeontologists’ (L69). This is a more accurate\nreflection of the state of the field and implicitly acknowledges the excellent efforts of other\npalaeontologists to move away from global analyses towards spatially sensitive frameworks.\nLine 137: WIS first occurred here. But the abbreviation was only explained in line 421.\nWe have ensured that the abbreviation is now provided at the first instance of the complete\nphase (L95).\nLines 144-155: When reading this paragraph, the question popped up: Would the bias of\nsampling efforts completely rule out?\nOur analytical strategy is designed to minimise various sources of sampling bias as much as\npossible, but we agree that it is probably impossible to eliminate its effects entirely. We do\nnot address this point within the lines referenced by the reviewer to keep the results free\nfrom any discussion but have referred to the issue of residual sampling bias as a new\nsection later in the manuscript (L324–340). We believe that this inclusion greatly improves\nthe quality of our work and thank the reviewer for suggesting that we consider the issue of\nsampling in more detail\nLines 242-243: the sentence “regression-driven loss of habitat space is more commonly\nassociated with elevated extinction rates” needs references.\nWe have added references to support the impact of regression on extinction, but we have\nalso caveated that this relationship is not necessarily consistent across space, which helps\nto clarify why rising sea levels might instead coincide with elevated extinction rates (L243–\n250).\nLines 248-250: difficult to follow, explain why prevalence in the smaller sampling regions and\nabsence in the larger occur.\nThis was a remark on the prevalence of significant correlations between sp/ex and\nammonoid diversity in the smaller sampling regions, so we have added a reference to the\nfigure and reworked the phrasing for clarity. We have also added a comment with references\non a possible explanation for this finding based on the literature – that RQ effects are\noverwhelmed by CJ effects at larger spatial scales (L260–264).\nLines 250-252: regional ammonoid origination and extinction rates showed marked\nheterogeneity. How much difference can be inferred as marked heterogeneity? Please\nexplain.\nWe have removed ‘marked’ as this phrasing was vague with no clear basis for determining if\ndifference is marked or otherwise (L265).\nLines 288-293: Not very clear here. Taxonomic assignments are not affected by\nmeasurement location ….\nWe were a little uncertain of whether the reviewer is giving their interpretation of what we\nsaid in the manuscript in relation to the next point raised, or simply stating the above point\nplainly. We aimed to make the same point in the manuscript (L304–305), so we hope that we\nhave understood the reviewer here and that our response below addresses their intended\nconcern.\nLines 288-293 (cont): …. This seems to contradict the idea of spatial heterogeneity\nmentioned earlier. And the following sentence “fossil occurrences associated with taxonomic\nassignments display geographically biased sampling” doesn’t seem to make sense. Explain\nin more detail and give some references.\nWe disagree that independence of taxonomic assignment from spatial location contradicts\nthe idea of spatial heterogeneity. Our point is that the sampling of fossils as data points is\nspatially biased, but that this spatial bias will not have any impact on the taxonomic identities\nassigned to those data points (L304–309). Baculites sp. identified in one location will belong\nto the same taxonomic unit as Baculites sp. from another location, even if one of those\nlocations is disproportionately better sampled than the other. We have revised our phrasing\nto help ensure that this point is clear.\nLines 309-312: What then was the reason (driver) for ammonoid extinction?\nSentence and references added for the current stance in the literature on ammonoid\nextinction at the K-Pg (L354–358). We did this at a different point in the manuscript to the\nlines where the reviewer raises their comment to keep our conclusions as free as possible\nfrom inferences made elsewhere in the literature.\nLine 311-312: the sentence “the hypothesis of long-term global diversity decline prior to their\nfinal extinction at the K-Pg boundary” needs references.\nReferences added (L351).\nReviewer #3:\nI laude the authors for their focus on regional diversity trends as a way to understand\n“global” diversification patterns of ammonoids during the Late Cretaceous. As Geerat\nVermeij stated in 2003, global diversity does not really mean what it implies given the well-\nestablished biases in the known fossil record at the global scale. This study is an interesting\nand important contribution to the understanding of a major group of animals as they\napproached their ultimate demise. The results indicate that despite parsing out global data\ninto distinct regions, ammonoid diversity was not waning during the Late Cretaceous as\nmany workers have suggested. Moreover, the results suggest that drivers of regional\ndiversity patterns are more likely driven by regional drivers, such as temperature and sea-\nlevel change.\nWe are pleased by reviewer’s positive response on our work and thank them for drawing our\nattention to Vermeij (2003), whose point of whether global trends are actually meaningful is\nvery pertinent to our work, so we have added this reference to the introduction (L65).\nOverall, this study highlights an exciting approach for assessing the evolutionary history of\nLate Cretaceous ammonoids and potential application to other major invertebrate groups in\nthe fossil record. The methods also seem sound to me, particularly in their approach for\ncontrolling sampling biases. I only have a few minor comments to make regarding the\ndiscussion and methodology, otherwise I believe this manuscript should be accepted once\nthese issues have been addressed.\nWe are pleased that the reviewer considers our approach to be robust and thank them for\ntheir additional comments, which we found to be very thought provoking.\nMinor comments:\nSea-level change has long been attributed to be a driver of both extinction and origination in\ndynamics in the fossil record. However, evidence suggests that sea-level change itself does\nnot meaningfully effect macroevolutionary dynamics, regardless of spatial scope (Holland\n2012, Paleobiology). Given that ammonoids are pelagic animals and even less likely to be\naffected by sea-level change, I fail to see the reason to include sea-level analysis in this\nstudy. Sea-level change will have a stronger relationship with the types of facies preserved,\nwith relatively deeper-water facies being deposited during times of sea level rise, which are\nmore likely to preserve ammonoid fossils than shallower water facies. Which leads into my\nnext comment.\nWe included sea level due to its effect on the extent of shallow marine shelf area. While\nammonites are pelagic, we nonetheless expect increased volume of pelagic habitat overlying\nproductive shallow marine shelf area to promote their diversification. Additionally, many\nheteromorph ammonitezs, which became increasingly diverse through the Late Cretaceous\nand are well represented in our dataset, were nektobenthic rather than pelagic and so\nshallow marine shelf area will be even more relevant to these taxa We chose to use sea\nlevel as a proxy rather than a direct measurement of shelf extent due to ongoing uncertainty\nbetween palaeogeographic reconstructions, which would only be compounded by the\nuncertainty in sea level change itself.\nHolland (2012) did find that that the relationship between sea level and shallow marine area\nwas complex and asserted that at local scales sea level will not be a simple predictor for\ndiversity under the species area effect. Holland (2012) also showed, however, that at\nregional scales that sea level and shallow marine area remain positively correlated (Fig. 4 in\npaper). Additionally, Holland (2012) did not actually conduct any formal analyses of empirical\ndiversity versus changing area under sea level and simply predicted diversity as a function of\nbenthic area using the Arrhenius equation. More recent studies of the empirical fossil record\ninstead provide support for long term sea level as having predictive power at larger spatial\nscales (e.g., Roberts and Mannion 2019, loc cit).\nUltimately, we assert that sea level remains a worthwhile predictor to test, but we accept that\nits relationship with diversity may be complex. We have added new material to the\ndiscussion to properly address this possibility (L241–250). At best, sea level had a genuine\neffect and is captured by our analyses. At worst, its effects could be ignored as excluding it\nfrom the PyRateMBD analyses would not impact the modelled effects of the other predictor\nvariables.\nAgain, while I appreciate the effort of the authors to consider the regional heterogeneity of\ndiversity dynamics, there is no mention regarding the nature of the stratigraphic record at the\nregional scale. While ammonoids are preserved across the majority of marine facies, they\nare most likely to be live and preserve in relatively deeper-water environments (80–400m;\nRitterbush 2014, J. of Zoology; Lukeneder et al. 2010, Earth Planet. Sci. Lett.). As such, the\nstructure of the stratigraphic record regionally would have a greater effect on diversity rather\nthan sampling biases. I wonder if stratigraphic packages regionally would potentially have a\n1st order control on diversity. I do realize, however, the difficulty in actually testing this, given\nthe poor state of facies/deposition environment data in the PBDB and lack of global\nstratigraphic packages in Macrostrat. Just a thought.\nIt has been previously determined that PyRate provides robust results even when confronted\nby strong preservation rate heterogeneity. As such, our regionalised results will account for\npreservation rate variation induced by regional stratigraphic architecture. We still agree,\nhowever, that there will be a trade off between biological signal versus stratigraphic noise in\nrelation to region size. As the reviewer points out, this would be challenging without proper\nstratigraphic data for each region, but it is still worth pointing this out in the manuscript. No\nwork has been done on how facies architecture might bias PyRate analyses and this will\ncertainly be important to tackle in future work, particularly given the reviewer’s pertinent point\nregarding sea level controlling facies types. We have added a paragraph to the discussion to\noutline these issues and believe that it significantly improves the presentation of our study\nlimitations (L324–340).\nLastly, while sea-surface temperatures might have mixed effects on diversity in this study, it\nis possible that temperature in conjunction with oxygen might be a better environmental\ncomponent to assess. While modern cephalopods are sensitive to temperature, they have\nbehavioural adaptations to compensate for increased temperatures, but not for temperature-\ndependent hypoxia.\nAs with the previous issue of the lack of stratigraphic databases to test for facies-driven\nbiases, our temperature and redox data are unfortunately insufficient to properly tackle this\nissue as they are one dimensional time series which summarise global environmental states,\nrather than providing the necessary information on how temperature and hypoxia spatially\nco-varied between regions. A couple of studies have convincingly demonstrated this effect in\nthe fossil record, but all have spatially resolved estimates of oxygen availability and\ntemperature from earth system models coupled with metabolic-limitation models of marine\ninvertebrates. Such analyses fall beyond the scope of this study, but they also provide\nexcellent examples of spatially explicit approaches to modelling diversification which can\nadditionally solve the issues we raise regarding spatially biased proxy records. We have\nadded reference to these model-based studies in our discussion to provide a clear pathway\nfor future work to pursue (L390–392), as well as using the empirical point of temperature-\ndependent hypoxia to present the potential limitations of our study more thoroughly (L251–\n258), although recent meta-analysis of cephalopod temperature tolerances found\nconsistently negative effects of warming across lineages and ontogenies (Borges et al.,\n2023, loc cit).\nOtherwise, great and interesting study!\nThanks again!"}]

Decision letter