Terrestrial forcing of marine biodiversification

The diversification of the three major marine faunas during the Phanerozoic was intimately coupled to the evolution of the biogeochemical cycles of carbon and nutrients via nutrient runoff from land and the diversification of phosphorus-rich phytoplankton. Nutrient input to the oceans has previously been demonstrated to have occurred in response to orogeny and fueling marine diversification. Although volcanism has typically been associated with extinction, the eruption of continental Large Igneous Provinces (LIPs) is also a very significant, but previously overlooked, source of phosphorus involved in the diversification of the marine biosphere. We demonstrate that phosphorus input to the oceans peaked repeatedly following the eruption and weathering of LIPs, stimulating the diversification of nutrient-rich calcareous and siliceous phytoplankton at the base of marine food webs that in turn helped fuel diversification at higher levels. These developments were likely furthered by the evolution of terrestrial floras. Results for the Meso-Cenozoic hold implications for the Paleozoic Era. Early-to-middle Paleozoic diversity was, in contrast to the Meso-Cenozoic, limited by nutrient-poor phytoplankton resulting from less frequent tectonism and poorly-developed terrestrial floras. Nutrient runoff and primary productivity during the Permo-Carboniferous likely increased, based on widespread orogeny, the spread of deeper-rooting forests, the fossil record of phytoplankton, and biogeochemical indices. Our results suggest that marine biodiversity on geologic time scales is unbounded (unlimited), provided sufficient habitat, nutrients, and nutrient-rich phytoplankton are also available in optimal amounts and on optimal timescales.


Methods
We assessed the roles of phosphorus availability, primary productivity, and nutrient recycling on rates of nannofossil and marine genera origination rates (GOR) for the last 159.5 Ma (Jurassic, ca. late Oxfordian) of the Meso-Cenozoic. Detailed records of phosphorus accumulation rates (PAR), calcareous nannoplankton, isotope data, and genera origination rates are available for this interval. Data was compiled from a range of sources, which used various time scales and bins (see Supplementary Table S1 and Supplementary Fig. S1 for their raw data with their originally assigned ages).
Besides strontium, detailed records of two other stable isotopes are also available for the Meso-Cenozoic: carbon (δ 13 C) and sulfur (δ 34 S). Both carbon and sulfur cycles control redox conditions at the Earth's surface by acting in a reciprocal manner. Positive carbon isotope ratios (δ 13 C) indicate enhanced marine and/or terrestrial photosynthesis, whereas negative ratios indicate decreased photosynthesis and/or input of isotopically-lighter 12 C from various sources such as the erosion and oxidation of terrestrial organic carbon and its input into the oceans. High positive sulfur isotope values (δ 34 S) are interpreted to indicate extensive sulfate (SO 4 2− ) reduction by sulfate-reducing bacteria, which are intolerant of oxygen and use dissolved SO 4 2− as an electron acceptor to oxidize organic matter under anoxic conditions. Conversely, negative δ 34 S values are interpreted to indicate lower rates of sulfate reduction.
Our analyses were conducted using previously-published 11-myr binned data for stable isotopes and GOR; we further establishied 5-myr bins at the approximate mid-points of 11-myr bins in an attempt to achieve finer temporal resolution within the temporal constraints of the raw data.
Data used by Cárdenas and Harries 12 were initially binned by them to 5-myr intervals using Linear Interpolation, including Genera Origination Rate (GOR) data of Alroy 38 reported by him in approximately 11-myr bins. www.nature.com/scientificreports/ All data were then re-binned to 11-myr bins by Cardenas and Harries 12 , mirroring those used in the calculation of origination rates by Alroy 38 , prior to undertaking their statistical analyses. Materials sources are indicated below.
1. Strontium isotope ratios (‰) 12  The total P dataset includes information on phosphate phases that were originally dissolved and bioavailable and those that were detrital, i.e., derived from continental weathering without intermittent biological involvement. Föllmi attempted to discriminate between these two phases by considering a subset of pelagic biogenic sediments 43 . Phosphorus phases in this data subset were considered by him to likely be mainly nondetrital and representing reactive bioavailable phosphorus. The similarity of the total PAR and biogenic PAR curves ( Fig. 2) indicates that changes in total P fluxes were closely tracked by changes in biogenic PAR fluxes or, perhaps more generally according to Föllmi, that changes in total continental weathering rates led to comparable changes in chemical weathering as the main long-term source of dissolved bioavailable phosphorus 43 .
4. Carbon isotope ratios (δ 13 C, ‰) 12 : Data originally from Veizer et al. 45 . Veizer et al. made carbon isotope measurements on calcitic and phosphatic shells, mainly brachiopods with some conodonts and belemnites, collected at high temporal resolution (up to 0.7 Ma or one biozone) from the stratotype sections of all continents but Antarctica and from many sedimentary basins, and mostly from paleotropical domains. Scanning electron microscopy, petrography, cathodoluminescence and trace element analyses of calcitic shells and the conodont alteration index (CAI) of the phosphatic shells indicated excellent preservation of their ultrastructure. These datasets were complemented by extensive literature compilations of Phanerozoic low-Mg calcitic, aragonitic and phosphatic isotope data of analogous skeletons. Data were treated by Cárdenas and Harries as for Sr isotope ratios 12 . 5. Sulfur isotope ratios (δ 34 S, ‰) 12 : Data originally from Kampschulte and Strauss in 5-myr intervals 46 . Sulfur isotopes were measured on biogenic calcite (brachiopods and belemnites) and micritic carbonates, all of Paleozoic and Mesozoic age. Results were supplemented with published sulfur isotope data for Neogene foraminifera, whole rock carbonates across the Cenomanian-Turonian boundary, and Cenozoic marine barites. Data were treated by Cárdenas and Harries as above 12 . 6. Sea-level (meters) 47 Table 1).

