Climate change during Late Ordovician mass extinction drove an ecologically selective decline in the abundance of deep-water macrozooplankton, accompanied by a shift to simpler, less even communities. These results indicate that the species abundance structure of planktic communities may be a leading indicator of the effects of climate change on biodiversity and a more sensitive gauge of those effects than taxonomic diversity change alone. Additionally, we present a Bayesian likelihood method of habitat inference that may be applicable to other similar cases.

Mass extinctions disrupt ecological communities. Although climate changes produce stress in ecological communities, few paleobiological studies have systematically addressed the impact of global climate changes on the fine details of community structure with a view to understanding how changes in community structure presage, or even cause, biodiversity decline during mass extinctions. Based on a novel Bayesian approach to biotope assessment, we present a study of changes in species abundance distribution patterns of macroplanktonic graptolite faunas (∼447–444 Ma) leading into the Late Ordovician mass extinction. Communities at two contrasting sites exhibit significant decreases in complexity and evenness as a consequence of the preferential decline in abundance of dysaerobic zone specialist species. The observed changes in community complexity and evenness commenced well before the dramatic population depletions that mark the tipping point of the extinction event. Initially, community changes tracked changes in the oceanic water masses, but these relations broke down during the onset of mass extinction. Environmental isotope and biomarker data suggest that sea surface temperature and nutrient cycling in the paleotropical oceans changed sharply during the latest Katian time, with consequent changes in the extent of the oxygen minimum zone and phytoplankton community composition. Although many impacted species persisted in ephemeral populations, increased extinction risk selectively depleted the diversity of paleotropical graptolite species during the latest Katian and early Hirnantian. The effects of long-term climate change on habitats can thus degrade populations in ways that cascade through communities, with effects that culminate in mass extinction.

Mass extinctions affect communities over a range of hierarchical levels of organization (1, 2). Droser et al. (1) described four levels of community impacts, ranging from the broad scale of the development or loss of entire ecosystems to the fine scale of abundance changes within a single community. Studies of community changes among fossil marine invertebrates based on the occurrence of genera or broad functional groups have yielded interesting insights about the nature of mass extinctions, particularly the apparent decoupling of taxonomic and ecological impacts of mass extinction (e.g., refs. 1 and 2). Discussions of the ecological impact of mass extinction typically focus on the idea that communities are disrupted by the loss of species diversity (3, 4), or keystone species (5). Few studies (e.g., refs. 6 and 7) have examined within-community changes using abundance data collected from fossil invertebrate communities, but those few have produced novel and unexpected insights. In their analyses of terrestrial plant communities, where species (or generic) abundance data are the norm, McElwain et al. (8) demonstrated that fine-scale community changes may occur before mass extinctions as well as in the wake of those events, and that relative abundance curves revealed changes in community structure over the time span of a mass extinction event. Similarly, recent work on the response of modern biotas to anthropogenic effects on habitats and species distributions indicates that striking changes in the species abundance structure of communities may precede, and possibly presage, species losses (9, 10).

We consider here the fine-scale impact of climate change on planktonic communities during the Late Ordovician mass extinction (LOME) and the immediately preceding interval (∼447–444 Ma). The loss of graptolite biodiversity in the LOME, accompanied by the wholesale extinction of the previously dominant Diplograptina (taxonomic use follows ref. 11) and their replacement by the previously marginal, high-latitude Neograptina (12⇓⇓⇓–16), provides an opportunity to study the impact of climate change on a macroplanktonic invertebrate fauna over several million years during an interval of unusual species turnover (17, 18). A focus on climate change over geological timescales as a driver of extinction dynamics leads us to ask whether there is evidence of ecological community changes in the interval leading up to mass extinction. Do ecological communities collapse as species go extinct, or do environmental conditions cause disruptions of community structures that result in species extinctions, or do both species and communities react in similar ways to environmental insults? Examination of marine planktonic community structures before the main pulse of a mass extinction may yield insights about both the dynamics of the extinction process and the impact of community structure on changes on biodiversity—insights that may be of use in our efforts to understand the current accelerated rates of community change and species loss.

Planktic graptolite species may be grouped into two major ecological groups based on their habitat affinities, corresponding patterns of resource utilization, and species turnover dynamics (19, 20). One group is thought to have inhabited a deep-water biotope within the twilight or mesopelagic zone below the surface mixed layer (Fig. S1), with their peak diversity and abundance associated with productivity driven by nutrient upwelling at the edges of continental shelves from the oxygen minimum zone of a highly stratified ocean. The majority of graptolite species, however, appear to have occupied the relatively sunlit, surface mixed layer (epipelagic zone) of the ocean. The mesopelagic zone today is commonly separated from the overlying surface mixed zone by the pycnocline. This hydrographical boundary may have contributed to the segregation of pelagic graptolites within two more-or-less discrete habitats: the mesopelagic and epipelagic biotopes. After death, all these graptolites settled to form mixed assemblages on the sea floor, with epipelagic graptolites present at relatively onshore as well as at more offshore localities, whereas mesopelagic species occur only in strata deposited near the shelf margin and beyond. Changes in deep-ocean circulation and oxygenation driven by climate change contributed to the Ordovician mass extinction (12, 13, 15, 21, 22), which suggests that the mesopelagic biotope may have been the most vulnerable and should exhibit disproportional extinction.

Conceptual ecological model of the distribution of planktic graptolite communities in the water column along the shelf margins of the Late Ordovician continents and their relationship to physical and biological properties of the oceans. The epipelagic species appear to have dwelt mainly in the surface mixed zone above the pycnocline, whereas another group of species was associated with oxygen minimum zones in the twilight, mesopelagic zone of the water column. Also shown is the resulting distribution of death assemblages on the sea floor beneath these zones. Modified with permission from ref. 20 .

The LOME is associated with a shift from the greenhouse climate that had dominated most of the Cambrian to Mid-Ordovician metazoan macroevolutionary radiations to icehouse conditions through an interval of increasingly cool and unstable conditions during the Late Ordovician (23⇓–25). Hirnantian (latest Ordovician) glaciation lowered sea levels by some 70–100 m (26, 27), reduced tropical sea surface temperatures by ∼6 °C (22, 26), and increased oxygenation of the deep oceans (22, 28, 29). In combination, these changes appear to have reduced the supply of nutrient-rich, denitrified waters, which led to changes in the composition of the phytoplankton community at most sites formerly occupied by the diverse shelf-edge graptolite communities (21, 22, 30).

In this study, we ask whether the ecological composition and structure of graptolite communities exhibited changes during the interval of climate deterioration and gradual species losses that preceded the onset of major Hirnantian glaciation, and presaged the dramatic taxonomic turnover that accompanied that event. Further, we ask whether those changes in community abundance structure provide novel insights into the mass extinction that differ from those suggested by the well-documented history of changes in the species composition of Late Ordovician faunas (17).

