Phanerozoic diversity and neutral theory

Abstract

Although Phanerozoic increases in the global richness, local richness, and evenness of marine invertebrates are well documented, a common explanation for these patterns has been difficult to identify. Evidence is presented here from marine invertebrate communities that there is a Phanerozoic increase in the fundamental biodiversity number (θ), which describes diversity and relative abundance distributions in neutral ecological theory. If marine ecosystems behave according to the rules of Hubbell’s Neutral Theory of Biodiversity and Biogeography, the Phanerozoic increase in θ suggests three possible mechanisms for the parallel increases in global richness, local richness, and evenness: (1) an increase in the per-individual probability of speciation, (2) an increase in the area occupied by marine metacommunities, and (3) an increase in the density (per-area abundance) of marine organisms. Because speciation rates have declined over time and because there is no clear evidence for an increase in metacommunity area through the Phanerozoic, the most likely of these is an increase in the spatial density of marine invertebrates over the Phanerozoic, an interpretation supported by previous studies of fossil abundance. This, coupled with a Phanerozoic rise in body size, suggests that an increase in primary productivity through time is the primary cause of Phanerozoic increases in θ, global richness, local richness, local evenness, abundance, and body size.

Introduction

The species richness and evenness of local marine invertebrate communities have increased through the Phanerozoic (Bambach 1977; Powell and Kowalewski 2002; Bush and Bambach 2004; Kowalewski et al. 2006; Alroy et al. 2008). 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. 1981; Alroy et al. 2008). The connection between these parallel trends has been unclear, in part because local richness is typically much less than 1% of global richness (Holland 2010). 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. 2001; Powell and Kowalewski 2002). 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 2001) 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 2001). 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 2005; Ricklefs 2006; Purves and Turnbull 2010) 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. 2011), including a Phanerozoic decline in speciation rates (Wang et al. 2013). Therefore, neutral theory serves as a useful baseline for understanding biodiversity and biogeography (Rosindell et al. 2012).

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 2001), 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 (2007) 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 (2009) 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. 1990), the end-Jurassic (Hallam 1986), and the Cenomanian/Turonian (Elder 1987), 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. 2000). 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 2001). In exceptional cases, θ can exceed 40 and approach 80, again, within the range of θ in modern settings (Hubbell 2001). Some of the variation in θ reflects metacommunity size, as predicted by neutral theory (Hubbell 2001), with spatially larger metacommunities having larger values of θ, as has been shown in the Ordovician of Laurentia (Sclafani and Holland 2013). 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 ²) 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 2012), 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 2009), 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 2000; Bush and Bambach 2004), 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. (2006) 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 2001). 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 2007; Etienne 2009). 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 2004; Volkov et al. 2007; Tomašových and Kidwell 2010), 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. 1981; Alroy et al. 2008) and local diversity (Bambach 1977; Powell and Kowalewski 2002; Bush and Bambach 2004; Kowalewski et al. 2006; Alroy et al. 2008). 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 2002). 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. 2006).

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 1998; Alroy 2008) 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 2011), 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 1994; Martin 2003; Finnegan and Droser 2008; Smith and McGowan 2008; Pruss et al. 2010; Li and Droser 1999).

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 1975; Martin 1996; Allmon and Martin 2014), indicated by Phanerozoic increases in the average body size of marine organisms, total marine biomass, and metabolic rates (Vermeij 1987; Bambach 1993; Finnegan et al. 2011; Payne et al. 2014). Similarly, increased productivity has been linked to higher biodiversity in modern environments (Chase 2010). 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 2008).

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. 2008; Sepkoski et al. 1981), local richness (Bambach 1977; Powell and Kowalewski 2002; Bush and Bambach 2004; Kowalewski et al. 2006), local evenness (Powell and Kowalewski 2002), abundance (Kidwell and Brenchley 1994; Li and Droser 1999; Martin 2003; Finnegan and Droser 2008; Smith and McGowan 2008; Pruss et al. 2010), and body size (Bambach 1993; Finnegan et al. 2011; Payne et al. 2014). 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.