Impact of global climate cooling on Ordovician marine biodiversity

Domain: earth_science
Published: 2023-10-10
Editorial outcome: minor_revision

Abstract
Global cooling has been proposed as a driver of the Great Ordovician Biodiversification Event, the largest radiation of Phanerozoic marine animal Life. Yet, mechanistic understanding of the underlying pathways is lacking and other possible causes are debated. Here we couple a global climate model with a macroecological model to reconstruct global biodiversity patterns during the Ordovician. In our simulations, an inverted latitudinal biodiversity gradient characterizes the late Cambrian and Early Ordovician when climate was much warmer than today. During the Mid-Late Ordovician, climate cooling simultaneously permits the development of a modern latitudinal biodiversity gradient and an increase in global biodiversity. This increase is a consequence of the ecophysiological limitations to marine Life and is robust to uncertainties in both proxy-derived temperature reconstructions and organism physiology. First-order model-data agreement suggests that the most conspicuous rise in biodiversity over Earth’s history – the Great Ordovician Biodiversification Event – was primarily driven by global cooling.

Introduction
Following the appearance of most animal phyla during the “Cambrian Explosion”1, the Ordovician reflects their rapid diversification at lower taxonomic levels during the “Great Ordovician Biodiversification Event” (GOBE)2. This radiation constitutes the most prominent increase in the biodiversity of marine organisms during the entire Phanerozoic3,4. It was already identified in the seminal work of Sepkoski5 and stands in most recent compilations3,4, although the timing of its onset is strongly dependent on the studied taxonomic group and geographical region considered6.

The causes of the GOBE remain largely unconstrained and several hypotheses are debated. A first hypothesis invokes a decrease in tropical ocean temperatures from levels deleterious to marine organisms7 (>40°C) in the Cambrian to modern-like values (ca. 30°C) in the Early and Middle Ordovician8. This hypothesis is supported by updated temperature proxy compilations9,10. A second hypothesis proposes that the fragmentation of supercontinent Rodinia and related eustatic sea-level rise may have created new ecological niches2,11. A third hypothesis, based on the analysis of O-dependent fractionation of carbon isotopes during photosynthesis, suggests a long-term oxygenation of the early Paleozoic atmosphere12. A fourth hypothesis suggests that frequent impacts on Earth of kilometer-sized asteroids may have created new ecological niches13. A last hypothesis proposes that the observed increase in phytoplankton biodiversity may have modified trophic structures and fostered the diversification of plankton-feeding groups14.

Recently, the extraterrestrial hypothesis has been rejected based on updated dating of geological strata15. In addition, evidence has accumulated that atmospheric oxygen concentrations may have already been relatively elevated during the Cambrian and Early Ordovician16–18, and proxies for the global extent of anoxia testify of stable ocean redox during the main phase of the GOBE19, together ruling out oxygenation as the main driver of biodiversification. Question remains, whether the GOBE may result from extrinsic (i.e., paleogeographical and/or climatic) or intrinsic (i.e., biological) drivers.

Here we use a coupled paleoclimatic and macroecological model20,21 to evaluate the individual contributions of extrinsic mechanisms, namely the long-term Ordovician cooling trend8–10 and paleogeographical evolution11,22, to the trajectory of marine biodiversity during the Ordovician (between 490 to 430 Ma). In contrast with previous attempts based on temporal correlations, our approach allows us to quantify processes and establish causal relationships.

Results
Global cooling triggers Ordovician biodiversification
Ocean temperatures were simulated using the Fast Ocean Atmosphere Model (FOAM), a mixed-resolution ocean-atmosphere general circulation model23 (Methods). We conducted one simulation every 10 Myrs from 490 Ma to 430 Ma. Solar luminosity follows the model of stellar physics of ref. 24 and, in our main experiments, we varied atmospheric pCO2 to reproduce the long-term Ordovician cooling trend of the oxygen isotope compilation of ref. 10. As a result, tropical sea-surface temperatures decrease from ca. 45°C at 490 Ma to ca. 30°C at 430 Ma. Importantly, FOAM has been shown to satisfactorily capture the latitudinal seawater temperature gradient during both Paleozoic icehouse periods (Ordovician glaciation25) and Mesozoic thermal maxima (Oceanic Anoxic Event 2, see Peer Review file published with ref. 26).

Biodiversity was simulated using a macroecological model based on the interaction between many modeled pseudo-species and their environment, which captures biodiversity patterns well in the modern ocean21 (Methods). Here we forced it with the ocean temperatures simulated with FOAM (Methods). A similar coupling has previously been shown to successfully capture the spatial-temporal biodiversity patterns of early Paleozoic acritarchs20. In the model, thousands of pseudo-species were generated, each one occupying a unique thermal niche sensu Hutchinson27 based on the principle of exclusive competition28. Although the number of ecological niches can be decoupled from the number of species under certain circumstances in reality29, we note that this and similar modeling approaches have previously been proven to satisfactorily simulate biodiversity at the global scale amongst groups of organisms of various complexity, both on land and in the ocean30,31. In each model grid point, biodiversity is defined as the number of pseudo-species whose thermal niche overlaps the mean annual ocean temperature. The lowermost (tmin = –1.8°C) and uppermost (tmax = 44°C) temperature bounds of the niches were chosen to provide the best simulation of biodiversity spatial patterns in the modern ocean21. We further note that the lower bound is physically imposed as the temperature at which seawater freezes, while the upper bound is close to the value of 41°C derived from the analysis of the temperature limits of tropical marine ectotherm species of ref. 7. We focused our analysis on the continental shelves and excluded both the deep ocean (here defined as all grid points deeper than –250 m) and polar (>60° latitude) regions since they are not well represented in paleontological databases20,32,33.

