Introduction
The species richness and evenness of local marine invertebrate communities have increased through the Phanerozoic (Bambach Reference Bambach1977; Powell and Kowalewski Reference Powell and Kowalewski2002; Bush and Bambach Reference Bush and Bambach2004; Kowalewski et al. Reference Kowalewski, Kiessling, Aberhan, Fürsich, Scarponi, Barbour Wood and Hoffmeister2006; Alroy et al. Reference Alroy, Aberhan, Bottjer, Foote, Fürsich, Harries, Hendy, Holland, Ivany, Kiessling, Kosnik, Marshall, McGowan, Miller, Olszewski, Patzkowsky, Peters, Villier, Wagner, Bonuso, Borkow, Brenneis, Clapham, Fall, Ferguson, Hanson, Krug, Layou, Leckey, Nürnberg, Powers, Sessa, Simpson, Tomašových and Visaggi2008). These trends mirror those of global marine biodiversity, a pattern that is robust to taxonomic level, methods of diversity tabulation, and sample standardization (Sepkoski et al. Reference Sepkoski, Bambach, Raup and Valentine1981; Alroy et al. Reference Alroy, Aberhan, Bottjer, Foote, Fürsich, Harries, Hendy, Holland, Ivany, Kiessling, Kosnik, Marshall, McGowan, Miller, Olszewski, Patzkowsky, Peters, Villier, Wagner, Bonuso, Borkow, Brenneis, Clapham, Fall, Ferguson, Hanson, Krug, Layou, Leckey, Nürnberg, Powers, Sessa, Simpson, Tomašových and Visaggi2008). The connection between these parallel trends has been unclear, in part because local richness is typically much less than 1% of global richness (Holland Reference Holland2010). Furthermore, the increase in evenness through the Phanerozoic has been viewed primarily as a problem for sample standardization of diversity, because the relative order of standardized diversity values can depend on the sample-size quota (Alroy et al. Reference Alroy, Marshall, Bambach, Bezusko, Foote, Fürsich, Hansen, Holland, Ivany, Jablonski, Jacobs, Jones, Kosnik, Lidgard, Low, Miller, Novack-Gottshall, Olszewski, Patzkowsky, Raup, Roy, Sepkoski, Sommers, Wagner and Webber2001; Powell and Kowalewski Reference Powell and Kowalewski2002). If estimated diversity depends on sample-size quota, then changes in evenness are a concern because they prevent a unique answer to how diversity has changed through the Phanerozoic.
The Unified Neutral Theory of Biodiversity and Biogeography (Hubbell Reference Hubbell2001) addresses diversity at scales from a local community to a province or metacommunity, and therefore has the potential to offer a unified explanation for changes in diversity at local, regional, and global scales. Neutral theory simulates diversity and relative abundance structure through models of birth, death, dispersal, and speciation within a single trophic level. In neutral theory, diversity and relative abundance are described over a wide range of spatial scales by θ, a recurring parameter in the models that is defined as two times the metacommunity size (measured in numbers of individuals) multiplied by the per-individual probability of speciation (Hubbell Reference Hubbell2001). Metacommunity size is the product of metacommunity area and the density (per-area abundance) of organisms in the metacommunity. Metacommunities and communities with larger values of θ have flatter relative abundance distributions, greater evenness, and greater richness than those with smaller values of θ.
Neutral theory has two critical assumptions. First, a metacommunity is assumed to have a fixed number of sites that can be occupied by organisms, and those sites are always occupied; this is known as the zero-sum rule. Second, neutral theory assumes that all individuals of all species are competitively equal, such that long-term changes in the abundance of any given species are controlled by ecological drift, not by niche characteristics. These assumptions have attracted much criticism (Chase Reference Chase2005; Ricklefs Reference Ricklefs2006; Purves and Turnbull Reference Purves and Turnbull2010) and are unlikely to be strictly true. Even so, neutral theory successfully predicts many aspects of biodiversity and biogeography even with modest departures from these assumptions (Rosindell et al. Reference Rosindell, Hubbell and Etienne2011), including a Phanerozoic decline in speciation rates (Wang et al. Reference Wang, Chen, Fang and Pacala2013). Therefore, neutral theory serves as a useful baseline for understanding biodiversity and biogeography (Rosindell et al. Reference Rosindell, Hubbell, He, Harmon and Etienne2012).
Materials and Methods
Relative abundance data from shallow marine fossil communities were obtained from the Paleobiology Database (paleobiodb.org, download June 2014). Supplemental data on depositional environment, lithology, lithification, and geologic age were also downloaded from the Paleobiology Database. Collections containing only a single species, fewer than five individuals, or no numerical abundance values were removed.
Collections were grouped into data sets, with each data set representing a single reference source and containing one or more collections from the same geographic region, geologic age, and depositional environment. Each data set therefore contains replicate collections from the same setting and is regarded as a sample of a metacommunity. Most data sets have five or fewer collections, but some have as many as 213. Overall, 1140 data sets with a total of 7916 collections were analyzed. Analyzed collections are included in Supplementary Appendix 1.
Because neutral theory is based on diversity dynamics at a single trophic level (Hubbell Reference Hubbell2001), this study focuses on first-order consumers, specifically suspension feeders and deposit feeders. For each data set, only species belonging to classes or orders that consist primarily of suspension and deposit feeders were included (Table 1). In most collections, this culling results in the removal of a few producers (algae) and predators (nautiloids and vertebrates, for example). Trophic information was determined from the Paleobiology Database.
Table 1 Classes and orders consisting primarily of suspension and deposit feeders that were included in this analysis.
An abundance matrix, with collections in rows and taxa in columns, was prepared for each data set, and θ was calculated from the species abundance distributions of each data set. For each data set, the best-fit θ was estimated using the Etienne (Reference Etienne2007) likelihood method, which produces a single estimate of θ when using all collections from a data set. This method also produces a single estimate of the migration parameter m, which was not used in this analysis. Tests using the Etienne (Reference Etienne2009) likelihood method, which allows for a different value of m for each collection, but a single overall value of θ, showed that the values of θ did not differ between the 2007 and 2009 method, and that the 2007 method was substantially faster. Etienne’s methods, which are available as an online supplement to his articles, run in the PARI/GP algebra system, available as a free download and run within a UNIX terminal.
Estimates of θ for all data sets are included in Supplementary Appendix 2, as are data on sample size, depositional environments, rock type, and lithification.
Results
The value of θ in marine invertebrate suspension-feeding and deposit-feeding metacommunities increases through the Phanerozoic, in both its median value and its variance (Fig. 1). Although data coverage is sparse during some periods, the intervals of relatively dense sampling indicate a first-order trend of an increase in θ through the Phanerozoic. Through the Silurian, median θ is generally less than 5, and slowly increases to values generally above 5 by the Recent, with a slope of 0.008 (95% bootstrapped confidence is 0.004 – 0.012). Variance likewise increases erratically through the Phanerozoic (Fig. 2).
Figure 1 Changes in the median fundamental biodiversity number (θ) through the Phanerozoic, plotted by the Paleobiology Database 10-Myr bins and based on 1140 data sets containing 7916 total collections. The bootstrap-based 95% confidence interval is shown in gray. Black arrows indicate the timing of the five major mass extinctions, with gray arrows indicating three other known extinctions.
Figure 2 The fundamental biodiversity number (θ) for each of the 1140 data sets containing 7916 collections in aggregate. Because the distributions are right-skewed, the base-10 logarithm of θ is plotted to illustrate the distributions better. Darker grays indicate overlapping data points.
The time series is marked by several abrupt drops in θ. Five of these correspond to well-known global mass extinctions in the Late Ordovician, Late Devonian, end-Permian, end-Triassic, and end-Cretaceous (Fig. 1). Three other pronounced drops in θ also correspond to extinction events in the early Carboniferous (Raymond et al. Reference Raymond, Kelley and Lutken1990), the end-Jurassic (Hallam Reference Hallam1986), and the Cenomanian/Turonian (Elder Reference Elder1987), although the last of these has also been interpreted as only an apparent decline in diversity caused by changes in the preserved stratigraphic record (Gale et al. Reference Gale, Smith, Monks, Young, Howard, Wray and Huggett2000). Values of θ typically continue to decline following extinction events, and pre-extinction intervals are commonly local maxima. Whether these latter two patterns are robust should be investigated in higher-resolution regional studies of these events.
Within any individual time interval, θ varies markedly and is right-skewed (Fig. 2). Estimates of θ therefore tend to be lower in intervals where data are sparse. θ is generally less than 40, as is common in many modern examples (Hubbell Reference Hubbell2001). In exceptional cases, θ can exceed 40 and approach 80, again, within the range of θ in modern settings (Hubbell Reference Hubbell2001). Some of the variation in θ reflects metacommunity size, as predicted by neutral theory (Hubbell Reference Hubbell2001), with spatially larger metacommunities having larger values of θ, as has been shown in the Ordovician of Laurentia (Sclafani and Holland Reference Sclafani and Holland2013). Variations in speciation rate might also contribute to the variation in θ, as may differences in the spatial density of organisms within those metacommunities, although the contributions of these two factors to regional variation in θ cannot be evaluated at present.
Several biases that might produce an apparent increase in θ over the Phanerozoic were tested. First, because the data sets vary in the number of collections, we tested whether the estimates of θ varied with the number of collections on which they are based (Fig. 3). Although the number of collections in each data set varies from 1 to 213, most data sets have five or fewer collections. Throughout the Phanerozoic, values of θ based on many collections are fully interspersed with values of θ based on few collections, and the value of θ does not increase with the number of collections. Overall, the coefficient of determination (R 2) between θ and the number of collections is 0.031, indicating that the number of collections has no substantial effect on the estimate of θ.
Figure 3 Estimates of log θ for all data sets in the study, coded by the number of collections in each data set.
Second, because taxonomic composition varies markedly across depositional environments (Patzkowsky and Holland Reference Patzkowsky and Holland2012), there is a possibility that θ might vary systematically among depositional environments and that systematic changes in these depositional environments through time might produce an apparent increase in θ. Coding data sets by the depositional environment suggests no systematic pattern through time (Fig. 4). Although there is a statistically significant difference in θ among environments (Kruskal-Wallis chi-squared=22.6, df=6, p-value=0.0009), it is driven by the low values of θ in the comparatively rare estuary and foreshore environments (61 of 1140 collections). The large number of collections (1140) also contributes to the low p-value. The two environments that show the highest values of θ (shallow subtidal and marine indeterminate) are visually no more common in the post-Jurassic than in Jurassic and earlier strata, suggesting that changes in depositional environment are not responsible for the increase in θ. Furthermore, the percentage of data sets from those two environments drops from 69% in the Permian–Jurassic to 55% in the Cretaceous–Cenozoic, the opposite of what would drive an increase in θ.
Figure 4 Estimates of log θ for all data sets in the study, coded by the depositional environment recorded in the Paleobiology Database.
Third, differences in fossil preservation between carbonate and siliciclastic lithologies might cause differences in θ. If this were true, and if the ratio of these lithologies changed systematically through time, it could produce an apparent temporal change in θ. Coding data sets by rock type indicates no visual relationship between lithology and θ (Fig. 5). Similarly, the two dominant lithologic types (carbonate, siliciclastic) do not have statistically distinguishable median θ (randomization p-value=0.13). Likewise, mixed lithologies in the post-Jurassic display a similar range to carbonate values in the pre-Jurassic, suggesting that changes in lithology do not generate the secular trend in θ.
Figure 5 Estimates of log θ for all data sets in the study, coded by primary lithology.
Fourth, because fossils are easier to extract from unlithified samples than lithified samples, and because unlithified samples are more common in the Cenozoic than the pre-Cenozoic (Hendy Reference Hendy2009), the rise in θ has the potential to be driven by changes in lithification over time. A secular trend in lithification is apparent, with poorly lithified and unlithified samples more common in the post-Jurassic (Fig. 6). Median θ is marginally different in the lithified versus the combined poorly lithified and unlithified samples (randomization p-value=0.013), but the difference in median θ is only 1.7, insufficient to drive all of the Phanerozoic trend.
Figure 6 Estimates of log θ for all data sets in the study, coded by degree of lithification.
Finally, enhanced preservation of aragonite can cause differences in diversity (Cherns and Wright Reference Cherns and Wright2000; Bush and Bambach Reference Bush and Bambach2004), and this might cause greater values of θ in samples with aragonite preservation than in those that lack it. Comprehensive data on the preservation of aragonite in these samples are lacking, as this was often not recorded in the original studies from which the data were entered into the Paleobiology Database. In many of the early Paleozoic collections with which the authors are familiar, aragonitic mollusks are preserved as molds, suggesting that changes in aragonitic preservation are unlikely to drive the Phanerozoic increase in θ.
Overall, the lack of correlation of θ with these confounding factors suggests that the Phanerozoic increase in θ is a true biological signal and not merely the result of fossil preservation or sampling. Using data that partly overlap those of this study, Wagner et al. (Reference Wagner, Kosnik and Lidgard2006) showed that relative abundance distributions shift through the Phanerozoic from relatively simple geometric to relatively complex lognormal distributions, and that this change is biologically real. Such a change is another manifestation of an increase in θ through the Phanerozoic (Hubbell Reference Hubbell2001). That these studies reached similar conclusions with overlapping data sets yet different methods suggests that the Phanerozoic trend in θ reflects changes in the actual diversity and relative abundance structure of marine invertebrates.
Discussion
The value of θ can be understood in two ways, one descriptive and one interpretive. A best-fit θ can be determined for any relative abundance distribution, and θ is usually estimated with likelihood methods (Etienne Reference Etienne2007; Etienne Reference Etienne2009). Like any metric used to describe the shape of a relative abundance distribution, θ could be thought of as simply a shape parameter not necessarily having any particular interpretive use or suggesting anything about the causes of a given diversity or relative abundance distribution. Thus, the secular rise in θ could be regarded purely as a description of changes in diversity and relative abundance through the Phanerozoic, much as evenness, local richness, and global richness are.
However, if communities and metacommunities behave as modeled in neutral theory, then θ gains interpretative value because it reflects the result of birth, death, immigration, and speciation in a saturated landscape where individuals of all species are competitively equal. If marine invertebrate communities behave according to these rules (Olszewski and Erwin Reference Olszewski and Erwin2004; Volkov et al. Reference Volkov, Banavar, Hubbell and Maritan2007; Tomašových and Kidwell Reference Tomašových and Kidwell2010), even with modest departures from them, the Phanerozoic increase in θ would have three important implications and provide a causal mechanism linking a broad suite of previously recognized Phanerozoic trends.
First, because θ reflects the product of speciation rate and metacommunity size, it describes their control on the diversity of metacommunities. Similarly, diversity in a local community is described by θ and a migration parameter, m, which is the probability that a death in a local community will be replaced from the metacommunity. Thus, θ links changes in diversity at the local scale to those at the metacommunity scale, and observed changes in θ imply a unified cause for changes in richness at all spatial scales. Changes in θ would therefore provide the critical link to unify explanations for the well-documented Phanerozoic increases in global diversity (Sepkoski et al. Reference Sepkoski, Bambach, Raup and Valentine1981; Alroy et al. Reference Alroy, Aberhan, Bottjer, Foote, Fürsich, Harries, Hendy, Holland, Ivany, Kiessling, Kosnik, Marshall, McGowan, Miller, Olszewski, Patzkowsky, Peters, Villier, Wagner, Bonuso, Borkow, Brenneis, Clapham, Fall, Ferguson, Hanson, Krug, Layou, Leckey, Nürnberg, Powers, Sessa, Simpson, Tomašových and Visaggi2008) and local diversity (Bambach Reference Bambach1977; Powell and Kowalewski Reference Powell and Kowalewski2002; Bush and Bambach Reference Bush and Bambach2004; Kowalewski et al. Reference Kowalewski, Kiessling, Aberhan, Fürsich, Scarponi, Barbour Wood and Hoffmeister2006; Alroy et al. Reference Alroy, Aberhan, Bottjer, Foote, Fürsich, Harries, Hendy, Holland, Ivany, Kiessling, Kosnik, Marshall, McGowan, Miller, Olszewski, Patzkowsky, Peters, Villier, Wagner, Bonuso, Borkow, Brenneis, Clapham, Fall, Ferguson, Hanson, Krug, Layou, Leckey, Nürnberg, Powers, Sessa, Simpson, Tomašových and Visaggi2008). Changes in diversity at both scales would thus have a common underlying origin in the factors that control θ, specifically speciation and metacommunity size.
Second, the factors that control θ drive not only diversity in neutral theory but also the shape of relative abundance distributions. Larger values of θ result in flatter relative abundance distributions characterized by greater evenness than those produced by smaller values of θ. Thus, the Phanerozoic rise in θ would also explain the previously documented parallel increase in evenness (Powell and Kowalewski Reference Powell and Kowalewski2002). Rather than the increase in evenness being a factor that complicates the interpretation of standardized diversity, the increase in evenness is another manifestation of an increasing θ and its effects on diversity and diversity structure. The Phanerozoic rise in θ would also explain the Phanerozoic shift from simple geometric to complex lognormal distributions (Wagner et al. Reference Wagner, Kosnik and Lidgard2006).
Third, the mathematical definition of θ identifies three possible causes for its increase over the Phanerozoic: (1) an increase in the per-individual probability of speciation, (2) an increase in the area of shallow-marine settings, and (3) an increase in the spatial density of organisms in shallow-marine ecosystems. It is difficult to compare per-species speciation rates measured from the fossil record with the per-individual speciation rate of neutral theory, but the Phanerozoic decline in speciation rates (Sepkoski Reference Sepkoski1998; Alroy Reference Alroy2008) makes it unlikely that the Phanerozoic rise in θ was caused by an increase in the per-individual speciation rate. Similarly, an increase in shallow-marine area is not a likely cause of the increase in θ. The area of shallow-marine ecosystems has waxed and waned over the Phanerozoic, it has not shown any long-term trend, and shallow-marine area is lower in the Neogene than in the Ordovician (Hannisdal and Peters Reference Hannisdal and Peters2011), opposite to the trend needed to generate increasing θ. Although absolute abundance of organisms is difficult to infer from the fossil record, a growing body of evidence suggests that the abundance of marine organisms has increased over the Phanerozoic (Kidwell and Brenchley Reference Kidwell and Brenchley1994; Martin Reference Martin2003; Finnegan and Droser Reference Finnegan and Droser2008; Smith and McGowan Reference Smith and McGowan2008; Pruss et al. Reference Pruss, Finnegan, Fischer and Knoll2010; Li and Droser Reference Li and Droser1999).
Of the three mechanisms possibly underlying the secular increase in θ, an increase in the spatial density of organisms is the best supported and therefore the most likely. Such an increase in the density and abundance of marine life is consistent with previous interpretations of increasing primary productivity through the Phanerozoic (Jackson Reference Jackson1975; Martin Reference Martin1996; Allmon and Martin Reference Allmon and Martin2014), indicated by Phanerozoic increases in the average body size of marine organisms, total marine biomass, and metabolic rates (Vermeij Reference Vermeij1987; Bambach Reference Bambach1993; Finnegan et al. Reference Finnegan, McClain, Kosnik and Payne2011; Payne et al. Reference Payne, Heim, Knope and McClain2014). Similarly, increased productivity has been linked to higher biodiversity in modern environments (Chase Reference Chase2010). If the Phanerozoic changes in diversity are driven primarily by changes in the spatial density or abundance of organisms, in turn caused by changes in primary productivity, it is noteworthy that all five major mass extinctions, plus three additional extinction events, are marked by abrupt drops in θ. This pattern suggests that mass extinctions were associated not only with a loss of diversity, but also with a loss of abundance, possibly triggered by a drop in primary productivity. This interpretation is complicated, however, by the common association of anoxia with mass extinction, which could reflect an increase in primary productivity (Meyer and Kump Reference Meyer and Kump2008).
Three decades of paleobiological research have documented a rich array of trends in marine systems through the Phanerozoic, including increases in global richness (Alroy et al. Reference Alroy, Aberhan, Bottjer, Foote, Fürsich, Harries, Hendy, Holland, Ivany, Kiessling, Kosnik, Marshall, McGowan, Miller, Olszewski, Patzkowsky, Peters, Villier, Wagner, Bonuso, Borkow, Brenneis, Clapham, Fall, Ferguson, Hanson, Krug, Layou, Leckey, Nürnberg, Powers, Sessa, Simpson, Tomašových and Visaggi2008; Sepkoski et al. Reference Sepkoski, Bambach, Raup and Valentine1981), local richness (Bambach Reference Bambach1977; Powell and Kowalewski Reference Powell and Kowalewski2002; Bush and Bambach Reference Bush and Bambach2004; Kowalewski et al. Reference Kowalewski, Kiessling, Aberhan, Fürsich, Scarponi, Barbour Wood and Hoffmeister2006), local evenness (Powell and Kowalewski Reference Powell and Kowalewski2002), abundance (Kidwell and Brenchley Reference Kidwell and Brenchley1994; Li and Droser Reference Li and Droser1999; Martin Reference Martin2003; Finnegan and Droser Reference Finnegan and Droser2008; Smith and McGowan Reference Smith and McGowan2008; Pruss et al. Reference Pruss, Finnegan, Fischer and Knoll2010), and body size (Bambach Reference Bambach1993; Finnegan et al. Reference Finnegan, McClain, Kosnik and Payne2011; Payne et al. Reference Payne, Heim, Knope and McClain2014). Understanding and demonstrating the causal connections among these patterns has been elusive, but the recognition of a Phanerozoic increase in θ provides that causal link. Taken as a whole, these patterns point to a dominant role for productivity in driving Phanerozoic changes in marine ecosystems.
Acknowledgments
We thank A. Platsky for assistance in data collection, R. Etienne for assistance with calculating θ, and M. Foote for comments on an early draft of the manuscript. We also thank J. Payne and an anonymous reviewer for their helpful suggestions. This research was supported by National Science Foundation grant EAR-0948895 (S.M.H.). This is Paleobiology Database publication no. 219.
Supplementary Material
Supplemental materials deposited at Dryad: doi:10.5061/dryad.n1d61.
