Biogeographic and bathymetric determinants of brachiopod extinction and survival during the Late Ordovician mass extinction

The Late Ordovician mass extinction (LOME) coincided with dramatic climate changes, but there are numerous ways in which these changes could have driven marine extinctions. We use a palaeobiogeographic database of rhynchonelliform brachiopods to examine the selectivity of Late Ordovician–Early Silurian genus extinctions and evaluate which extinction drivers are best supported by the data. The first (latest Katian) pulse of the LOME preferentially affected genera restricted to deeper waters or to relatively narrow (less than 35°) palaeolatitudinal ranges. This pattern is only observed in the latest Katian, suggesting that it reflects drivers unique to this interval. Extinction of exclusively deeper-water genera implies that changes in water mass properties such as dissolved oxygen content played an important role. Extinction of genera with narrow latitudinal ranges suggests that interactions between shifting climate zones and palaeobiogeography may also have been important. We test the latter hypothesis by estimating whether each genus would have been able to track habitats within its thermal tolerance range during the greenhouse–icehouse climate transition. Models including these estimates are favoured over alternative models. We argue that the LOME, long regarded as non-selective, is highly selective along biogeographic and bathymetric axes that are not closely correlated with taxonomic identity.


Introduction
The Late Ordovician mass extinction (LOME) was one of the largest extinctions of the past 500 million years [1,2], involving the extinction of almost half of marine invertebrate genera and an estimated approximately 85% of species [3]. Major extinction pulses occurred at the boundary between the Katian and Hirnantian stages, and in the mid-Hirnantian (at the base of the Metabolograptus persculptus graptolite Biozone) [3,4]. Both these boundaries are associated with major climatic and oceanographic transitions: the former with global cooling, expansion of south polar ice sheets and falling sea levels, and the latter with warming, melting of ice sheets and continental flooding [3][4][5][6][7][8][9].
Because of the apparent correspondence between climatic transitions and extinction pulses, the LOME has long been recognized as a climatically driven event [3][4][5][6][7][8]10]. However, substantial uncertainty remains regarding the mechanisms of extinction-there are many ways in which climatic transitions and associated oceanographic changes could potentially lead to elevated extinctions, particularly through habitat loss [5,8,11]. Previous studies of the LOME have implicated two major drivers of extinction: (i) loss of species inhabiting shallow cratonic seaways that drained as Gondwanan glaciers grew; and (ii) loss of species with narrow and/or relatively warm thermal tolerance ranges as the polar front advanced and the latitudinal temperature gradient steepened [12,13]. Late Ordovician -Early Silurian climatic events seem also to have been associated with major changes in oceanographic circulation, productivity and oxygenation of outer shelf and slope settings [6,14]. Changes in oxygenation, in particular, have been suggested as a major additional agent of extinction, particularly for taxa with planktonic life stages [15 -17]. Finally, the fact that the Hirnantian stage is represented by hiatuses, lacunae or major facies shifts in many regions [6,13,18 -20] and the preponderance of 'Lazarus' genera that are missing from the Hirnantian and earliest Silurian records, but later reappear [21], suggest that some proportion of apparent losses during the LOME may be artefacts of record failure.
A recent analysis of sequence stratigraphic architecture in Upper Ordovician sections in Gondwana and tropical Laurentia suggested that the first pulse of the LOME does not in fact coincide with cooling and sea level fall but rather with the first major interglacial episode of the Hirnantian stage [19]. This interpretation, if supported, would require a radical reconsideration of the nature of the LOME and the lessons that can be drawn from it regarding general relationships between climate change and extinction in the fossil record.
A comprehensive analysis of latest Ordovician extinction selectivity patterns and their context in broader Late Ordovician -Early Silurian selectivity patterns is therefore needed to evaluate relative support for different causal models of the LOME. Here we use a large and taxonomically standardized occurrence database [22 -24] to examine Late Ordovician -Early Silurian extinction patterns within one of the most diverse and well-preserved Palaeozoic clades, the rhynchonelliform brachiopods.