Analysis of the spatial patterns reveals that the maximum of simulated biodiversity (or equivalently, maximum number of niches available) shifts from high latitudes (60° S) between 490 Ma and 470 Ma, to middle latitudes at 460 Ma, and finally to low latitudes (from 450 Ma onward) (Fig. 1). This spatial-temporal pattern arises from the interaction between the modeled pseudo-species and their environment, with lethally high ocean temperatures making the low-latitude shelves hostile to marine life in the late Cambrian and gradually more suitable as global climate cools. This result, i.e., that global cooling triggers biodiversification in the model, arises in part from the niche-temperature interaction and a thermal mid-domain effect34,35; ecological niches being defined between minimum and maximum temperature bounds, a maximum number of niches overlaps intermediate temperatures, resulting in maximum biodiversity around 23°C (ref. 35). In the absence of better constraints on the ecophysiology of marine organisms in the deep past, our approach necessarily relies on the assumption that this biodiversity–temperature relationship did not change over geological time. While such assumption may seem limiting, we note that the biodiversity–temperature relationship emerging from our model is supported by the mechanistic model of ref. 36 that demonstrates that the diversity-temperature relationship is expected, on grounds of thermodynamics and physiology principles, to peak at moderate temperatures in the range of 15–25°C and decrease outside of this range.

The increase in marine biodiversity simulated during the Ordovician stands when the cooling scenario is varied in the model, although the specific shape of the temporal trend in biodiversity depends on the details of the temperature reconstruction9,37 (Supplementary Figs. 1 and 2). The possible underestimation of the magnitude of Ordovician cooling arising from the heterogeneous spatial distribution of temperature proxies38 suggests that the simulated increase in biodiversity may be underestimated as well. Our model does not capture the Late Ordovician Mass Extinction (ca. 445 to 443 Ma), which was caused by a combination of factors that are not accounted for here, including Hirnantian sudden cooling, sea-level fall and possibly anoxia25,39,40.

Model results were compared with three different fossil global biodiversity curves, mainly composed of brachiopods, trilobites and conodonts (Fig. 2). Differences between these curves arise from the variable temporal resolution, with the seminal curve of Sepkoski et al. (ref. 5) exhibiting the lowest resolution (>5 Myrs). Differences also reflect the use of different datasets. Rasmussen et al. (ref. 41) based their compilation on the Paleobiology Database (PBDB), which mostly includes data from Laurentia and Baltica (>50% in 2019), while the Geobiodiversity Database3 (GBDB) mainly includes Chinese data. In our simulations, species richness monotonically increases in response to global climate cooling from the late Cambrian (490 Ma) to the Katian (Late Ordovician, 450 Ma) and subsequently stabilizes between 450 Ma and the early Silurian (430 Ma), in first-order agreement with the reference biodiversity curves (Fig. 2). Biological evolution may have played a role during the GOBE as well14,42, but it is not necessary to invoke this process to simulate a first-order rise in marine biodiversity during the Ordovician in our model.

In order to facilitate the attribution of causes and quantify the respective contributions of global climate change and paleogeographical evolution to model biodiversification, we conducted additional simulations characterized by a constant tropical sea-surface temperature of ca. 30°C (Fig. 2). These simulations constitute a sensitivity test to the continental configuration. They show that continental drift alone triggers no substantial increase in marine biodiversity and confirm that global climate cooling is the driver of the simulated marine biodiversification. The same conclusions are reached when conducting simulations under a warmer global climate with tropical sea-surface temperatures of ca. 40°C (Supplementary Fig. 3).

Impact of the thermal upper limit of model pseudo-species
Our main experiments are based on contemporaneous thermal limits for marine Life21. This possibly constitutes an actualistic bias, since Ordovician organisms may have been adapted to the warmer early Paleozoic climate. To determine to what extent model biodiversification is affected by the thermal upper limit tmax, we conducted additional simulations for alternative values of tmax. We quantified the degree of similarity between the biodiversity spatial patterns simulated under warm and cool climates, characterized by tropical sea-surface temperatures of ca. 40°C and 30°C. These values were chosen to represent the end-member (respectively warm and cool) climatic states of the long-term Ordovician cooling trend. The closer the results obtained under these contrasting global climatic states, the weakest the impact of global climate cooling on biodiversity. Simulations were performed using the continental configuration for 460 Ma, which represents an intermediate state between 490 Ma and 430 Ma. This sensitivity analysis is presented in a Taylor diagram (Fig. 3). Correlation is the highest for tmax = 65°C. Standard deviation is the closest to the reference simulation for tmax = 78°C. The root mean square deviation is the lowest for values of tmax ranging between 68°C and 78°C. The correlation decreases very fast when tmax decreases, and only the simulation using tmax = 65°C shows a correlation exceeding 0.95. Altogether, these results demonstrate that a thermal upper limit exceeding 65°C is necessary to suppress the impact of global climate cooling on biodiversity. This value is >20°C above the (44°C) modern limit used in the main experiments.

Physiological adaptation of the Paleozoic Life to such extent is unlikely. The compilation of ref. 43 shows a very limited change in the thermal tolerance of marine organisms over the last 500 million years. In addition, several lines of evidence suggest that the atmospheric oxygen concentration was, if anything16,17, lower than today during the early Paleozoic44,45. Due to the increase in organism metabolism (thus dissolved oxygen consumption) with increasing temperatures, a reduced ocean oxygenation would lower the thermal upper limit for marine ectotherms compared to modern33. It is therefore likely that tmax during the early Paleozoic would have been lower, rather than higher, than modern. Therefore, the increase in biodiversity simulated in response to global cooling stands as a robust model result that is poorly dependent on (reasonable) changes in organism physiology.