We address these questions by analysis of species abundance patterns through measured stratigraphic sections at two localities from the northwestern passive margin of Laurentia (Fig. S2 and Datasets S1 and S2). The Vinini Creek site (Nevada; hereafter VC) records deposition in a bathyal continental slope to ocean floor setting, slightly south of the paleoequator (15, 31). The Blackstone River site (Ogilvie Mountains, Yukon, Canada; hereafter BR) was a shallower site located ∼3,000 km away, on the continental shelf, north of the paleoequator (32, 33). Additionally, paleoenvironmental proxy data are available at these sites (21, 22, 30, 34) that permit independent characterization of paleooceanographic changes. These data include ε Nd values (Fig. 1), from which we reconstruct relative sea level changes at BR and VC (see ref. 34 and SI Text for further discussion) that we denoted as cycles O4a to O6b, following refs. 35 and 36.

Species biotope assignments, abundances and community structure obtained from bulk sampled Late Ordovician graptolite faunas plotted against schematic representations of the stratigraphic successions ( 21 , 34 , 41 ) and geochemical proxies for environmental change at (A) BR (note expanded vertical scale in shaded interval at this site) and (B) VC. Carbon isotope values (δ 13 C organic and δ 13 C carbonate ) ( 21 ) and hopane/sterane biomarker ratios ( 30 ) document the timing of faunal change relative to changes in phytoplankton communities and the basal Hirnantian Elkhorn and the mid-Hirnantian (HICE) carbon isotopic excursions ( 22 ). T/R cycles are relative sea level change: transgressions and regressions, respectively, separated by inferred sequence boundaries (SB) based on ε Nd values and interpretation in ref. 34 . Values for ε Nd at graptolite sample levels are based on piece-wise linear interpolation between flanking measured values. Symbols for Shannon information, evenness, percent of mesopelagic species, and interpolated ε Nd values correspond to the SAD type for each sample (key at top of figure). Species abundance data and key to numbering are presented in Datasets S4 and S5 . Species are grouped according to inferred biofacies affiliation and clade membership; sizes of symbols for species occurrences are proportional to the log of the number of attributed specimens in each collection. See ref. 34 for explanation of lithological symbols.

The onset of the glaciation early in M. extraordinarius Zone time is marked at both VC and BR by a decrease in ε Nd , suggesting relative sea level fall, a strong positive excursion in δ 13 C, and the abrupt replacement of the diplograptine fauna by neograptine species. At both sites, the replacement of Diplograptina by Neograptina occurs during the early phase of relative sea level fall. After being absent from the paleotropics throughout the late Katian ( 14 ), the Neograptina invaded at the onset of the glacial epoch and immediately became the dominant species: Neograptines comprise 90–99% of specimens recovered in this interval at both sites, despite their different depositional settings. At VC, the appearance of the neograptine community is associated with hopane−sterane ratios that suggest the phytoplankton community shifted from cyanobacteria in O4 to dominantly green algae in O5. Concomitantly, observed δ 14 N values suggest reduced denitrification in the local water column ( 21 , 22 , 30 ). During the Hirnantian at the bathyal VC site, 10 diplograptine species reappeared briefly in association with intervals of relative sea level rise, but diplograptines remained minor elements (0–8% of specimens) within communities dominated by Neograptina ( Fig. 1 ; see also ref. 41 ). This pattern suggests that ecological stress on graptolite communities arising from changes in productivity patterns and displacement by the invading Neograptina exceeded a critical tipping point that led to the widespread collapse of the previous diplograptine community structure. The diplograptine species that persisted through the glacial interval as minor constituents of the reorganized neograptine communities were essentially “dead clades walking” ( 42 ). All of them went extinct during the postglacial extinction pulse ( 12 , 13 ).

The late Katian O4b−c interval within the P. pacificus Zone corresponds globally to an interval of relatively warm sea surface temperatures ( 24 ) and climate amelioration that is often referred to as the Boda Event ( 37 ⇓ – 39 ). Deposition of black shales, especially in regions of oceanic upwelling, was widespread ( 22 ), and global graptolite species diversity reached a Late Ordovician peak at this time ( 25 , 40 ). As sea surface temperatures subsequently declined into the Hirnantian glaciation, extinction rates also rose and graptolite species diversity declined ( 13 , 17 ). At the present study sites, most graptolite communities of the early P. pacificus interval are best fit by complex SAD models ( Datasets S4 and S5 ). Samples from the subsequent D. mirus interval at VC exhibit lower evenness, and complex SAD models no longer match the observed patterns. Values of H decline toward the end of the D. mirus interval, and individuals of the remaining mesopelagic and epipelagic species become codominant. All of these signals appear to indicate that the community structure at VC became simpler, largely due to the decline in species diversity and concomitant preferential loss of deep-water taxa, so that there was less vertical niche subdivision during this interval at this site. Community changes at BR were more limited but again result in smaller communities dominated by a few epipelagic species joined by the invading Neograptina. Thus, the basic character of each site (a predominantly mesopelagic, oceanic community at VC and a predominantly epipelagic mixed shelf community at BR) persisted, but the fine structure of community composition became increasingly different during the latest Katian, preglacial interval ( Fig. 1 ). This contrast raises the possibility that ocean conditions above the abyssal VC site responded differently to Late Ordovician climate change than did the shelf seas above the BR site.

The differences in pattern of changing composition between these two sites partly reflects their different ecological settings relative to the graptolite habitats. As the evenness (E) and diversity (H) declined through the latest Katian, the preglacial communities at the midshelf BR site became ever more strongly dominated by a smaller number of epipelagic species as the oceanic water mass moved offshore [indicated by the upsection shallowing recorded in Neodymium isotope ratio (ε Nd ) values] ( 34 ). In contrast, specimen counts from the bathyal VC site reveal a more subdued shift toward faunas, with an increased proportion of specimens attributed to epipelagic species ( Table S3 ). Those shelf-edge communities maintained their more dominantly mesopelagic species composition (60–80%) until nearly the end of the Katian, even as E and H declined through the D. mirus interval. Local depositional conditions also affected the record of faunal response. In particular, the O5a transgressive interval is thicker and more completely sampled at BR ( Fig. 1 ), whereas this interval is strongly condensed at VC ( 34 ). At both sites, however, the available data from the O5a transgressive interval indicate that faunas failed to recover as they did in the preceding transgressive to highstand intervals: H, E, and proportion of mesopelagic species remain like those reached during the relative sea level fall and reach their lowest levels coincident with the O5a minimum in ε Nd values. The correspondence between ε Nd values and community diversity measures (H, E, and proportion of mesopelagic species) in samples from O5 succession are markedly different from those of the preceding O4 samples ( Fig. 1 and Fig. S8 ). These relationships indicate that the graptolite faunas tracked shifting water masses during the Late Ordovician relative sea level changes (as recorded in the changing ε Nd values) and that these relationships altered significantly during the course of the latest Katian environmental change.

Graptolite communities at the study sites also exhibit a significant change in biotope structure over the course of the late Katian. Faunas at both sites exhibit small fluctuations in the proportion of epipelagic and mesopelagic species through the late Katian O4 interval ( Fig. 1 ; 100–54% mesopelagic at VC and 70–44% at BR), roughly in step with relative water depth changes and changes in H and E. Contingency tests ( vassarstats.net/ ) for changing proportions based on specimen counts reveal sharp and highly significant declines in the dominance of mesopelagic species (from 35% to 11%; Table S2 ) and the abundance of mesopelagic specimens over the course of the pre-D. mirus interval at the BR site (from 51.6% to 4.7% mesopelagic; Table S3 ).