Statistical analysis.
We first conducted statistical analyses in Paleontological Analysis Statistical Software (PAST) version 4.03 and reaffirmed the results in R 52,53 . Linear Interpolation was used by us to re-bin all original data to 5-myr intervals approximately midway between 11-myr bins intervals in an attempt to achieve greater temporal resolution while remaining within the temporal constraints of the raw data. In several cases for data not taken from Cárdenas and Harries 12 , we readjusted the youngest ages to exactly 17 Ma so as to obtain exactly uniform 11-myr bins and from which we produced uniform 5-myr bins via linear interpolation. Readjusted ages (see below for original sources) are for: CO 2 (originally 17.5 Ma), Phosphorus Accumulation Rate (PAR, total and biogenic; originally 17.5 Ma) and nannofossil diversification rates (originally 17.788 Ma; Supplementary  Tables S2-S3; Supplementary Figs. S2, S3). We first tested for the normal distribution of the undifferenced 11-myr and 5-myr binned data sets (Normality Tests option under PAST Univariate menu), one or more tests of which (Shapiro-Wilk, Anderson-Darling, χ 2 , and Jacque-Bera) indicated that the undifferenced 11-myr binned GOR were not normally-distributed (p < 0.05 for non-normality) and the p value of the δ 34 S approached that of non-normality (p < 0.05). One or more the same tests indicated that the undifferenced 5-myr binned CO 2 , Sr, δ 34 S, GOR, and sea level data were also nonnormally distributed. Pearson's r therefore proved inapplicable for correlation. www.nature.com/scientificreports/ We therefore correlated using Spearman's correlation coefficient (ρ), which is a non-parametric test that rank orders interval data when it is non-normally distributed. Spearman's ρ has less sensitivity (power, 1-β, the likelihood of accepting the null hypothesis when it should have been rejected) than does Spearman's r; the higher a statistical test's power, the greater the probability that a small difference or correlation will be found to be significant. Unlike Pearson's r, however, Spearman's ρ makes no assumptions regarding the parameter mean or standard deviation of the sampled population. We used a significance level of p < 0.05. We also noted certain correlations which were considered marginally insignificant (0.05 > p < 0.1); these particular correlations may reflect the binning of the original time series, which exhibited different scales of temporal resolution.
Individual time series varied slightly in length; consequently, correlations were conducted only as far back in age as comparisons with the PAR time series permitted: ~ 159 Ma ( Supplementary Fig.S4). Data were differenced to eliminate possible spurious correlations using the Evaluation Expression (u-d) of the Transform Data  Tables S4, S5). We correlated both 11-and 5-myr binned uncorrected data and the same data corrected with the Bonferroni procedure, which has been advocated by some but not all investigators in the natural sciences. Because of these contrasting views, we correlated with and without the Bonferroni correction for the sake of comparison of our results.
The Bonferroni correction has been widely recommended when conducting multiple comparisons to cull null hypotheses which should be rejected. The Bonferroni correction is based on the equation p (corrected) = α/m, where α = desired uncorrected p at the outset (e.g., 0.05 in our study) and m = number of comparisons ("hypotheses" to be tested). The correction has therefore also been severely criticized, as it is extremely conservative. This is because as the number of variables used in multiple comparisons increases, the required p value with the Bonferroni correction for the rejection of the null hypothesis becomes increasingly small to the point that it may begin to produce false negatives (Type I error: null hypothesis rejection when the hypothesis is true) rather than omitting false positives (Type II error: null hypothesis acceptance when the hypothesis is false) [54][55][56] . Some workers therefore advocate using only uncorrected data [54][55][56] . The theoretical basis for advocating an adjustment for multiple comparisons is that "chance" serves as the first-order explanation for observed phenomena. This particular assertion is thought by some workers to undermine the basic premises of empirical research: that nature obeys regular laws that can be studied through observation; omitting the Bonferroni adjustment is therefore thought to be preferable because it will lead to fewer errors of interpretation when the data under evaluation are not randomly distributed but actual observations.