The temperature-biodiversity relationship also depends on the distribution of the niches in the temperature dimension. In our simulations, the niches are uniformly distributed along the temperature axis, between minimum and maximum temperatures. For a same tmax, however, an asymmetric distribution of the niches skewed towards the high temperatures would result in more species adapted to warmer climates. It could be argued that such adaptations of marine Life to warmer temperatures may have developed in early Paleozoic oceans, significantly reducing the impact of global cooling on marine biodiversity in our simulations. However, such alteration of the temperature-biodiversity relationship would increase the number, packing and ultimately overlapping of niches at high temperatures. It would conflict with the fundamental principle of competitive exclusion28. Therefore, we consider such changes in the temperature-biodiversity relationship as unlikely.

Establishment of a modern latitudinal diversity gradient
The latitudinal biodiversity gradient (LBG) – the increase in biodiversity from the poles to the low latitudes – is a pervasive macroecological pattern of life on Earth21,46. The distribution of fossil occurrences shows that this pattern has changed over long periods of time20,47,48. Our simulations suggest that, following the appearance of most animal phyla during the Cambrian Explosion1, a modern-like LBG first established during the GOBE in response to global climate cooling. During the late Cambrian and early Ordovician (490 to 470 Ma), the maximum in modeled biodiversity is located over mid to high latitudes, contrasting with modern LBGs, while the lower latitudes are too warm to host diverse and abundant marine life (Fig. 4). In response to global climate cooling, the maximum in biodiversity gradually shifts to lower latitudes as large habitats gradually become livable for the modeled pseudo-species. The maximum in biodiversity reaches ca. 40° S at 460 Ma and covers the 40° S–20° N latitudinal band from 450 Ma onwards. The evolution of the LBG simulated during the Ordovician is robust when the model is run over the global ocean (Supplementary Fig. 4). Therefore, results are not overly dependent on the considered continental reconstruction and position of the shallow-water shelves (see Fig. 1).

Although sampling bias may distort the LBGs reconstructed based on paleontological data and obscure model-data comparison46, simulated spatial-temporal patterns of biodiversity are supported by the LBG of Cambrian and Ordovician acritarchs, which show a minimum in biodiversity along the equator at 490 Ma and a shift of maximum biodiversity from ~60° S at 480 Ma to the low latitudes during the Late Ordovician20. Our results are also supported by the LBG independently reconstructed during the Ordovician by Kröger et al. (ref. 32) based on the Paleobiology Database4, which shows a transition from a bimodal Early Ordovician LBG to a unimodal Late Ordovician LBG. The authors note that this evolution results from the large increase in low-latitude biodiversity during the GOBE, in agreement with the results of our coupled climate-macroecological model (Figs. 1 and 4).

Our results, obtained with a coupled climate-macroecological model featuring minimal biological assumptions, shed new light on the mechanisms that allowed the establishment of a modern-like LBG for the first time in Earth’s history. They suggest that global climate change, rather than the two other usual suspects, paleogeography or history48, was decisive. They also reaffirm the contribution of a cool global climatic state to the development of a steep LBG featuring a maximum in biodiversity at the low latitudes in line with paleontological data of the Permian-Triassic transition47.

Regional biodiversification patterns
Our spatially-resolved simulations allow biodiversification signals to be compared in different paleocontinents. For comparisons, we used the brachiopod database of ref. 6. (Fig. 5). Our model captures the increase in fossil brachiopod biodiversity in Laurentia, Siberia and South China (Fig. 5a–c), although with varying degrees of model-data agreement. In Laurentia and Siberia, simulated biodiversity rises too early compared to data. In South China, the general timing of biodiversity increase matches the temporal trends reconstructed based on data from the PBDB between 490 Ma and 450 Ma6. From 440 Ma to 430 Ma, model biodiversity stabilizes while paleontological data show a large drop in biodiversity. This discrepancy in South China probably reflects the even larger uncertainties in the latitudinal position of Baltica due to the lower number of paleomagnetic poles in use (3 in ref. 50. vs. 24 in the reconstructions used in our main experiments51). Model-data disagreement for Baltica may also reflect the inability of our model to capture finer, regional-scale environmental changes contributing to biodiversification. Indeed, quantitative analysis of brachiopod occurrences from the eastern Baltic paleobasin suggests that carbonate shelf development played a key role in regional diversification by enhancing substrate heterogeneity52, in line with the vision that reefs are biodiversity hotspots53.

Extrinsic vs. intrinsic drivers of GOBE
Our results have implications for our understanding of the mechanisms that drove the GOBE. They highlight the decisive role of global climate change, although large uncertainties remain in the reconstruction of early Paleozoic ocean temperatures. Indeed, the low oxygen isotope values typifying the early Paleozoic are sometimes interpreted as reflecting very high ocean temperatures (>40°C8,9 at low latitudes, possibly exceeding 50°C10,37) under the assumption of a constant oxygen seawater composition during the Phanerozoic – an assumption supported by coupled oxygen isotope and clumped isotope measurements that show no long-term trend in seawater δ18O over the last 500 Myrs54. Other studies argue that the long-term trend in δ18O represents variations in the oxygen isotope composition of Phanerozoic oceans55,56. Correction for this temporal drift results in early Paleozoic temperatures that are in the same order of magnitude as Mesozoic and early Cenozoic greenhouse climates55 or even modern56. This second view is supported by the climate model-based interpretation of the distribution of the lithological indicators of climate of ref. 57. Our coupled climate-biodiversity modeling shows that a long-term cooling from very high temperature in the Cambrian to modern-like levels in the Ordovician provides a conservative, minimal scenario allowing the GOBE to take place in the model, and a modern-like LBG to emerge in agreement with paleontological evidence, even without invoking additional, biological drivers. We do not intend to mean that evolutionary innovation and new ecospace occupation did not contribute to biodiversification2, and actually suggest that biotic factors may have played an important role at the regional scale in Baltica, but the fact that a major radiation takes place globally in the model without the need to invoke such ecological changes suggests that major global climate cooling may have played a leading role.