Scatter plot of paleoecological data from the (Left) VC and (Right) BR samples, with sample markers color-coded according to whether they are from late Katian O4a−c relative sea level cycles or the Hirnantian O5a−c cycles. Note the different relationships between paleoecological data and ε Nd values between these two intervals at both sites. High (relatively less negative) ε Nd values correspond to carbonates formed under relatively more oceanic (probably deeper) water mass conditions, whereas lower ε Nd values represent sediments formed in more-shelf-dominated water masses ( 34 ).

Scatter plot of paleoecological data from the VC and BR samples, combined. (Left) The proportion of mesopelagic species in each sample versus the Shannon−Weaver information (H). Sample markers are color-coded to correspond to the fitted SAD; below are regression statistics for these data. (Right) Interpolated ε Nd values versus the Shannon−Weaver diversity (H). Samples are color-coded to correspond to the fitted SAD; below are regression statistics for these data. Although percent variance explained by these relationships is modest, they are nonetheless statistically highly significant, suggesting that high proportions of mesopelagic species within samples occur more frequently in samples with high H (indicative of high evenness) than expected at random and, similarly, that more positive ε Nd values (indicative of high relative water depth) occur more frequently in association with high H than expected at random. The ε Nd values are from ref. 34 , and statistical analysis is from vassarstats.net/ , accessed April 19, 2016.

Rarefaction curves for the (A and C) BR and (B and D) VC collections. In A and B, the curve for each sample is identified by meterage and color-coded by biozone. Curves that rise steeply reflect relatively more-even SADs, whereas those that rise slowly reflect more uneven distributions. Those that plateau are relatively completely sampled. In C and D, the same data and curves are show as in A and B, but curves are labeled by the sampled evenness and colored according to the best-fit community type. Although curve shapes reflect both E and community type, there is no significant relationship between community type and sampled diversity among the studied collections.

The community structure statistics for VC and BR ( Fig. 1 ) exhibit peak values of Shannon information (H) and evenness (E) in the Dicellograptus ornatus to mid-Paraorthograptus pacificus zones (see also Datasets S4 and S5 ), during the interval of the O4 relative water depth cycle ( 34 ). The highest recorded H and E values occur in the early P. pacificus Zone interval (O4b) when abundant black shales were deposited beneath a strongly developed denitrification zone at the site ( 21 , 22 ). Biomarkers ( Fig. 1B ) indicate the presence of a predominantly cyanobacterial phytoplankton community in this interval ( 30 ). Graptolite communities at this oceanic site exhibit generally lower H and E during times of lower relative water depth ( Fig. 1 ; see also Figs. S6 – S8 ). Evenness values at both sites are significantly lower during the Diceratograptus mirus interval (mean 0.293 at VC, 0.338 at BR), before the first Hirnantian glacial advance, compared with pre-mirus values (pre-mirus mean 0.4372 at VC, 0.54556 at BR, one-tailed t-Test, df = 9; t = 2.11; P = 0.032 at VC and df = 13, t = 3.15, P = 0.004 at BR). The Aikake Information Criteria based model choice procedure used to determine the species abundance distribution (SAD) model that best describes the graptolite abundance data for each sample ( Dataset S6 ) produced clear assignments to either a simple or a complex community structure for all but one of the sampled VC assemblages. Complex community structures are absent from VC within the D. mirus Subzone. Community structure and evenness reach their lowest levels coincident with the major faunal turnover at the beginning of the Hirnantian Age (base of the Metabolograptus extraordinarius Zone) during the O5a maximum flooding interval.

Stacked plots of the number of species classified as mesopelagic, epipelagic, and indeterminate (colored bands) and the percent of mesopelagic species [N(meso)/[N(meso)+N(epi)]] sample by sample at VC and BR. Data for VC are from Štorch et al. ( 41 ), which includes both the 16 bulk samples for which abundance data have been analyzed here and 30 additional samples with species presence-only data. Temporal spacing between samples is assumed to be the total duration of the set of zones spanned divided by the number of samples; estimates of zone durations are from Cooper et al. ( 94 ).

The species compositions of the late Katian graptolite faunas at the BR and VC sites are very similar. At both sites, a small subset of species are regularly among the most abundant; however, the particular species that dominate each collection varies substantially among successive collections. Furthermore, although nearly all of the species present at BR are also present at VC, the dominance structure of communities differs substantially between the two sites. Epipelagic species are significantly more abundant than mesopelagic species (higher total specimen counts; Table S1 and Fig. 1A ) at BR and are more frequently among the most abundant species in nearly all BR samples, whereas, at VC, it is the mesopelagic species that have the highest total specimen counts. These differences in biotope dominance match expectations from the depositional setting of the two sites. Geological context, lithofacies, and associated fauna ( 19 , 20 , 28 , 29 ) suggest that BR represents a midshelf site, whereas the VC site was on the continental slope to ocean floor. Because the species occurrences at BR were not used in the Bayesian biotope inference process, this result serves as a cross-validation of the biotope assignments themselves. We also examined species dominance patterns at five coeval sites in South China, which also were not used as part of the biotope inference process, and they too exhibit faunal contrasts that match expectations from their geological setting ( Figs. S4 and S5 ).

With the aid of a Bayesian biotope assignment procedure with a 90% confidence level cutoff, we identified 12 late Katian diplograptine species at BR and VC as mesopelagic ( Fig. 1 ) and 9 species as epipelagic ( Dataset S3 ). The remaining 22 diplograptine species had indeterminate biotope affinities. These results were consistent with a detrended correspondence analysis (DCA) based depth assignment of the same species and localities ( Dataset S3 and Fig. S3 ). Neograptine species invade the paleotropics immediately before the start of the Hirnantian Age ( 14 ), and, thus, only four of these species occur within the latest Katian at the study sites ( Dataset S2 ); two emerge as epipelagic, but the other two are too rare to classify.

Graptolite community structure reveals evidence of environmentally induced stress on these organisms before the main pulse of mass extinction during the early Hirnantian. The changes in late Katian graptolite communities reflect the combined effects of species losses due to extinction and ecological reorganization driven by the effects of climate change on the distribution of nutrients, water column oxygenation, and phytoplankton community composition. Even the apparently abrupt change in communities at the start of the Hirnantian was predominantly a consequence of dramatic changes in community structure as opposed to imminent species extinctions, because many of these extinctions occurred hundreds of thousands of years after this ecological turnover event and neograptine invasion. Ecological stresses on communities, therefore, may not be obvious from species diversity patterns alone. Closer examination of community structures, habitat preferences, and abundance patterns provides unique information about the impacts of climate change before the onset of elevated rates of extinction.

The dominantly epipelagic community at BR was consistently less complex and had lower evenness than VC communities, which were relatively even and complex with large mesopelagic components before the D. mirus Subzone interval. The graptolite communities at both study sites experienced significant change during the latest Katian Age, but the different responses at each site led to a divergence in community composition and abundances structures. Furthermore, relative sea level rise during the latest D. mirus interval did not restore the previous community structures. These features indicate that the observed changes in graptolite community structure were not a sampling artifact caused by shifting facies and eustatic sea level fall ( 43 ) and reinforce suggestions that carefully controlled analysis across multiple depositional environments is essential to understand the common cause effects of sea level and environmental changes on mass extinction ( 43 , 44 ).