Material and methods
See the electronic supplementary material for additional details on data and analyses.

(a) Database
Our dataset was compiled from the literature and from ongoing research programmes in Durham and Copenhagen. Local stratigraphic ranges of rhynchonelliform brachiopods were recorded to the species level where possible. In the literature, though, only generic lists are often published. Thus, we omitted species-level data for this study and focused on genera. Stratigraphically, all data have been correlated with the British stage system using Bergströ m et al. [25] in order to standardize correlation. For this study, we have translated the regional stages into global units, again using Bergströ m et al. [25]. Taxonomic revisions dealing with specific genera have been assessed and their recommendations implemented wherever possible to avoid synonyms. Brachiopod species, and genera exhibit strong depth preferences and benthic assemblage (BA) zones have long been used as depth indicators in the early Palaeozoic [26 -28]. Each genus in our dataset was assigned to a range of BAs based on information in the literature, often with reference to associated fauna or lithology. In cases for which literature data were not available, assignments were based on the authors' experience with faunal and facies associations in the field. Further details regarding the database are available in previous publications [22 -24].
(b) Quantifying extinction and risk predictors All data manipulations and analyses were carried out in the R programming environment [29]. For each of 23 Late Ordovician -Early Silurian timeslices resolvable in our database (electronic supplementary material), we tallied all genera that were sampled in at least one locality (i.e. the timeslice falls within the local stratigraphic range of the genus at least one locality). To avoid edge effects, we excluded the Sandbian (early Late Ordovician) and Telychian (late Early Silurian) intervals. In the remaining 19 timeslices, the number of extant genera ranges from 121 (earliest Rhuddanian) to 238 (earliest Katian). Within each timeslice, we count genera as survivors if they are known to occur in younger timeslices and extinctions if they are not. We further divided survivors into genera that are sampled in at least one region in the immediately succeeding timeslice and 'Lazarus' genera [21,30] that disappear from the record for at least one timeslice but subsequently reappear. In keeping with common palaeobiological practice, we analysed geographical and stratigraphic ranges at the genus level, because species are inconsistently identified in the literature. If the determinants of extinction risk differ between the species and genus levels, then species-poor genera are expected to be at higher risk than relatively speciose genera. To account for such differences, we tallied the number of named species assigned to each genus during each timeslice (species richness). Genera lacking named species (e.g. all occurrences in that timeslice were recorded as 'sp.') were assigned a minimal value of 1.
Within each timeslice, we tabulated several aspects of geographical distribution for each genus. These included great circle distance (the maximum distance between any two localities occupied by a genus during a given timeslice; genera occurring at only a single locality were assigned a value of 1), grid cell occupancy (the number of 108 latitude by 308 longitude grid cells that contain at least one locality occupied by a genus during a given timeslice) and number of palaeocontinents (the number of palaeocontinents and terranes occupied by a genus during a given timeslice). We quantified latitudinal distribution by tabulating the minimum and maximum absolute palaeolatitude at which each genus occurs in each timeslice and its absolute palaeolatitudinal range (e.g. from 08 to 908). We quantified bathymetric distribution by tabulating the minimum and maximum depth (as measured by BA membership) and the depth range of each genus. To evaluate whether genera entirely or largely confined to cratonic seaways were harder hit than open-shelf genera we tabulated % cratonic localities (the percentage of all localities occupied by each genus in each timeslice that were located in cratonic seaways).
The sharp decline in the amount of preserved sedimentary rock between the Katian and the Hirnantian [13] raises the possibility that some of the apparent extinction during the first pulse of the LOME reflects record failure rather than reflecting true extinction. No existing database captures the global distribution of sedimentary rock with spatio-temporal resolution sufficient to compare directly with our database. Instead, in each timeslice, we categorized each locality in our database as continuous if at least one brachiopod genus occurred at that locality in the immediately succeeding timeslice, and discontinuous if no brachiopod genera are documented from that locality in the succeeding timeslice. A given locality may appear to be discontinuous in the immediately succeeding timeslice for several reasons, including (i) absence of sedimentary rock owing to non-deposition and/or subsequent erosion, (ii) lack of exposure, (iii) lack of expert collecting effort and (iv) lack of appropriate palaeoenvironments (e.g. rocks are terrestrial, or anoxic marine black shales), but all of these scenarios represent failures of preservation and/or collection that have the effect of truncating true stratigraphic ranges. To determine whether such truncations were an important predictor of apparent extinction risk, we rspb.royalsocietypublishing.org Proc. R. Soc. B 283: 20160007 tabulated % discontinuous localities (the percentage of all localities occupied by a genus during a given timeslice that have no records in the succeeding timeslice).
For the latest Katian timeslice (Katian 5), we undertook a further set of analyses to evaluate whether extinction risk was influenced by interactions between biogeographic distribution and shifting climate zones (see the electronic supplementary material). Using published climate simulations [31], palaeogeographic reconstructions [32,33] and observed palaeogeographic distributions, we estimated minimal thermal tolerance ranges for all genera extant during Katian 5. We then determined whether each regional population in a given genus would have been able to access habitat within this thermal range by shifting its range equatorward during the greenhouse -icehouse transition (electronic supplementary material, figure S1). Genera with at least one population that should have been able to remain within its thermal tolerance range without crossing an ocean basin were counted as predicted survivors. Those for which no suitable habit would have been available without dispersing across open oceans were counted as predicted extinctions. For comparison, we also made predictions based on the opposite scenario: an assumed icehouse state in the latest Katian transitioning to a Hirnantian greenhouse state.