Methods
Climatic simulations
Model description. FOAM version 1.523 is a mixed-resolution ocean-atmosphere general circulation model. The atmosphere is a parallelized version of the National Center for Atmospheric Research’s (NCAR) Community Climate Model 2 (CCM2), with radiative and hydrologic physics upgraded after CCM3 version 3.258. It was run using a R15 spectral resolution (4.5°×7.5°) and 18 levels in the vertical dimension. The oceanic component is the Ocean Model version 3 (OM3), a 24-level z-coordinate ocean general circulation model benefiting of a higher resolution than the atmosphere (1.4°×2.8°). The simulation of sea ice uses the thermodynamic component of the CSM1.4 sea-ice model, which is based on the Semtner 3-layer thermodynamic snow/ice model59. The coupled model, FOAM, is well designed for paleoclimate studies. It has no flux correction and its quick turnaround time allows for long millennium-scale integrations. The higher execution speed compared to higher-resolution Earth system models, e.g. of the Coupled Model Intercomparison Project phase 6 generation, permits exploring the parameter space considered here (more than twenty 2000-year ocean-atmosphere simulations). FOAM has been widely used in paleoclimate studies26,55,60, including in the Ordovician25,61–64.

Description of the experiments. We ran equilibrium FOAM simulations every 10 Myrs between 490 Ma and 430 Ma, both included (for a total of 7 time slices). In our main experiments, we used the Ordovician continental reconstructions of ref. 22. and a rocky desert vegetation type in every land grid point. Solar luminosity was calculated for each time slice after the stellar physics of ref. 24, increasing from 1313.24 Wm−2 at 490 Ma to 1319.7 Wm−2 at 430 Ma. In terms of astronomical forcing, we used a circular (null eccentricity) Earth’s orbit and the obliquity was set to 22°. These parameters provide an equal, annual insolation for both hemispheres with minimal seasonal contrast. Simulations were run for a wide range of atmospheric pCO2 values permitting to capture the magnitude of Ordovician cooling (Supplementary Fig. 1). In our main experiments, we reproduced the long-term Ordovician climatic evolution of ref. 10. by imposing atmospheric pCO2 values of 192 PAL (preindustrial atmospheric CO2 level; 1 PAL = 280 ppm CO2) at 490 Ma, 96 PAL at 480 Ma, 48 PAL at 470 Ma, 24 PAL at 460 Ma, 8 PAL at 450 Ma, 6.1 PAL at 440 Ma and 4 PAL at 430 Ma. While pCO2 values used for 460–430 Ma are well in line with values reconstructed for the early Paleozoic using long-term carbon cycle models63,65,66, pCO2 values used for older time slices may seem particularly elevated. Yet, we note that CO2 is the only greenhouse gas varied in our model while variations in the concentration of other gases, and notably the reduced atmospheric oxygen concentration67, may have also contributed to global climate warming. Therefore, the pCO2 forcing used should be considered as a pCO2-equivalent forcing rather than strictly pCO2. In addition, the temperature increase resulting from CO2 forcing (i.e., climatic sensitivity) is highly model-dependent. For instance, even most recent general circulation models of the Coupled Model Intercomparison Project phase 6 (CMIP6) generation lead to a large spread of values ranging from 1.8 to 5.6 K per pCO2 doubling68. The climatic sensitivity of FOAM is ca. 3.5 K26; lower pCO2 values would permit simulating identical early Paleozoic temperatures using a model with a higher climatic sensitivity. Finally, the exceptionally high climatic sensitivity and/or pCO2 increase required to simulate early Paleozoic climate is not unexpected. It was previously noted by ref. 37 and it is not our objective here to provide climatic mechanisms. Instead, we focus on the consequence of resulting ocean temperatures for marine habitability and the temporal trajectory of marine biodiversity in the deep geological past.

Simulations were initialized with a spatially uniform field of ocean salinity (35‰). In order to ensure intense convection and rapid deep-ocean equilibrium, simulations were initialized with a warm ocean. Simulations were run for 2000 years until deep-ocean thermal equilibrium (median ocean temperature drift of 0.0035°C yr−1 over the last 100 years at 3700 m depth). Mean annual results from the last 50 years were used to derive the temperature forcing fields used in the biodiversity model.

Macroecological modeling
Model description. We modeled species richness in the Ordovician from 490 Ma to 430 Ma with a 10-Myr interval by means of a model based on the MacroEcological Theory on the Arrangement of Life (METAL)35,69; the version of the model used here was Species Niche and Climate Integration (SNCI) model35. The METAL theory is described in https://biodiversite.macroecologie.climat.cnrs.fr. We created pseudo-species (i.e. virtual species), each having a unique thermal niche with distinct degrees of thermophily and eurythermy. No hypothesis was made about species distribution. We allowed pseudo-species to settle in any given region of the global ocean so long as they could withstand the mean annual upper-ocean temperatures averaged of the first 100 meters of the water column. Thermal niches within the pool ranged from tmin = –1.8°C, the temperature at which seawater freezes, to tmax = 44°C, the maximum temperature under which a metazoan can reproduce70. The breadth of the thermal niche varied from n=1°C (stenothermic species) to n=45.8°C (universal eurythermic species). In total, 101,397 thermal niches were generated, as a function of niche breadth and the increment between niche breadth. To model species richness, we adapted the procedure described in Zacaï et al. (ref. 20). A total of 10,000 pseudo-species were randomly selected (without replacement) among the 101,397 thermal niches available. The random extraction was repeated 100 times. We used a lower number of repetitions compared to previous work (i.e., 100 vs. 1000 in ref. 20) to reduce computational cost but results were poorly impacted by the number of repetitions.