The picture that emerges of graptolite community structure from these two sites in the interval leading up to the mass extinction is of planktonic communities under stress due to the deteriorating conditions for mesopelagic taxa. Although sea level actually rose near the end of the interval leading up to the glaciation and the mass extinction interval, the oxygen minimum zone contracted and the supply of nutrients from denitrified deep ocean water decreased as a consequence of changes in deep-water ocean circulation patterns. These changes (and probably others, most notably lower sea surface temperature and increased seasonality) stressed graptolite communities, leading to reduced community complexity and evenness at both oceanic and midshelf sites. The species that most strongly declined, in terms of their individual contributions to community composition, were preferentially those adapted to the mesopelagic habitat. Brachiopod species reveal a similar selectivity ( 18 ). Extinction of the majority of the graptolite species making up these stressed communities did not occur, however, until the onset of the glacial interval (or even later), when diplograptines were replaced by the neograptine species that migrated into the paleotropics and diversified during the Hirnantian Age. This pattern of selectivity may have contributed to the unusually high losses of old, well-established graptolite species (positive age-dependent risk of extinction) during the LOME ( 17 ).

Additionally, we fit a series of SAD models (geometric, log series, and lognormal) to the data. The lognormal model differs from the geometric and log-series models in the number of species of intermediate abundances, and may be indicative of relatively complex communities with substantial species interactions and species interdependence ( 6 ). The geometric and log-series models have been interpreted as indicating relatively simple community structure, without extensive interaction among the species within the community. AIC model ranking ( 46 ) was used to determine which of the fit models should be considered as viable at each horizon ( Dataset S6 ). We considered horizons to have a complex distribution if the weight of the lognormal model was above 90% and to have simple distribution if the combined weight of the simple distributions (geometric and log series) was above 90%. All other horizons are indeterminate.

We inferred biotope affinity through a novel Bayesian approach. We treated each collection horizon as an individual case, so that there were 20 shallow collections, 37 deep collections, 25 intermediate collections, and, excluding taxa seen in only one section, a total of 43 species in the analysis. The basic approach is to build a maximum likelihood model of finding each species at a shallow or deep horizon, given the observed patterns of presence and absence at each horizon at each locality and the biotope affinity of each species (see SI Text for analytical details). A Markov chain Monte Carlo (MCMC) search process was used to determine the posterior probability distribution of biotope affinity, starting from a prior probability of 0.5 for the biotope affinities of each species. This statistical modeling procedure assigns species with strong evidence of restriction to deep water a posterior probability of being mesopelagic that is close to 1, whereas those that appear regularly in shallow water sites will have probabilities of being mesopelagic that are close to 0 (and vice versa). Intermediate values indicate contradictory or limited evidence. We categorized as mesopelagic those species that received a Bayesian posterior probability of being mesopelagic of >90%, and classed as epipelagic those with a mesopelagic probability of <10%; all others are indeterminate ( Dataset S2 ). We chose a 10% confidence level because of the limited sample sizes. Twenty-three of the 40 Katian species present at VC and BR were confidently classified with this protocol (9 epipelagic and 14 mesopelagic) with the other 17 left in the indeterminate category ( Datasets S3 and S4 ). In light of the very clear segregation of our sites and taxa into binary categories (shallow/epipelagic and deep/mesopelagic), we used DCA ( 45 ) to analyze the same set of data that we used in the Bayesian analysis as a means to assess the plausibility of the outcome of this new procedure.

We used statistical procedures to estimate the biotope affinity (mesopelagic or epipelagic) of each of the species of Diplograptina present within late Katian strata at VC and BR. In addition to the species abundance data available for VC, we also used horizon-by-horizon records of species presence/absence from the P. pacificus Zone at 18 other localities to provide a strong contrast between deep-water and shallow water localities ( Fig. S2 and Dataset S1 ). The Trail Creek, Idaho, and Dob’s Linn, southern Scotland, sites provide additional data from deep oceanic settings, whereas sites at Truro Island, Eleanor Lake, and Cape Manning in the Canadian Arctic provide species occurrences from relatively shallow, midshelf to platform margin settings (ref. 22 and Dataset S1 ). A set of five sections in the Omulev Mountains of Siberia and six sections in the Chu-Ili Range of Kazakhstan expand the sample set beyond Laurentia and encompass a range of site conditions from deep back-arc to midshelf ( Dataset S1 ). We used data from a set of five sections from South China to cross-validate the biotope assignments ( SI Text ).

We gathered species abundance data from bulk samples at the well-studied, deep-water Vinini Formation in central Nevada. All identifiable specimens of each species were counted, resulting in collection sizes from 166 to 606 identified specimens per horizon ( Dataset S3 ). Bulk samples at the BR site in the Ogilvie Mountains of Yukon (Canada) were larger than at VC, and ranged from 124 to 3,072 identified specimens per sample. We preferentially sampled well-laminated, fine-grained sediments throughout both sections, seeking, as much as possible, to hold preservational and facies-related variance in assemblage composition to a minimum. Variation in ε Nd values through the sections indicates that most samples represent the deeper parts of the series of relative water depth cycles present through the studied interval, adding to their taphonomic consistency.

SI Text