(c) Modelling relationships between predictors and extinction risk
We examined relationships between the predictors described above and extinction patterns in the nine timeslices (out of 19 total) that had at least 10 observed genus extinctions. We used generalized boosted regression models (GBMs) [34,35] to determine which predictors were most strongly associated with extinction risk in each interval and the marginal effect (the effect with all other variables held constant) of varying each predictor on extinction risk. GBM models are useful for delineating extinction selectivity patterns throughout the study interval, but ensemble-based models of this kind cannot easily be compared in a maximum-likelihood framework. To make explicit models comparisons for the Katian 5 timeslice, we constructed multivariate binomial logistic regression models [36] using the subset of predictors that were included in 75% or more of the top-ranked models (those that had AICc scores within four points of the 'best' model). These predictors were palaeolatitude range, minimum depth, % cratonic and % discontinuous. We then evaluated relative support for a model including only these predictors and the intercept, these predictors and the intercept plus predicted extinctions and survivors in a greenhouse-icehouse transition, and these predictors and the intercept plus predicted extinctions and survivors in an icehouse-greenhouse transition. We used a classification tree to illustrate how the predictors in the best-supported model partition Katian 5 genera into extinctions and survivors.

Results and discussion
The performance of GBM models in accurately classifying genera as extinctions or survivors varied across intervals (figure 1), with the Katian 1, Katian 3 and Hirnantian models being notably poor. The marginal effects of predictors on extinction risk in each interval, and their relative influence on predicting extinction risk/survival, are shown as partial dependence plots in figure 1. Relative influence of predictors varied among intervals, but compared with other timeslices, the timeslice coinciding with the major pulse of the LOME (Katian 5) stands out for the unusual relative influence and marginal effects of palaeolatitude range, minimum depth, % discontinuous and % cratonic. Great circle distance is also an important predictor of extinction risk during Katian 5, but similar relationships between great circle distance and marginal risk are apparent in several other timeslices. In contrast, although palaeolatitudinal range is an important predictor of extinction risk in some other timeslices, no other timeslices show a similar relationship between palaeolatitudinal range and marginal risk: genera with extremely restricted palaeolatitudinal ranges (often those restricted to a single locality or region) exhibit elevated marginal risk in some timeslices, but only Katian 5 exhibits a broad plateau of elevated marginal risk among genera with palaeolatitudinal ranges less than about 358. Katian 5 brachiopod extinctions are also strikingly selective with respect to bathymetric distribution. Genera that ranged into shallower waters (BA 1 -2) experienced much lower extinction rates than those restricted to deeper waters (BA 3-6). The pattern is not limited to a single region but is apparent across a broad palaeolatitudinal range (figure 2), implying that it reflects the operation of a globalscale environmental driver and that the previously noted demise of the widespread deep water (BA 5-6) Foliomena fauna [14,22] is only the most extreme manifestation of a broader selective sweep. Interpretation of this pattern is complicated by uncertainty regarding the actual depth ranges represented by BAs (which are restricted to shelf and slope settings) [26], but the sign of the depth signal is informative. Cooling alone would be expected to steepen bathymetric temperature profiles and preferentially cause the extinction of genera restricted to warmer surface waters rather than those restricted to cooler deeper water. It is more likely, therefore, that the bathymetric extinction gradient reflects some other change in water mass characteristics.
Several lines of evidence, including biomarkers [37], nitrogen isotopes [38,39], molybdenum isotopes [14,40], iron speciation [14] and black shale distributions [39], suggest that the onset of the Hirnantian icehouse climate state was accompanied by increased oxygenation of shelf environments. Many of the species found in deeper-water environments during the late Katian were small and thin-shelled [41], like modern low-oxygen specialists [42]. Benthic species that had adapted to relatively low-oxygen conditions prevailing in deeper waters during the Katian greenhouse climate state may have been driven extinct by expansion of more oxygenated waters with their incumbent shallow-water faunas. This scenario is consistent with the contemporaneous extinction of multiple graptolite lineages thought to have inhabited denitrifying waters on the margins of oxygen minimum zones [16,39,43,44], including many older and previously extinction-resistant lineages [45]. An alternative interpretation is that the extinction of deeper-water genera reflects expansion of oxygen minimum zones and associated sulfidic conditions [14], but recent studies [46] do not support this hypothesis. Changes in oxygenation would not necessarily affect all taxa similarly, and whether such changes can explain other selective patterns such as the disproportionate extinction of trilobites with planktonic life stages [17] is not clear.
Both the first (end Katian) and second (late Hirnantian) pulses of the LOME exhibit slightly elevated extinction risk for genera occurring primarily in cratonic seaways relative to those occurring primarily in open-shelf settings (figure 1). Although this effect is weak, it is notable because it is opposite to that observed in most other intervals, during which genera restricted to open-shelf environments exhibit slightly greater rspb.royalsocietypublishing.org Proc. R. Soc. B 283: 20160007 risk than cratonic genera. Preferential extinction of cratonic genera contrasts with preferential survival of such taxa during Meso-Caenozoic mass extinctions [47], but is consistent with previous suggestions that some latest Katian extinctions were driven by eustatically forced regression and draining of cratonic seaways [7,48]. It should be noted, however, that % discontinuous is also weakly predictive of extinction risk in both Katian 5 and (to a lesser extent) the Hirnantian. Moreover, in Katian 5, the marginal risks associated with % discontinuous and % cratonic are rather similar (figure 1).
The positive association between % discontinuous and extinction risk, though weak, raises the possibility that some proportion of the extinctions that appear to occur at the Katian-Hirnantian boundary represent artificial range truncations resulting from reduced preservation probability. Such truncations would not be expected to affect all genera equally. The draining of cratonic seaways may have led to genuine extinctions of genera with dominantly cratonic distributions, but these genera might also be expected to have experienced disproportionate decline in preservation potential. In shelf areas most preserved sequences record substantial shoaling [3,6,19,49,50], raising the possibility that deeper-water genera may also have experienced a disproportionate decline in preservation potential owing to reduced sampling of appropriate environments.
We evaluated the likelihood that observed extinction selectivity patterns are driven by stratigraphic range truncations by comparing mean predictor values in genera that appear to go extinct in Katian 5, 'range-through' genera that survive but have a post-Katian 5 sampling gap, and genera that survive and occur in the subsequent (early Hirnantian) timeslice. In the extreme case that all apparent extinctions actually represented stratigraphic range truncations, the factors  The unusual importance of palaeolatitudinal range in the latest Katian implies that changing sea surface temperatures played a major role in driving extinctions. Latitudinal ranges are closely linked to thermal tolerance ranges among extant aquatic ectotherms [51][52][53], and sea surface temperatures are the single most important determinant of biogeographic structure in the coastal oceans [54,55]. The thermal tolerance ranges of some extant marine invertebrate species have been conserved on million-year time scales [56], and many species shifted their ranges in response to Caenozoic and Quaternary climate changes [55,57,58].
Preferential survival of genera with wide palaeolatitudinal ranges does not necessarily implicate cooling; broad thermal tolerance range would be expected to buffer against extinction during any rapid climate change, and the Latest Ordovician world likely experienced glacial-interglacial oscillations analogous to those of the Caenozoic. It has been argued that the first pulse of the LOME corresponds not to the greenhouse-icehouse transition but rather to the first interglacial within the icehouse state [19], and the selective extinction of narrowranging genera is also consistent with this hypothesis. However, a multivariate logistic regression model that includes the four selected predictors plus predictions about which genera would be expected to go extinct owing to habitat loss during a greenhouse-icehouse transition strongly outperforms both the model including only the four original predictors and the model including predictions about which Genera were assumed to be present throughout their recorded palaeolatitudinal and depth ranges. Depth ranges are based on brachiopod depth assemblage zones ranging from 1 (shallowest) to 6 (deepest). Only palaeolatitude depth cells occupied by at least one genus are outlined; those occupied by fewer than three genera are indicated by grey fill. The first (Katian 5) pulse of the Late Ordovician mass extinction is distinguished by the overall intensity of extinction, but also by the strongly selective extinction of deeper-water genera.
rspb.royalsocietypublishing.org Proc. R. Soc. B 283: 20160007 genera would be expected to go extinct in a hypothetical icehouse-greenhouse transition (table 1). Consequently, we conclude that cooling and equatorward shift of climate zones led to the extinction of many species (and thus genera) that had narrow thermal tolerance ranges and inhabited coastlines or cratonic seaways of limited latitudinal extent. Extinction selectivity patterns during Katian 5 can be visualized with a classification tree ( figure 3). Unsurprisingly, given that our predictions for extinction versus survival in a greenhouse-icehouse transition are necessarily crude estimates, palaeolatitudinal range remains an important predictor of extinction risk even in models that include these predictions. Only 17% of genera with palaeolatitudinal ranges more than 358 go extinct (figure 3), but of those with ranges less than 358, extinction is substantially higher in genera confined to depths of BA3 or deeper than in genera that ranged into shallower waters. Predicted extinction in a greenhouse-icehouse transition is selected as a further split in both the exclusively deep water and shallower subsets of the narrow-ranging genera, with further splits on % cratonic and % discontinuous among narrow-ranging deep-water genera that are not otherwise predicted to go extinct ( figure 3). This analysis thus confirms that aspects of thermal tolerance range and bathymetric distribution seem to have been more important than occupying cratonic seaways and/or regions with widespread stratigraphic truncations in determining which genera went extinct and which survived.   Table 1. Comparison of three different multiple logistic regression models using the modified Akaike information criterion (AICc): 'no prediction' (extinction palaeolatitudinal range þ minimum depth þ % cratonic þ % discontinuous), 'greenhouse-icehouse prediction' (extinction palaeolatitudinal range þ minimum depth þ % cratonic þ % discontinuous þ predicted extinction/survival in a greenhouse-icehouse transition) and 'icehouse-greenhouse prediction' (extinction palaeolatitudinal range þ minimum depth þ % cratonic þ % discontinuous þ predicted extinction/survival in an icehouse-greenhouse transition). The 'greenhouse-icehouse prediction' model, in which the expected survival or extinction of a genus in an icehouse climate state is predicted based on its distribution in an assumed greenhouse climate state, is favoured over the other models by evidence ratios The poor performance of the Hirnantian GBM model (figure 1) indicates that none of the variables considered here are particularly effective predictors of extinction risk during the smaller second pulse of the LOME in the mid-to late Hirnantian. A previous analysis of the Laurentian record also failed to find effective predictors of extinction risk for the second LOME pulse [13]. It has been suggested that extinctions during this interval are related to warming and shoaling of oxygen-poor waters during partial deglaciation [14,39], but we observe no strong depth or palaeolatitudinal signal in the distribution of extinction risk. Future analyses incorporating additional environmental factors (e.g. the percentage of sites occupied by a given genus in the late Hirnantian that are overlain by laminated black shales in the Rhuddanian) may be able to resolve a clearer selective signal.
Our results suggest an explanation for one of the enduring questions about the LOME: why did the Late Ordovician greenhouse -icehouse transition cause a major mass extinction when subsequent greenhouse -icehouse transitions did not? In fact, both the late Eocene-Oligocene transition and the Plio-Pleistocene descent into full icehouse conditions were associated with significant extinction pulses [59][60][61][62][63], but these pulses were predominantly focused on tropical and subtropical regions and were not of comparable magnitude. The contrast may reflect three important differences in starting state between Late Ordovician and Caenozoic greenhouse -icehouse transitions.
(1) Early Palaeozoic dissolved oxygen levels may have been substantially below typical modern levels, especially in outer shelf and slope settings [64]. Whereas the Late Ordovician greenhouse -icehouse transition seems to have been accompanied by a substantial increase in shelf oxygenation [37 -39,43], there is little evidence for a comparable increase during the Caenozoic transition. There are many potential explanations for this contrast, but one key factor may be Phanerozoic evolutionary trends that have acted to deepen the e-folding depth of organic remineralization [65]. Such a shift would deepen oxygen minimum zones and reduce the exposure of shallow benthic habitats to substantial fluctuations in dissolved oxygen content.
(2) The late Katian world was an extreme highstand with extensive flooding of continents, but cratonic seaways were of limited extent in the Caenozoic and almost entirely gone by the Pliocene [66] (electronic supplementary material, figure S3). Loss of cratonic faunas, many of which probably experienced substantial oxygen fluctuations [67] and would have had difficulty dispersing into open-shelf settings [68], was almost certainly a more important driver of extinctions during Late Ordovician cooling and sea level fall than during Caenozoic cooling episodes.
(3) Except for the supercontinent of Gondwana, the Late Ordovician world was characterized by widely dispersed island continents and an abundance of isolated terranes [24,33], whereas the Caenozoic world was generally characterized by extensive continuous north -south coastlines [69] (electronic supplementary material, figure S3). Caenozoic continental configurations would thus have allowed many marine species to shift their ranges equatorward in response to shifting climate zones without having to disperse across open ocean basins. The Late Ordovician configuration may have had a greater tendency to 'trap' species on coastlines from which they could not easily shift their ranges equatorward, accounting for the selective extinction of genera with narrow palaeolatitudinal ranges.
Our finding may also help to explain the observation that the LOME was taxonomically relatively non-selective and seems to have had far less long-term impact on the taxonomic and ecological structure of marine ecosystems than other major mass extinction events [70,71]. Whereas both the Permian-Triassic and Cretaceous-Palaeogene events removed groups that had previously been diverse and ecologically dominant, very few higher taxa disappeared during the LOME. Our analyses show that among rhynchonelliform brachiopods extinction was in fact strongly selective, but along bathymetric and biogeographic gradients that would have affected most major taxa. Consequently, although taxonomic losses were very high, these losses were relatively evenly distributed across ecospace. With the re-establishment of background environmental conditions surviving lineages rapidly diversified [72] into similar regions of ecospace, giving rise to ecosystems structurally similar to those that preceded the LOME.
Data accessibility. Data used in the analyses reported here are deposited at Dryad: http://dx.doi.org/10.5061/dryad.66vb3. The data can also be downloaded, along with R scripts to reproduce all analyses, at GitHub: https://github.com/sethfinnegan/Ord.Brach.Extinctions. Code.git. See the electronic supplementary material for further instructions.