Results
Both 5-and 11-myr binned uncorrected data yielded a significant Spearman's correlation (p < 0.05) between 87 Sr/ 86 Sr and GOR, corroborating the role of nutrient runoff in GOR, as previously reported for the entire Phanerozoic (Table 1) 12 . Strong significant correlations for both bin intervals were also found between biogenic PAR and GOR, directly implicating the bioavailability of phosphorus and its transfer along food chains in GOR.
The increased temporal resolution of five-myr binned uncorrected data further yielded significant positive correlations: biogenic PAR with both total PAR and nannofossil diversification; total PAR with GOR; and nannofossil diversification with GOR; a significant negative correlation between δ 13 C and δ 34 S was also found. These particular correlations for 5-myr binned data were found to be either insignificant (p > 0.1) or, suggestively, marginally insignificant (p > 0.05 but < 0.1) for 11-myr binned data. Other marginally insignificant correlations were also found: negative correlations between δ 34 S and nannofossil diversification for both 5-and 11-myr binned data, and positive correlations for 5-myr binned data between δ 13 C and both nannofossil diversification (p < 0.085) and GOR (p < 0.077).
Sea level exhibited no impact on GOR with either 11 or 5-myr binned uncorrected data, similar to earlier findings for 11-myr time bins (Table 1) 12 . Both binning methods, however, implicated sea level in total PAR and www.nature.com/scientificreports/ carbon burial (δ 13 C) at significant to marginally insignificant levels. Sea level also exhibited a positive correlation with CO 2 in 5-myr binned data, which was only very marginally insignificant (ρ = 0.37, p < 0.508). We conducted a second set of correlations using the Bonferroni correction (see Methods). We chose to eliminate CO 2 and sea level because they did not correlate with GOR using uncorrected data; this allowed us to further test the strength of correlations which had previously yielded significant correlations with uncorrected 11-and 5-myr binned data ( Table 1). Elimination of CO 2 and sea level resulted in a somewhat higher p value of 0.0024 (p = 0.05/21 comparisons). No significant correlations were found for 11-myr binned corrected data, whereas five-myr corrected data again yielded significant correlations between 87 Sr/ 86 Sr, biogenic PAR, total PAR, and nannofossil diversification with GOR (Table 2). These correlations therefore appear robust, given the previous results for 5-myr binned, uncorrected data (cf. Tables 1, 2).
Cross-correlations of lagged indices against GOR weakened Spearman's correlations of both uncorrected and corrected data to insignificance, however, suggesting that the interactions of environmental forcings and the responses of biodiversification to them are occurring primarily on relatively geologically-short time scales at least down to the range of 5-11 myr. The temporal spacing of data points of some time series did not permit us to further refine the time range without the questionable interpolation of further data points between more widely-spaced intervals.