Graptolite Paleoecology. An understanding of the paleoecology of graptolites forms a critical predicate to the present study. It underlies the basic hypothesis that we have sought to test and the methods that we devised to test our hypotheses. Because this background is not likely to be common knowledge among the biological, or even paleobiological, community, we provide a brief summary of that consensus here and describe its connection to our approach to community change during the lead-up to the LOME. This summary is based largely on the Cooper et al. (19) synthesis and the works cited therein. A considerable body of observational data on planktic communities in the modern oceans together with data on the occurrence of graptolite species in the rock record suggests that planktic graptolites dwelt primarily over the outer shelf and adjacent shelf edge. Some species were endemic to particular shelves or epicontinental seas on particular continents, but most were relatively widely distributed—some were cosmopolitan within the tropical or subtropical realms that they inhabited. These shelf endemics and many of their more nearly cosmopolitan relatives occur in shelf sediments deposited above storm wave base together with shelly benthic fossil assemblages (or at least are interbedded with them—preservational conditions tend to affect graptolites differently than calcareous macrofauna). These species then appear to have dwelt in the surface mixed layer of the ocean at relatively shallow water depths, with high sunlight and oxygen levels: the epipelagic biotope (or epipelagic zone; Fig. S1). These near-surface dwellers are referred to as epipelagic species, and we refer to this habitat preference as “epipelagy.” Along the shelf edge in many regions around the globe, wind or current-driven upwelling leads to enhanced surface productivity in the phytoplankton and their various consumers in response to the concomitant resupply of otherwise limiting nutrients. This surface productivity also supplies enhanced fluxes of particulate and dissolved organic matter to the water beneath the epipelagic zone, where it may stimulate high productivity among detritivores and decomposers and, as decomposition consumes the oxygen available in the upwelling currents, denitrifying bacteria. Many graptolites appear to have been specialist consumers within this twilight, or mesopelagic zone, beneath the surface mixed layer in regions of periodic upwelling or sustained high surface productivity (Fig. S1). We refer to this habitat preference among graptolites as “mesopelagy.” The physical water mass boundary at the base of the surface mixed zone (which is typically associated with rapid changes in salinity and temperature, among other properties, leading to a sharp density gradient, the pycnocline) and the shift from photosynthetic to decomposition-based food chains most likely divided graptolites into two distinct adaptive groups: the epipelagic biotope and mesopelagic biotope communities. This situation contrasts with the adaptive array of species one commonly finds along extended ecological gradients, such as those that organize marine benthic animal or terrestrial plant communities (see ref. 19 for more detailed comparison with modern analogs). For taphonomic reasons, graptolites are recovered most commonly from laminated, organic-rich mud rocks; however, they occur in a very wide range of shelf and deep oceanic facies, as one expects from pelagic organisms. Sediments deposited in shelf settings within the epipelagic zone of the oceans primarily receive graptolites settled from the overlying waters, and, accordingly, they record a fauna that is strongly dominated by epipelagic species. Deep oceanic facies, on the other hand, include both epipelagic and mesopelagic species in substantial numbers. Thus, epipelagic species occur in both onshore and offshore settings. Only the mesopelagic species exhibit any significant degree of facies limitation, and only in that the depositional site lies within the mesopelagic zone or deeper. Mesopelagic species probably were occasionally transported onto the shelf by currents before burial, and so may occur rarely in assemblages that otherwise would contain only epipelagic species. Relative to the epipelagic species, the mesopelagic species generally had more limited geographic distributions, shorter species durations, and higher risk of extinction (17, 20, 47). Disentangling the cause and effect relations among these three properties is difficult but has been addressed to some degree in the studies cited above. Considerable further work is needed on this subject, however. As a result of these ecological and sampling issues, it is clear that the determination of epipelagy or mesopelagy from graptolite species occurrences is an asymmetrical inference process that is dependent on the absence of species from depositional sites regarded a priori as within the epipelagic zone of the oceans. As such, reliable conclusions depend critically on the quality of collections and prior knowledge about depositional conditions, and we herein take a conservative Bayesian approach to this process that allows us to explicitly take account of the particular requirements of graptolite community analysis.

Collection and Locality Data. Material for this study is a combination of material collected expressly for the purpose of analyzing species relative abundance patterns at VC (41, 48), and BR (21, 32), and previous data collected for biostratigraphic analysis at sites 2 and 4–10 that are suitable for documenting species’ presence or absence but lack abundance information. Information on the present and paleogeographic locations of the study sites is summarized in Dataset S1 and Fig. S2. Data from Trail Creek were previously collected by several others (49⇓⇓⇓–53) and restudied by C.E.M.; those from Truro Island, Eleanor Lake, and Cape Manning were gathered by M.J.M. (54, 55) (M.J.M. additional collections as shown in Dataset S2); data from Dob’s Linn are based on restudy of material collected by Williams and coworkers (56⇓–58) and new collections by M.J.M. and J.L. (data included in Dataset S2). Data were updated from the literature for five sections in the Omulev uplift region, Kolyma, Siberia (59, 60), six sections in the Chu-Ili Range, Kazakhstan (61, 62), and five sections from the Yangtze Platform region, South China (63⇓⇓–66). Systematic treatment of all of the collections follows that of Štorch et al. (41), and the complete data used are included in Dataset S2. We used presence/absence data from the P. pacificus Zone interval at sites 1, 2, and 4–9 (Fig. S2 and Datasets S1 and S2) in both the Bayesian and DCA approaches to depth affinity estimation. Sites were coded as shallow or deep based on previous geological studies that determined them as midshelf to inner shelf (shallow) vs. slope-rise or back arc settings (deep), based on a combination of sedimentological features (e.g., common storm- and wave-generated sedimentary structures vs. debris flow and turbidite/contourite deposits), regional basin reconstructions, and tectonic studies. We omitted sites that were poorly constrained or that occupied intermediate depositional settings. This set treated each sample as an independent data point and, thus, comprises data on 43 species in 82 collections. Presence/absence data from the Yangtze Platform sites (128 collections) were combined with those from BR (18 collections) in a cross-validation analysis of the Bayesian biofacies inference procedure. Upon conclusion of the Bayesian biofacies inference procedure, we turned to a detailed analysis of the biotope and species abundance structure of graptolite communities through the late Katian at VC and BR based on counts of specimens from 34 bulk samples at those two sites. That graptolites and other organisms experienced mass extinction during the Hirnantian glacial episode and its immediate aftermath has been demonstrated by many previous studies (13, 22, 44, 67, 68) and is not itself the subject of this study. Rather, our goal in this last phase of the present study was to use inferences about species habitat preference to examine whether these two sites record temporal changes in graptolite community composition and selective loss of mesopelagic species during the period of Late Ordovician climate deterioration that led up to the mass extinctions. The VC site has received considerable previous study (summarized in refs. 22, 34, and 41). The interbedded carbonate mudstones, grain flows, and siliceous mudstones are generally laminated and contain no benthic fauna. In addition to graptolites, the rocks yield conodonts, various organic-walled microfossils, sponge spicules, and radiolarians. The rocks are part of the upper plate of the Roberts Mountains allochthon, and appear to have been deposited at abyssal depths in an outer slope and rise setting within a back-arc basin system before being transported onto the Panthalassic shelf of Laurentia during the Antler Orogeny. Thus, eustatic sea level fluctuations during the late Katian and Hirnantian were probably small relative to total local water depth at this site. The principle facies changes are limited to changes in the frequency of grain flows, which decrease upward. The rocks remain well laminated and continue to lack benthic fossil burrows. There are no obvious omission or erosion surfaces within the studied portion of the Vinini Formation, although one might expect to find such associated with the lowstand fans that ought to have accompanied the glacial low stands during the Hirnantian interval. Thus, we see no evidence of facies change or hiatuses in the VC section of a sort that should significantly alter the preserved graptolite community. Instead, observed changes should primarily reflect changing ecological conditions in the overlying water column. The BR site records a markedly farther upslope depositional setting. In contrast to the VC site, the BR section does record significant changes in bottom facies, including several gradational alternations between laminated, graptolite-rich mud rocks and argillaceous carbonate mudstones. Benthic fauna (trilobites and brachiopods) are also sometimes present, in both the carbonate mudstones and black shales, but, with the exception of abundant trilobite material that occurs at the base of the Hirnantian, this is rare within the study interval. Although most strata from the pre-Hirnantian produced graptolites, the carbonate mudstones were difficult to sample adequately, and, thus, our samples focus on the laminated mudstone intervals that occur intermittently (primarily in the relatively deeper water intervals as indicated by the ε Nd data) throughout the late Katian. Thus, the BR samples track a consistent depositional facies through the studied interval. Eustatic sea level fall during the early Hirnantian extinguished this laminated mudstone facies, however. As a result, our sampling for this study ended very low in the Hirnantian, Metabolograptus extraordinarius Zone interval. Graptolites return in abundance in the postglacial interval at this site, whereas that interval is marked by a starvation surface (base of the Elder Formation) at the deep-water VC site. The previous data on the carbon isotopic composition of organic matter and bulk carbonate at these sites (ref. 21 and Fig. 1) provide independent means to identify the onset of the Hirnantian, which, globally, is marked by a prominent positive carbon isotopic excursion (22). In addition, Holmden et al. (34) reported variation in Nd concentration and ε Nd values that they interpreted to reflect changes in relative water depth (Fig. 1). As a result of mixing between continental and midocean ridge sources, shelf seas may exhibit a gradient from moderately negative ε Nd values in oceanic sites to strongly negative values in more interior sites. Accordingly, stratigraphic change in ε Nd values may reflect the changing composition of the overlying water mass location in response to increase or decrease of relative water depth at the site. The broadly synchronous pattern of changing ε Nd values between VC and BR, despite their abyssal and shelf setting, suggests that the shared component of these patterns tracked Katian to Hirnantian eustatic sea level cycles (labeled O4a−O6b in Fig. 1). This interpretation is supported by the very similar pattern and timing of sea level change recorded in the sequence stratigraphic architecture of the late Katian to Hirnantian Anticosti Island succession on the other side of the continent (details presented in refs. 22, 34, and 69). Thus, we conclude that the VC and BR sites are appropriate sites from which to address Late Ordovician community change among graptolites. The sites fit well the conditions described by Holland and Patzkowsky (43) as essential to discriminate true species turnover patterns from artifacts driven by sequence stratigraphic architecture. In particular, their different positions relative to the edge of the continental shelf suggest that facies change at these sites may be expected to respond differently to eustatic sea level. Consequently, shared patterns of community change through the late Katian are likely to reflect shared regional environmental features in the overlying water column rather than local facies-driven artifacts of species distribution.