Description of the experiments. We ran the SNCI model on temperatures simulated with the ocean-atmosphere general circulation model FOAM reproducing, in our main experiments, the temperature reconstruction of Song et al. (ref. 10) (Supplementary Fig. 1a). In order to best capture the environmental conditions characterizing Ordovician shallow-water shelves, we used the upper-ocean temperatures averaged over the first 100 m. We mapped the results from 490 Ma to 430 Ma by calculating an average biodiversity based on the results of all members of the ensemble of 100 model realizations (Fig. 1). We then derived a global biodiversity estimate by calculating the median of all biodiversity values for each map between 490 Ma and 430 Ma. The median was chosen to avoid any strong alteration of the central tendency by potential outliers. We calculated the standard deviation of simulated biodiversity in the spatial domain to represent the likely range of values found in each time slice and represented the result in the form of an envelope (Fig. 2).

Statistics & reproducibility
No data was collected for this study.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The paper “Cooling Oceans Did Trigger Ordovician Biodiversification” is a well written, concise report of the results of a modelling of Ordovician marine biodiversity based on a combination of climate models, explicit paleogeographical reconstructions, and spatially explicit sea water temperature-species richness models. The results of the modelling largely reproduce global temporal diversity curves and Ordovician latitudinal biodiversity gradients, which were previously and independently estimated from empirical data.\nThe authors interpret their results as support for the hypothesis that global cooling triggered the biodiversification during the Ordovician. The presented results support their interpretation.\nThis is an interesting, novel, and significant approach to test among the hypotheses about the mechanism behind the Ordovician Biodiversification. It should be published in Nature Communications, not only because it discusses the processes that led to this specific radiation in the marine realm, but because it is a contribution toward the wider discussions about the geological history and evolution of the latitudinal diversity gradient.\nThat said, I would suggest a few improvements and I found some minor flaws in the text which are listed below. I would suggest, to refer more explicitly to the wider discussions about the nature and evolution of the latitudinal diversity gradient (see listed references below). Also, I would suggest a more informed discussion about potential processes that led to discrepancies between the model results and empirical data.\np. 4 line 83: Here it reads: “This spatial-temporal pattern arises from the niche-environment interaction”. Can you explain what you specifically mean with this? What interacts with the environment. I guess it is not the niche, because the niche is an abstract concept, if anything interacts with the environment, then the organism. So, in your model: what interacts?\np. 5 lines 99-101. It is said that the simulations result in a monotone richness increase in response to climate cooling. In the context of your model, can you say a few words about the causes of the increase in your model? Do more nodes simply correspond to livable thermal niches?\np. 7. Lines 143-144: It is stated that “The distribution of fossil occurrences shows that this pattern has changed over long periods of time.” With a single reference that only gives an example from Ordovician microfossils. I think, there are more appropriate references, e.g.,:\nBadgley, C., T.M. Smiley, R. Terry, E.B. Davis, L.R.G. Desantis, D.L. Fox, S.S.B. Hopkins, T. Jezkova, M.D. Matocq, N. Matzke, J.L. Mcguire, A. Mulch, B.R. Riddle, V.L. Roth, J.X. Samuels, C.A.E. Strömberg, and B.J. Yanites. 2017. Biodiversity and Topographic Complexity: Modern and Geohistorical Perspectives. Trends in Ecology & Evolution, 32: 211–226.\nCrame, J.A. 2002. Evolution of taxonomic diversity gradients in the marine realm: a comparison of Late Jurassic and Recent bivalve faunas. Paleobiology, 28: 184–207.\nJablonski, D., S. Huang, K. Roy, and J.W. Valentine. 2017. Shaping the latitudinal diversity gradient: new perspectives from a synthesis of paleobiology and biogeography. The American Naturalist, 189: 1–12.\nMannion, P. D., Upchurch, P., Benson, R. B., & Goswami, A. (2014). The latitudinal biodiversity gradient through deep time. Trends in ecology & evolution, 29(1), 42-50.\nPowell, M. 2009. The Latitudinal Diversity Gradient of Brachiopods over the Past 530 Million Years. The Journal of Geology, 117: 585–594.\nSaupe, E.E., C.E. Myers, A. Townsend Peterson, J. Soberón, J. Singarayer, P. Valdes, and H. Qiao. 2019. Spatio-temporal climate change contributes to latitudinal diversity gradients. Nature ecology & evolution, 3: 1419–1429.\nSong, H., S. Huang, E. Jia, X. Dai, P.B. Wignall, and A.M. Dunhill. 2020. Flat latitudinal diversity gradient caused by the Permian–Triassic mass extinction. Proceedings of the National Academy of Sciences, 117: 17578–17583.\nIt would be good to discuss the results presented in the manuscript with the main topics discussed in these reviews: Phanerozoic temperature change vs paleogeography. What is the contribution of the current manuscript to these discussions?\np. 8. Line 170. There is written that the timing of the biodiversity increase “matches the paleontological database 490 Ma and 450 Ma.” What does this mean? Do you mean \"matches that of results from estimations based on data from the PBDB\". If yes: reference is needed. Also: What means the “Afterwards” in the same line. After what?\np. 9 line 184-187: Here I would refer to Penny, A. M., Hints, O., & Kröger, B. (2021). Carbonate shelf development and early Paleozoic benthic diversity in Baltica: a hierarchical diversity partitioning approach using brachiopod data. Paleobiology, 1-21. doi:10.1017/pab.2021.3.\nWhy: Penny et al. (2021) give a hypothesis why the diversity rose so extremely in Baltica: “The results are consistent with carbonate shelf development as a driver of Ordovician regional diversification, both by enhancing broadscale differentiation between shallow- and deep-marine environments and by generating heterogeneous carbonate environments that allowed increasing numbers of brachiopod genera to coexist. However, temporal turnover also contributed significantly to apparent regional diversity, particularly in the Middle–Late Ordovician.”\nIt would be interesting to discuss if this process also potentially contributed to the discrepancies you see between the model results and estimations from empirical data in Laurentia and Siberia. In this context see also:\nClose, R. A., Benson, R. B., Saupe, E. E., Clapham, M. E., & Butler, R. J. (2020). The spatial structure of Phanerozoic marine animal diversity. Science, 368(6489), 420-424.\np. 9. Line 200: Reference needed.\np. 10. Lines 206-207. I am not sure if I understand that sentence. Does this suggest that the authors hypothesize that biological drivers were absent (but see above). Or that they are not necessary to explain the diversification in the light of the model (but see lines 173, 184). I would suggest to improve the final statement.\np. 17 line 377: is written: The procedure was repeated 100 times, ten times less than the previous procedure”. It is not clear what is the previous procedure. Please improve.", "reviewer_2": "In the submitted manuscript \"Cooling Oceans Did Trigger Ordovician Biodiversification\" Ontiveros et al. use simulations to reconstruct potential biodiversity patterns during the Great Ordovician Biodiversification Event (GOBE). A Fast Ocean Atmosphere Model was employed first to reconstruct spatially explicit patterns of paleotemperature, which then defined habitats that were allowed to be inhabited by virtual species. The results of the authors suggest that the decreasing temperature creates the modern-like latitudinal diversity gradient that did not exist previously - increasing overall diversity, which matches the GOBE observed in the fossil record.\nThese results are excellent in highlighting the patterns under the conditions and processes that the authors assumed, and they definitely need to be published and shared with the scientific community. Yet, these results imply specific assumptions (see below), that need to be stated explicitly and might indicate that the primary interpretation of the results (i.e. climate changes caused the GOBE) needs to be toned down (e.g. climate alone could cause the GOBE). Although there are reasons suggesting that the assumptions behind the results hold - in my opinion - they do not conclusively falsify other potential drivers of the diversification event (e.g. evolutionary innovations). Nevertheless, the results of the study represent a significant advancement and (with some additional analyses and clarifications) will be important contributions to the study of the early Phanerozoic life and the GOBE.\nThe results of the study are based on two important aspects of the study: 1) the distribution of virtual species in the temperature spectrum and 2) the overall decrease of global temperature, which was tracked by the authors using the CO2 concentration in the climate modelling.\n1) Virtual species distributions\nAs the authors noted, the results are entirely dependent on the distribution of species (niches) in the temperature spectrum. Based on the description and the supplied code (Step5_niches.m), the method used to generate virtual species has an implicit assumption that all possible species niches are equally likely to exist. In this model, the expected number of species niches is distributed in the temperature spectrum similar to a Beta distribution, with most species occurring around the mid domain of ~23°C, with richness declining towards the colder and warmer end symmetrically. The result is that at temperatures close to the maximum limit the expected richness will be lower, which is projected to geographic space - thus creating the pattern of lower diversity in the tropical Cambrian.\nIn the modern oceans, however, the spatial diversity of taxa varies extensively over latitude and temperature (especially across taxa), and, as mentioned by Beaugrand et al (2018, cited in Zacaï et al 2021, on which the authors rely) a considerable number of niches can be vacant/unsaturated - decoupling diversity from the number of potential niches that are possible in an area. Species also do not occupy their entire thermal niches, making this connection even looser. The method therefore illustrates the number of potential niches, and I would state this explicitly for clarity and unambiguity (e.g. the cooling climate makes a higher number of niches available and therefore could cause the increase in diversity).\nAlthough showing the results based on the assumed richness ~ temperature relationship is very important, they are, nevertheless, based on hard assumptions. The authors also assumed that that this relationship did not change over time through evolution, but it is easy to imagine a scenario where more species adapt to the high temperatures (i.e. more colder niches remain unfilled) when such conditions are more widely distributed. It is also important to note that model was assessed for fishes, cetaceceans, pinnipeds and planktonic organisms - and its applicability to well-fossilizable, benthic organisms that make up the majority of the fossil record might be different than for planktonic lifeforms (such as in Zacaï et al 2021).\nThe authors can consider assessing the sensitivity of the results to differences in the distribution of species in the temperature dimension, not only to the changes of the maximum limit, but to the shape and balance of the niche-distribution as well (e.g. maintaining the same t-max but making the distribution of the number of niches at a temperature more assymmetric, i.e. shift the mode). In addition to the through-time comparison of simulated and observed patterns of diversity, the analyses of the latitudinal diversity gradient in each time interval (as in Zacaï et al 2021) might help with the assessment of the model's adequacy and make its applicability to the scenarios more credible. For the sake of clarity and more accurate interpretation by the reader, the authors can use \"number of sampled niches\" or something similar throughout the manuscript to better illustrate the difference between actual diversity and what the models imply. The virtual species (i.e. niches) are also not limited in dispersion (e.g in Hagen et al. 2021), which is necessary for the continent configuration to have a serious effect on diversity patterns via allopatric speciation. This aspect (i.e. dispersal limitation) can also be added to the model if the effect of continent configuration is to be ruled out as a causal factor in the diversification.\n2) Temperature and climate modelling\nGiven the sensitivity tests run by the authors (stable climate, changing continental configuration), the long-term increase in diversity and the change of the LDG shape is entirely dependent on the temperature decrease. Since I am neither a climate modeller nor a geochemist, I cannot assess the robustness of the paleoclimate reconstructions themselves. The authors also tried an alternative cooling scenario (Extended Data Figure 2), in which the presence of the diversity increase seems to be dependent on a single time point (480Ma) and the pseudodiversity at 490Ma is comparable to those later in the Ordovician. This indicates that the sensitivity to this parameter might need to be tested more, maybe by considering additional reconstructions (see Grossman and Joachimski, 2022). In any case , this clearly illustrates the importance of the climate reconstruction in producing these results.\nIn the paragraph between lines 194 to 207 the authors argue that their results have implications on the climate, which I find somewhat circular: the presented results are entirely dependent on both whether there is and overall decrease in the global temperature and whether the assumptions about the number of niches ~ temperature relationship holds. In my opinion, if we did not know anything about the history of the climate, the increase in global diversity and the model's implications would yield only a very weak suggestion about a potentially cooling climate.\nIn conclusion, the results themselves are interesting, and these points do not interfere with their soundness, merely with their interpretation. It is important to add to the discussion of early Phanerozoic life that decreasing temperature is capable of increasing global diversity when these assumptions hold, but the results are not conclusive evidence that cooling is the single or a more likely cause, and that the increase of global diversity by purely biotic processes can be ruled out.\nThe manuscript is written in clear and concise language. Besides these major points, I have only a couple of suggestions for the figures:\nFig. 1. Consider reversing the order of the panels (going from older to younger)\nFig. 2. Some color would be nice for this figure - and a legend would also be useful. Please add to the methods how the standard deviations of simulations was calculated. I assume that these come from the 100 replicates, but I am not sure (also the authors did not indicate whether the sampling from the virtual species pool was using replacement or not).\nFig. 3 The understanding of this figure is a bit difficult because it is not clearly stated between which variables was the correlation assessed. What is the observed/predicted variable? Is it the pseudodiversity values in every raster cell? Or the average for the latitude?\nFig. 5 Could also do better with a bit of color.\nNotes on links and deposited material:\n- I am very glad that the authors decided to deposit the code they used in the analysis.\n- Even though the climate models do not form the backbone of the study, I would recommend the deposition all necessary information to rerun the climate models upon the eventual publication of these results.\n- Consider the long-term deposition of this information on Line 366. My browser indicated that this link might be a security threat and there is no guarantee that this link will stay live for good (it is very annoying to click on link in a published paper just to find out that it is dead.)\nReferences:\nBeaugrand, G., Luczak, C., Goberville, E., & Kirby, R. R. (2018). Marine biodiversity and the chessboard of life. PLOS ONE, 13(3), e0194006. https://doi.org/10.1371/journal.pone.0194006\nGrossman, E. L., & Joachimski, M. M. (2022). Ocean temperatures through the Phanerozoic reassessed. Scientific Reports, 12(1), Article 1. https://doi.org/10.1038/s41598-022-11493-1\nHagen, O., Flück, B., Fopp, F., Cabral, J. S., Hartig, F., Pontarp, M., Rangel, T. F., & Pellissier, L. (2021). gen3sis: A general engine for eco-evolutionary simulations of the processes that shape Earth’s biodiversity. PLoS Biology, 19(7), e3001340. https://doi.org/10.1371/journal.pbio.3001340\nZacaï, A., Monnet, C., Pohl, A., Beaugrand, G., Mullins, G., Kroeck, D. M., & Servais, T. (2021). Truncated bimodal latitudinal diversity gradient in early Paleozoic phytoplankton. Science Advances, 7(15), eabd6709. https://doi.org/10.1126/sciadv.abd6709", "reviewer_3": "The manuscript by Ontiveros et al. presents a modelling study of the relationship between Ordovician global climate cooling and biodiversification.\nAs I cannot comment on the validity of the macroecological modelling, my comments focus mainly on the paleotemperature record.\nFor my evaluation, a more detailed description of the modelling parameters and modelling results is needed (see below). In particular, the modelled latitudinal temperature gradients which are essential for the modelled latitudinal biodiversity gradients should be documented and evaluated. In addition, some misunderstandings need be clarified.\nLine 41: „ to modern-like values (ca. 30°C) in the Early and Middle Ordovician”. This statement is not correct, because oxygen isotope temperatures decrease during the Early Ordovician and most of the Middle Ordovician with temperatures during the Late Ordovician > 30°C, at least in case of the temperature reconstructions published by Grossman & Joachimski 2022 (Scientific Reports 12:8938). This recent compilation of oxygen isotope temperatures shows higher temperatures compared to Song et al. because the latter authors didn’t apply a latitude correction for seawater d18O.\nLine 66-67: pCO2 has been adjusted to model the Ordovician cooling trend. Unfortunately, the assumed change in pCO2 is not documented. It would be interesting to assess how realistic the assumed atmospheric CO2 levels or the Ordovician decrease in pCO2 are.\nLines 80-81: it is argued that the maximum simulated biodiversity shifted from 60° S at 490 to 470 Ma to low latitudes at 450 Ma. What this means to me is that at 490 to 470 Ma, temperatures at 60°S and 40° N must have been considerably colder than in the hot tropics. What about the modelled latitudinal temperature gradient from the equator to 40° N or 60° S. Is this temperature gradient consistent with a relatively flat latitudinal temperature gradient of a greenhouse climate?\nLine 86-88: Lewis & Eichenseer (2021) argue that a variable spatial coverage of data can have an impact on the reconstruction of GLOBAL temperature. However, the temperature records of Trotter et al. (2008) and Song et al. (2019) are low latitude temperature records (Song et al: 0 to 40° latitude). Low latitudes do not show large sea surface temperature differences, so variations in the spatial coverage of the data should not have a significant effect. In addition, the combined Ordovician carbonate (brachiopod) and apatite (conodont) d18O records (Grossman & Joachimski 2022), records restricted to 0 to 30° latitude, show a comparable temperature decrease. Thus, the argument that the low-latitude temperature cooling is underestimated does not seem justified. And citing Lewis and Eichenseer 2021 is not supportive.\nLine 124 to 130: The experiments indicate that thermal limits > 65 °C suppress the impact of global climate cooling on biodiversity. I have problems follow this, as upper thermal limits of 65 to 78 °C are completely unrealistic for metazoans (except certain archaea and bacteria). Upper thermal limits may be at 47° C, possibly reflecting an oxygen limitation mechanism of thermal tolerance and a thermal sensitivity of macromolecular structures (Pörtner 2002). So what do these experiments indicate? Or have I missed something?\nLine 147 to 149: The authors state that in the late Cambrian to Early Ordovician, the maximum modelled biodiversity is at mid to high latitudes. This implies that the tropics were hot and the mid to high latitudes were colder. How steep is the modeled latitudinal temperature gradient? Is this gradient realistic for 0 to 40° N or 0 to 60° S in a greenhouse world with ice-free poles (see comment above)? Unfortunately no information is given. This is crucial because in case of a too steep or unrealistic latitudinal temperature gradient, the whole modeling approach becomes questionable. I recommend that the authors provide this information in the supplementary material.\nLine 197 to 199: “an assumption supported by clumped isotope measurements that show no long-term trend in seawater δ18O over the last 500 Myrs”. This sentence is misleading because clumped isotopes are independent of the oxygen isotope composition of seawater. Clumped isotopes only depend on temperature, e.g. during calcite precipitation.\nFigure 4a and extended data Fig. 3: Why is the Katian missing on the x-axis?"}, "author_response": "Response to the reviewers [NCOMMS-22-52383-T]\nCooling Oceans Did Trigger Ordovician Biodiversification\n[...]"}, {"round": 2, "reviewer_comments": {"reviewer_2": "Upon a second review, I still consider this paper very interesting and worthy of publication.\nOverall, I am happy with the authors' replies to my comments and glad to see them\nimplemented.\nI would still argue that there is a chance that the high number of species in the 15-25°C\ntemperature range (same as the mid domain outlined by the min-maximum range) might\nstill be the result of adaptation to the post-Cambrian climate and could be different in the\npast, potentially changing throughout the Cambrian-Ordovician interval. Again, this can\nhappen without changing the maximum limit and by left-skewing the distribution (the 65°C\nupper limit, when the pattern breaks down, corresponds to a mode of ~33˚C, which is warm,\nbut perhaps not impossible). As the authors state in the rebuttal, this would mean a higher\nconcentration of ranges in warmer temperatures. Competitive exclusion, as the authors\nstate in their rebuttal, is only a problem if species have to occupy the same physical space,\nwhich (as the results of the authors demonstrate this) was larger in high-temperature areas.\nI simply disagree with the authors on the possibility of such a pattern, which is indeed less\nuniformitarian, but given that we are talking about totally different organisms, I think it\nremains an inexcludable eventuality.\nThere is no need to do any more analyses, but for the sake of completeness, I would\nwelcome if the authors mentioned that this effect has the potential to influence the\noutcome of their study (the paragraph from Line 152 seems like a good place to do this).\nI only have a couple of quick comments:\nLine 62: \"between 490 Ma 430 Ma, both included\" - perhaps this can be simplified with\n\"from 490 to 430 Ma\".\nLine 69: \"which captures well biodiversity patterns\" - consider reordering to \"which captures\nbiodiversity patterns well\"\nFigure 4. Note that the figure coloration itself is based on specific levels between contours,\nand the legend is based on continuous colors. The temporal continuity of the plot\nmisleadingly suggests a time-continuous simulation, whereas the original results are 7\nsnapshots. Given that the quality of the patterns on the figure is not superb and the\ninterpolated contours are very noisy, perhaps consider showing rectangular bands instead.\nFigure 5. The lack of a visual legend makes this figure very difficult to read. Please also note\nthat the envelopes that the authors mentioned are not present in the lower-resolution\nfigure (they are present in the high-resolution ones at the end of the merged document).\nAlso consider decreasing the point size.\nI was also asked to provide feedback on the authors' response to Reviewer 1. I am glad that\nReviewer 1 raised the awareness of discussion the latitudinal diversity gradient and I am\ncontent with the replies of the authors.\nIn any case, I would like to congratulate to the authors for the very exciting manuscript,\nwhich will be worthwhile contribution to the literature."}, "author_response": "Response to the reviewers [NCOMMS-22-52383A]\nImpact of global climate cooling on Ordovician marine biodiversity\n[...]"}]

Decision letter