Discussion
Meso-Cenozoic phytoplankton stoichiometry fueled diversification of the modern fauna. The dominant eukaryotic phytoplankton of the Meso-Cenozoic have been allied with so-called "red lineages" (coccolithophorids, dinoflagellates diatoms). These taxa are characterized by the accessory pigment chlorophyll c, specific trace elements employed in their plastids (such as those mentioned above), and low carbon:phosphorus (C:P) ratios 20,37 . Although their stoichiometry varies in response to natural conditions of temperature and nutrient availability, culture studies of modern representatives have been determined to be relatively phosphorusrich and carbon-poor 20,37 . Modeling of C:P ratios occurring at the time of each major taxon's appearance in the geologic record resembles that of their cultured modern representatives, suggesting that the nutrient preferences and stoichiometric compositions of modern representatives are evolutionarily conserved and reflect ancestral conditions rather than modern ones 20,37 . Ecologic stoichiometric theory predicts that increasing the phosphorus content of phytoplankton decreases the amount of energy that consumers must expend to respire excess carbon to obtain inorganic macronutrients like phosphorus; this potentially leaves excess resources like energy to be devoted to metabolic activity, reproduction, and changes in life cycles that could potentially impact diversification 21 .
The appearance of different major taxa of phytoplankton significantly impacted the deposition of deep-sea oozes and the carbon cycle. The advent of deep-sea calcareous oozes with the appearance of calcareous nannoplankton linked tectonism to atmospheric CO 2 concentrations and weathering rates. A higher recycling efficiency between subducted carbon and CO 2 return flux to the atmosphere via volcanism, as occurred during the Meso-Cenozoic in the form of calcareous oozes, would have provided positive feedback on volcanism and CO 2 flux to the atmosphere, enhancing weathering rates, nutrient input to the oceans, and the primary productivity of calcareous phytoplankton 34,57,58 . Based on the fossil record and sedimentary markers, the apparent steep expansion of the Modern Fauna during the Cenozoic era was paralleled by the tremendous expansion of diatoms which are especially phosphorus-rich (via luxury storage), more so than calcareous nannoplankton 14,20,59-62 . The Table 2. Comparison of correlations omitting CO 2 and sea level using Bonferroni correction*. *Bonferroni correction set to p < 0.0024 (see Methods). Upper half of each binning method (above diagonally-arranged blank boxes): p (uncorrelated); lower half: Spearman's correlation coefficient (ρ). Italics: significant correlations (p < 0.0024). www.nature.com/scientificreports/ rain of phosphorus-rich dead organic matter likely promoted phosphorus bio-limitation in the marine realm. These trends are accompanied by a similar rise of both biogenic and total PAR in DSDP and ODP cores; total PAR may also partly reflect authigenic phosphorus precipitation during the last 15 myr of the Neogene (Fig. 2) 63 . The expansion of diatoms during the Cenozoic may have been aided by the evolution and diversification of terrestrial angiosperm floras, as gymnosperm and calcareous phytoplankton diversity began to decline [60][61][62]64,65 . Angiosperm leaf litter, especially that of woody deciduous species, tends to be relatively nutrient and phosphorusrich and decay relatively rapidly as compared to that of gymnosperms [66][67][68] . Diatoms increase their proportions of total phytoplankton biomass and primary production whenever silica is not limiting; grasses, which are silicarich, became widespread during the last half of the Cenozoic in response to decreased atmospheric moisture resulting from late Eocene-Oligocene glaciation, and are considered to have increased silica input to the oceans via their highly soluble phytoliths 61,62 .
Rates of weathering are nevertheless geologically slow, limiting rates of nutrient runoff to the oceans. Nutrient cycling therefore likely remained critical to continued marine diversification and the biogeochemical cycles of phosphorus. This is especially true for the Meso-Cenozoic, as evidenced by increasing rates and depths of bioturbation [11][12][13][14]69 . Rapid decay of dead organic matter and nutrient recycling via sulfate reduction resulting from primary productivity and the secondary productivity of GOR is suggested by the near mirror-image relation between δ 34 S and δ 13 C and their negative correlation ( Fig. 2; Table 1). Rates of weathering appear to have varied, however, on shorter time scales. Biogenic PAR peaked in association with initial LIP emplacement and peak CO 2 concentrations, but they began to decline almost immediately during the succeeding ~30 myr, as weathering inexorably slowed in response to CO 2 drawdown (Fig. 2).
Implications for Paleozoic diversity. We have concentrated on the Meso-Cenozoic because of its more precise chronologies and completeness of the stratigraphic record, especially that of deep-sea cores. Our results nevertheless hold implications for Paleozoic marine biodiversity: we hypothesize that Paleozoic marine biodiversity remained relatively subdued precisely for the same reasons that it rose so dramatically during the Meso-Cenozoic. In terms of ecological theory, marine biodiversity of the Meso-Cenozoic was "unbounded" (unlimited) whereas that of the Paleozoic, although variable, was much more subdued, or "bounded" because of nutrient limitation 70,71 .
In contrast to the relatively continuous and widespread orogeny beginning in the Mesozoic, and especially during the Cenozoic, the early-to-middle Paleozoic portion of the diversity curve was punctuated by widelyseparated peaks associated with peak strontium isotope ratios and orogeny (Fig. 1). Strontium isotope ratios are significantly greater for the Meso-Cenozoic than for the Paleozoic (Mann-Whitney U, p<0.0001, one-tailed), as are primary productivity and sulfate reduction, presumably in response to greater nutrient runoff (both Mann-Whitney U, p<0.0268, one-tailed; data from 12 ). The eruption of continental LIPs is also thought to have been less frequent during the Paleozoic 5,34 . Although the accuracy of the LIP record before ~200 Ma is limited by the much greater duration of time available for erosion, Paleozoic LIPs are still recognizable by such features as dike swarms, and phosphorus input has been implicated in Late Ordovician marine biotic turnover 23,29 . The earlier Paleozoic was also characterized by terrestrial floras consisting of relatively primitive, rootless and shallowrooting taxa. The litter of modern representatives of these taxa decays relatively slowly [64][65][66][67][68] .
These developments were again paralleled by the phytoplankton. Acritarchs were the dominant phytoplankton of the early-to-middle Paleozoic according to the fossil record and represent the organic-walled fraction of the phytoplankton preserved as cysts (resting stages resistant to inimical conditions). Acritarchs have been allied with eukaryotic "green" phytoplankton lineages. Presumed modern representatives of acritarchs are characterized by plastid trace elements which differ from those of red lineages, and by high C:P ratios (low phosphorus content, carbon-rich) that would have limited biodiversification according to ecologic stoichiometric theory 20,21 .
The Permo-Carboniferous appears to have been transitional between the earlier Paleozoic and the Mesozoic in terms of nutrient availability, primary productivity, and the diversification of marine biotas by increasing the ability of ecosystems to sustain higher but still optimal levels of diversity 72 . Acritarchs largely disappeared from the fossil record during the Permo-Carboniferous. Widespread orogenies associated with the formation of Pangea and the spread of deeper-rooting terrestrial floras inland are thought to have encouraged an increase of weathering rates and nutrient runoff, as indicated by strontium and selenium isotopes and phosphorus concentrations in shales, despite increased carbon and nutrient sequestration on land (Fig. 2) 19,[73][74][75] . Acritarchs may actually represent pre-dinoflagellate lineages or perhaps even dinoflagellates themselves based on biomarkers, whereas modern dinoflagellates appear to be intermediate between red and green lineages with regard to C:P ratios and biomarkers 20,76 . As compared to nutrient-enrichment experiments in modern ecosystems and mass and minor extinctions, which may initially lower diversity via eutrophication, geologically-slow nutrient inputs maintain the relative stability of trophic resources in sufficient quantity thought critical to biodiversity and biodiversification 12-14,72,77,78 . Mechanisms of biodiversification. The exact pathways by which trophic resources are transmuted into biological diversity remain poorly-understood 79 . Nutrient availability and productivity nevertheless undoubtedly play significant roles in biodiversification on different time scales, like those associated with stoichometric theory: enhanced metabolism, increased grazing and predation, all of which would have impacted biogeochemical cycles; enhanced resources available for reproduction, population increase and potential dispersal leading to genetic isolation; and life cycle changes 8,9,13,21 . Productivity, which would influence oxygenation, total PAR, carbon burial and sulfate reduction, has been found to increase niche diversity more than area alone. Higher resource density may allow for increased specialization along a resource axis while still maintaining minimum viable population sizes 80  www.nature.com/scientificreports/ Area in the marine realm corresponds to sea level, which appears to be a function in part of tectonism, given its very nearly significant (but still marginally insignificant) correlation with CO 2 . Previous studies have found sea level to either be associated with biodiversity or not 12,81-83 . Comparison of the sea level curve with those of other indices (Fig. 2) indicates that the flooding of shelves undoubtedly plays significant roles in biodiversification akin to the temporal scales of tectonically-driven Sloss seismic sequences (ca. tens-of-millions of years or more) and longer tectonic cycles of ~250-300 myr duration (as indicated by the broad patterns of marine fossil biodiversity documented by previous studies) [81][82][83][84] . Similar to water depth in the open oceans beyond the shelf edge, sea level likely establishes the broad constraints of habitat availability and environmental conditions on primary and secondary productivity (δ 13 C) and microbial sulfate reduction (δ 34 S) via substrate, sediment accumulation rates, carbon and phosphorus burial, water column stratification and its oxygenation, depth and areal extent of the photic zone, and competition for light and nutrients 7,8,85 .

Conclusion
The diversification of the Modern Fauna is coupled to the appearance and diversification of new and more phosphorus-rich phytoplankton taxa. The evolution of the major phytoplankton taxa of the Meso-Cenozoic in turn broadly parallels that of tectonism, the evolution of terrestrial floras, and the evolution of the biogeochemical cycles carbon and nutrients toward the present. Volcanism, long associated with extinction, serves as a rich source of phosphorus, especially with the eruption of mafic Large Igneous Provinces and the spread of angiosperms. The fossil and biogeochemical records for the Meso-Cenozoic suggest that enhanced nutrient availability and nutrient recycling relaxed the constraints of nutrient limitation, allowing marine diversity to remain more-or-less unbounded. In contrast, the Paleozoic was much more nutrient-limited overall, with the Permo-Carboniferous transitional between the early-to-middle Paleozoic and the Meso-Cenozoic based on phytoplankton and biogeochemical records.

Data availability
All data are available in the main text or the supplementary materials online.