Depth Affinity Estimation. Here we describe an approach to estimating biotope preferences of species (or other taxa) from their observed pattern of occurrence (biofacies). Probabilistic modeling of habitat preference has been done in both Bayesian and frequentist terms, using binomial models of habitat preference (70⇓⇓–73) as either the conditional probability of the data given the prior probability (71, 72), as a way of computing relative probabilities in frequentist form (70), or as a null hypothesis of equal preference (73) that could be rejected to assign a genus or species to one of two competing habitat groupings. These binomial models are based on the premise that we have some large number of collections of data, of which a portion p 1 are of type 1, and the remaining proportion p 2 = 1 – p 1 are of type 2. If there are n observations total, with n 1 in habitat 1, and n 2 = N − n 1 in habitat of type 2, then, under a hypothesis of no habitat preference, the binomial probability of n 1 observations of type 1 can be determined. This can be used as a Bayesian conditional probability under a prior of equal habitat preferences, or it may be used as a hypothesis test directly. Clearly, larger sample sizes (N) will allow for stronger probability statements. The epipelagic and mesopelagic categories within graptolites are somewhat different from a simple binary preference for one habitat over another, in that epipelagics are expected in both environments whereas mesopelagics are confined to deep-water sites. In the context of the binary models discussed above, one could simply consider epipelagics as species with no habitat preference, and mesopelagics as those with strong evidence for deep-water habitat preference. The approach developed and used here, however, is derived from the maximum likelihood methods used in capture−mark−recapture (CMR) study of wildlife and applied to biodiversity patterns in the fossil record (13, 74, 75). In this conceptualization of the problem, we consider that we have collected some number of samples at a series of locations (horizons within sections) of differing habitat types. At each site, there is a probability of finding each species, provided that the site is located within the temporal duration of the species. This corresponds to a sighting or recovery probability in the CMR framework. Each species has a sighting probability, which may be conditional on its habitat preference. Given the pattern of observance and nonobservance of the species at a number of locations of differing types, the probability (or likelihood) of the observed pattern may be calculated under the different competing hypotheses about habitat preference. Here, we used the competing hypotheses of mesopelagy and epipelagy, but the method could be applied to simpler or more complex patterns of habitat preferences. As will be discussed below, the method could also be used in cases in which the prior hypotheses about the nature of the habitats are not conclusive, with the prior probability of sites being deep being allowed to vary during the fitting process and a posterior probability of site condition as an outcome of the analysis. This approach is based on a series of likelihood functions that describe the probability of finding a particular species at a specific site, given that the species is either E, epipelagic (dependent on the photic zone), or M, mesopelagic (dependent on conditions below the surface mixed layer). The model may also be extended to include species that are restricted to shelf habitats [OS, obligate shallow, Cooper and Sadler’s (20) Group 3 taxa]. For the present study, however, the simpler two-category model was used. Each site is classified as either deep (D) or shallow (S), and, thus, this approach can be used to produce joint estimates of depth affinities of the species and depth category of the sites. This approach also allows for the incorporation of prior information about both sites and taxa into the procedure. In this model, each taxon has a probability of being collected in a deep-water or shallow water site that is dependent on the classification of the site and the taxon. The following variables are defined: N(D), total number of deep horizons at which the species could be recovered given that the species duration overlaps the time of deposition of the horizon—this value does depend on the species in question; N(S), total number of shallow sites horizons at which the species could be recovered given that the species duration overlaps the time of deposition of the horizon—this value does depend on the species in question; N i (D), number of deep sites at which the ith taxon is recovered/observed; N i (S), number of shallow sites at which the ith taxon is recovered; P i (D|E), probability that the ith taxon is seen at a deep site, given that the taxon is epipelagic; P i (S|E), probability that the ith taxon is seen at a shallow site, given that the taxon is epipelagic; P i (D|M), probability that the ith taxon is seen at a deep site, given that the taxon is mesopelagic; and P i (S|M), probability that the ith taxon is seen at a shallow site, given that the taxon is mesopelagic. Based on the definition of mesopelagic, it might seem that P i (S|M) would be zero. However, zero-probability events pose difficulties for numerical algorithms, as well as being unrealistic in geological terms, as rare appearances of mesopelagic species at shallow sites may occur by transport or as a result of a misidentification. Setting P i (S|M) = P e , where P e is on the order 0.001 or less, allows for some probability of identification errors of taxa or transport of a specimen from a different environment (or both), as well as allowing for smoother numerical solutions. In most situations, it is wise to run a sensitivity analysis on the value of P e to ensure that the answer is not strongly dependent of the value used for P e . A further simplification is to assume that epipelagics are recovered with equal probability in both deep and shallow sites P i (S|E) = P i (D|E) = P i (E). We can then write a likelihood function for each taxon under each possible assumption of its biotope. The likelihood function is used here in a Bayesian context, although it could also be used for maximum likelihood modeling. For the ith taxon, the probability of the observed pattern of recovery and nonrecovery at shallow and deep sites, we get four terms corresponding to being found at deep sites, not being found at deep sites, being found at shallow sites, and not being found at shallow sites. There are two different likelihood functions for each taxon, depending on in which of the two categories the taxon belongs. L i ( M ) = P i ( D | M ) Ni ( D ) × ( 1 - P i ( D | M ) ) ( N ( D ) -Ni ( D ) ) × P e ( Ni ( s ) ) × ( 1 - P e ) ( N ( S ) -Ni ( S ) ) L i ( E ) = P i ( E ) Ni ( D ) × ( 1 - P i ( E ) ) ( N ( D ) -Ni ( D ) ) × P i ( E ) ( Ni ( s ) ) × ( 1- P i ( E ) ) ( N ( S ) -Ni ( S ) ) . The total likelihood of all taxa at all sites is the product of these terms over all species. It is common to work with the log likelihood, as the terms in the total log likelihood then become a sum over the log likelihood of each taxon, which greatly simplifies the calculations. Probability values that maximize the summed log likelihood will also maximize the likelihood function itself. The total likelihood function is then dependent on the assignment of each species to one of the two biotope categories (M, E) and of each locality to one of two categories (D, S), as well as the fitted value of each P i . Given any set of biotope assignments for the species and the set of depth assignments for each locality, the likelihood of the observed pattern of presence/absence data may be calculated. It turns out that the most obvious estimates of the P i values do serve to maximize the likelihood and log likelihood functions, P i ( D | E ) = P i ( S | E ) = ( total sites recovered ) / ( all possible sites ) = ( N i ( S ) + N i ( D ) ) / ( N ( S ) + N ( D ) ) P i ( D | E ) = ( number of deep sites recovered ) / ( all deep sites ) = ( N i ( D ) ) / ( N ( D ) ) . The log likelihood function obtained by summing the log likelihoods of all of the taxa may be used to obtain either maximum likelihood or Bayesian estimates of each species’ biotope affinity, given the site categories. It also may be used to obtain Bayesian estimates of the site classifications given the species biotope affinity or to obtain joint estimates of both species and site classifications. In the maximum likelihood approach, the site categories are taken as givens, and the probability is calculated for each taxon under each possible assumption of biotope affinity. The species is then assigned the biotope affinity that maximizes the log likelihood. For each species, it is possible to form likelihood ratios for the different biotope assignments as one way of determining the relative strength of evidence for that biotope assignment presented by the data. In the current study, we regarded the assignment of site as deep or shallow as robust, and thus we treated these as known constant values, and the species affinities were regarded as unknown. We chose to use a Bayesian procedure based on flat species prior probabilities in which each species was equally likely to be initially assigned to each of the two biotopes. An automated optimization approach, MCMC (76, 77), was then used to arrive at Bayesian posterior estimates of the probability of each species being in each biotope. In the MCMC estimation procedure, using the Metropolis−Hastings algorithm (76, 77), a species is randomly chosen and randomly assigned a different proposed biotope. The change in the likelihood function is calculated. If the proposed biotope change results in a higher likelihood solution, then the proposed change is selected. If the proposed biotope change decreases the likelihood of the solution, the proposed assignment is accepted if a random number generator returns a value lower than the fractional decrease in the likelihood. The MCMC procedure is “burned in” for some relatively large number of trials, perhaps 100 times the number of species used, so that each species is changed 100 times on average. After the burn-in procedure (used so that the search algorithm is working at or near the optimal solution), the MCMC process is allowed to run for hundreds of additional trials per species, recording the species at some large fixed number of changes. The percentage of time over all MCMC trials the species is observed to be in a given category (E, M) is then interpreted as a Bayesian posterior estimate of the probability of that species belonging to each biotope. The Bayesian approach has the advantage of not simply assigning the biotope but also estimating the probability of that assignment. The larger the amount of data supporting the assignment, the higher the Bayesian posterior. Species with few occurrences will not have high Bayesian posteriors for any biotope, as they have little impact on the overall likelihood. In our analysis of the present dataset, it appears that species need to be observed at five or six collection horizons to have any chance to a reach 90% posterior probability of being in a particular biotope. These methods can also be extended to situations where some or all of the site depth assignments are not well known. In the MCMC search, the algorithm can also be allowed to randomly alter the site depth assignments (starting with some prior probability for each site). This approach would allow for simultaneous estimation of categories for both sites and taxa. The dual estimation approach was also checked for this dataset and did not produce substantially different results for the a posteriori species biotope affinity estimates. The results of the Bayesian modeling are shown in Dataset S3 for the taxa used as inputs to the estimation. This procedure resulted in confident biofacies estimates for 23 of the 43 species included in the analysis: 12 as mesopelagic and 9 as epipelagic. The other 22 could not be confidently classified by this approach. Examination of the number of horizons at which the indeterminate species were observed (Dataset S2) indicates that most of the indeterminate species were rare (seen at fewer than five horizons). Thus, the Bayesian approach indicates that the limited information available about these species’ occurrence means that they cannot be reliably assigned a depth affinity with a high level of confidence. We discuss these results in more detail below in comparison with results from DCA and from application of a simple binomial model using epipelagy as a null model.

A Binomial Model Using Epipelagy as the Null. We also applied the binomial model of habitat preference as used by Foote (73) to this dataset (Dataset S2). In this model, epipelagy implies no preference for deep sites over shallow sites, so that epipelagy becomes the null model. The probability of observing a species in k or more deep sites out of N trials when the proportion of deep sites is p can then be calculated. If this probability is low (below 5%), then the species must be assumed to be mesopelagic; if not, then the null of epipelagy cannot be rejected. Note that this binomial model does result in an acceptance of the null hypothesis, because the alternative to mesopelagy is epipelagy rather than no preference among sites. Thus, the binomial model classed as epipelagic several species for which we had relatively few observations, which the Bayesian approach left as indeterminate. This test also assigns specimens to mesopelagic if there is a statistically significant departure from equal preference for shallow and deep sites, whereas, in the more complex Bayesian model discussed above, the mesopelagics are expected to be largely restricted to deep sites, rather than simply preferring deeps. Thus, the Bayesian model has a somewhat stricter requirement and consequently found fewer taxa to be likely to be mesopelagic. Thus, the binomial test is less effective in the current setting than in situations where there may be a more appropriate null of truly no facies preference. DCA. We assessed the outcome of the Bayesian biofacies inference approach by two independent methods: DCA and cross-validation. The DCA ordination method has been criticized on theoretical grounds (78) but has also been used extensively for determining depth gradients of fossil species (79). An attractive feature of the method is that attempts to place localities and species in a common ordination, without any a priori assumption about the nature of the cline or gradient in the data other than that the great preponderance of the variance of interest in the dataset, can be accounted for by ordination along a single factor (45)—in our case, the epipelagic to mesopelagic ecological gradient. We conducted the analyses in Past 3.0 (80), which is freely available for download (folk.uio.no/ohammer/past/). The analysis of species occurrence patterns among deep to shallow sites arrayed the sites primarily along the first DCA axis (Fig. S3A), but also separated some deep sites along the second DCA axis based on their unique species occurrence patterns. A DCA Axis 1 score of 400 correctly separates 90% of the 77 samples that were confidently categorized by our Bayesian approach as either shallow or deep. Samples left as indeterminate by the Bayesian approach straddle that cutoff but fall mainly among the Bayesian shallow samples. Thus, the DCA results and those of our novel Bayesian method appear to be capturing similar information. In particular, the DCA analysis results, which had no a priori information about site properties or species biotope affiliation, unlike the Bayesian approach, agree well with the site conditions derived from published information on the geological setting and sedimentological characteristics of the sites (22). In the present results, species plot in a somewhat different portion of the DCA ordination space than do sites (Fig. S3). Overall, they occupy a narrower range, and their mean score is shifted to smaller values (272 vs. 354), but, again, high values on Axis 1 correlate with epipelagic status whereas low values suggest mesopelagic affinities. All but five taxa, however, have Axis 1 scores below 400 and so are associated in DCA ordination space with sites classed as deep by the Bayesian method and prior geological knowledge. On the other hand, a score of 318 on DCA Axis 1 (400 minus the difference in mean scores) separates the 23 Bayesian epipelagic and mesopelagic species with only one epipelagic species misclassified (96% correct) (Fig. S3B and Dataset S2). Mesopelagic species of graptolites are expected to be confined to deep-water sites, but most epipelagic species are expected to be present at both deep and shallow sites, and only shelf endemic species are thought to have been restricted to shallow sites (19, 20). DCA scores for species reflect the mean affinity of species, so that we might expect many epipelagic species to have intermediate scores, reflecting equal affinity to both deep and shallow sites, and thus to not be cleanly identified by the DCA. Furthermore, taxa that are unique to particular sites tend to have extreme values on one or the other of the first two DCA axes, which indicates that the ordination is also influenced by biogeographic differences in faunal makeup, i.e., that the species composition of samples is not a simple product of site depth alone. Although the DCA placed many species in the mesopelagic field that remained indeterminate by our Bayesian approach, we regard the latter as a more conservative approach that, unlike DCA, is supported by estimates of the reliability of the biofacies inference. In addition to these advantageous statistical features, the Bayesian approach appears to be more effective at disentangling biotope from regional faunal differences than DCA. This difference probably reflects the incorporation of a specific biofacies model, the method’s grounding in a probabilistic framework, and its ability to include prior information about site properties. Cross-validation of the Bayesian depth inference results. In addition to comparing results of DCA and Bayesian biotope classifications, which both use the same data, we have also taken a cross-validation approach to assessing our modeling results. The collections from BR and the Yangtze Platform region were not used in the depth affinity estimation discussed above, but many of the species present at these sites are shared with those at the training sites, and therefore, can be assigned based on the Bayesian results. Armed with those data we can determine the biotope composition of the shared fauna at BR, and the five Yangtze Platform sites. More specifically, prior geological, sedimentological and biofacies analyses at these sites allows us to independently classify each of them as relatively deep (and, therefore, expected to have a diverse fauna with many mesopelagic graptolites) or as relatively shallow (and, therefore, expected to have lower diversity faunas more strongly dominated by epipelagic species). Plots of the number of mesopelagic, epipelagic, and indeterminate species at VC and BR (Fig. S4) show precisely the expected difference: mesopelagic species make up a smaller fraction of the overall less diverse fauna (see discussion in the main text and Dataset 1). Among the Yangtze Platform region, graptolite faunas from five sections that span the late Katian to early Rhuddanian have been studied in detail (Dataset S1 and Fig. S5). Among these the Hirnantian stratotype section at Wangjiawan-North, its very nearby counterpart, Wangjiawan-South, and the Fenxiang section are from near the center of the basin (81, 82). The data for these sections are species presence/absence only. Despite this and the fact that a relatively high proportion of the Yangtze Katian graptolite fauna is endemic, again the fauna at these sections is strongly dominated by mesopelagic species. The near shore sites at Honghuayuan and Ludiping show the predicted shift to faunas with lower proportions of mesopelagic species and overall lower species diversity. These results provide an independent gauge of the applicability and validity of the Bayesian biotope assignments made from the set of training sites. Were they artifacts of the particular local depth assignments of the training sites, we would not expect to see this coherent and predictable depth-related difference in proportion of mesopelagic species at the test sites. Recent geochemical studies of the laminated and carbon-rich Wufeng Shale, which was deposited in the central Yangtze Platform region during the late Katian interval examined here, indicate that the local environment was persistently dysaerobic to anaerobic (28, 83). Thus, the occurrence of mesopelagic dominated faunas at these central Yangtze Platform sites also further strengthens the association between mesopelagic graptolite species and strongly developed denitrification zones rather than with upwelling zones per se because this was a partially restricted deep basin rather than a shelf-edge upwelling site (83). Additionally, it is noteworthy that all these sites, including the main study sites at VC and BR, as well as the Yangtze Platform sites, where we have a long, well sampled graptolite record, exhibit a consistent pattern of changing biotope dominance: deep-water and midshelf sites show declining proportions of mesopelagic species through the late Katian and early Hirnantian, interrupted by a short recovery associated with the sea level rise at the base of the D. mirus Subzone, whereas inner shelf sites show no significant change in the already low proportion of mesopelagics.

Simple Biodiversity Measures. For each horizon at VC and BR, we calculated a series of simple biodiversity measures (Dataset S3). These measures are the classic Shannon information H (84), H = −Σf i × ln(f i ), as well as the raw species count S and an evenness measure derived from H and S, E= eH/S. This measure of evenness is also known as Pielou’s J (e.g., ref. 85). Finally, we determined the ratio of deep to sallow species in each sample based on their biotope affinities, inferred via the Bayesian posterior probabilities described above.

Raw Species Counts and Rarefaction. Rarefaction curves (Fig. S5) are a simple and especially clear means to examine how complete the sampling of a community is. They also serve to illustrate the interaction of community evenness and species relative abundance with sample size. If a community has been well sampled, the rarefaction curve will exhibit a plateau, indicating that new species are unlikely to be found if further sampling is carried out. When two or more communities with unequal sample size are compared, it is often advantageous to compare their respective rarefaction curves to assess the degree to which each of the communities has been thoroughly sampled. Additionally, examination of the shape of a resampling curve is also informative about the evenness of the community. Samples of communities with high evenness exhibit a more rapid rise toward the plateau than do samples of less even communities, because high evenness exists when species abundances are similar, and all species, therefore, should be well represented at relatively low sample sizes. When evenness is low, several species dominate, and many are rare. Consequently, the rarefaction curve rises slowly toward its plateau, as most collected specimens will be from the limited number of dominant species. Detecting the rarer species in low-evenness communities may require a large sample size. Thus, the shapes of the rarefaction curves reflect the same population properties that are captured by the simple biodiversity metrics tabulated in Datasets S4 and S5 and plotted in Fig. 1. Species relative abundance distributions also affect the shape of rarefaction curves along with evenness, because, like evenness, these distributions describe differences in the number and commonness of different species within communities, especially in the number and rank abundance of the many uncommon but not truly rare species that are members of a community (e.g., ref. 86). As we describe below, the graptolite communities studied here exhibit both complex and simple SADs. These differences are also manifest in the differing rarefaction curve shapes (Fig. S6) present among our samples. We emphasize here two features of these curves. First, with a few exceptions among the relatively low evenness communities recovered at BR and VC, most of the curves clearly reach or nearly reach plateaus. Secondly, all have relatively low species diversities (both VC and BR average 12 in the Katian interval examined), despite their comparatively large sample sizes (most are in the range of 300–2,000 identified specimens; Datasets S4 and S5). Apart from the incompletely sampled, low-evenness communities, the Katian graptolite samples examined here are unlikely to gain more than a rare species or two by further sampling. Finally, we note that our analyses rely on relative abundance of epipelagic and mesopelagic species and specimens rather than on species diversity counts per se, and, as a result, we conclude that the available samples are adequate for the purpose at hand.