Long duration of benthic ecological recovery from the early Toarcian (Early Jurassic) mass extinction event in the Cleveland Basin, UK

Recovery intervals are crucial periods of biotic change that allow new groups to radiate into the emptied ecospace after extinction events (Erwin et al. 1987; Kauffman and Erwin 1995; Harries et al. 1996; Kauffman and Harries 1996; Solé et al. 2002; Twitchett 2006). However, the definition of this post-extinction period is highly dependent on the scale and groups of organisms studied (Danise et al. 2015; Dai et al. 2018; Song et al. 2018; Atkinson and Wignall 2019). Different groups display different post-extinction diversity dynamics due to their intrinsic evolutionary rates, which means that recovery of the whole ecosystem will take longer than the recovery of one clade alone (Stanley 2009; Guex et al. 2012; Song et al. 2016; Dai et al. 2018). The complexity of ecosystems is also important because local environments may be more susceptible to environmental change and may perpetuate adverse conditions (Sheehan 1985; Stanley 1988; Atkinson and Wignall 2019; Müller et al. 2020).

It is important to define what is meant by ecological recovery and several models have been developed to do so. Such models use a variety of metrices to document the overall ‘healthiness’ of a community, typically including a measure of the number of species present in different intervals. This can be measured through horizon-level or time-binned species counts (e.g. Hallam 1987; Mander et al. 2008; Pugh et al. 2014). This is a distinctly different measure of diversity from another commonly used metric – evenness, or conversely, dominance – which measures the ratio between each species to determine whether an assemblage is numerically dominated by a small number of species (Schubert and Bottjer 1995; Twitchett 2006; Danise et al. 2013). A third common component of recovery models is to include a measure of ecological diversity. This may be the number of different life modes present and the number of species or specimens that occupy it.

Twitchett (2006) found that the levels inhabited by organisms relative to the seafloor (ecological tiering) was a good indicator for recovery (Fig. 1). This includes burrowers at different depths (deep or shallow infaunal), organisms that sit low on the surface of the sediment (low-level epifaunal) and organisms that reach up into the water column from the seafloor, such as corals or crinoids (erect epifaunal). Following a mass extinction event, ecological tiering collapses, leaving communities dominated by low-level epifaunal organisms (Twitchett 1999). Another method often used in recovery studies is a measure of body size (Urbanek 1993; Payne 2005; Atkinson and Wignall 2020). In the aftermath of mass extinctions, organisms tend to be of diminutive stature due to either high rates of juvenile mortality or small adult body sizes; the return to large body sizes can be seen as a sign that the ecosystem has recovered (e.g. Atkinson and Wignall 2020).

Kauffman and Erwin (1995) divided the post-extinction interval into two distinct phases: the survival interval and the recovery interval (Fig. 1). The former is characterized by low taxonomic richness with fossil communities dominated, in terms of specimens, by a handful of small-sized species that occupy low levels of ecological tiering (e.g. low-level epifaunal). Such species are usually inferred to have been disaster taxa and opportunists that were likely R strategists – attaining sexual maturity rapidly and giving rise to large numbers of offspring (Guex 1992; Harries et al. 1996). The succeeding recovery interval of Kauffman and Erwin (1995) is characterized by a period of increasing taxonomic richness, with new species arising from surviving lineages and species migrating into the area from theorized refugia.

The model for post-extinction recovery of Twitchett (2006), which is based on the end-Permian mass extinction event, incorporates key ichnotaxa and has additional ecological parameters in four post-extinction phases (Fig. 1). The first equates to the survival interval of Kauffman and Erwin (1995), but includes a depauperate trace fossil assemblage. It is less clear how the second stage of recovery in the Twitchett (2006) model may relate to the Kauffman and Erwin (1995) model; it may be equivalent to the latest survival interval or the earliest recovery interval (Fig. 1). Subsequent phases are equivalent to the recovery interval of Kauffman and Erwin (1995), but Twitchett (2006) dilutes the importance of species richness alone by including the community structure, body size and morphological complexity (Fig. 1). This later model includes the expansion of seafloor tiering and the repopulation of ecological niches. Ecosystems in the aftermath of a mass extinction event have been likened to a ship sailed by a skeleton crew (Foster and Twitchett 2014), in that all ecological roles are present, but with each life mode represented by only one or two species. Consequently, species diversity remains low, while ecological diversity appears to be unaffected.

Post-extinction recovery has been best studied for the Permian–Triassic, end-Triassic and Cretaceous–Paleogene extinction events (Hallam 1987; Hansen 1988; Nützel 2005; Barras and Twitchett 2007; Hautmann et al. 2008; Sessa et al. 2009; Hu et al. 2011; Ros et al. 2011; Chen and Benton 2012; Aberhan and Kiessling 2014; Pugh et al. 2014; Damborenea et al. 2017; Foster and Sebe 2017; Dunhill et al. 2018; Song et al. 2018; Alvarez et al. 2019; Atkinson and Wignall 2019, 2020; Whittle et al. 2019). By comparison, recovery following second-order mass extinction events, such as that of the early Toarcian (Early Jurassic), is less well understood. The early Toarcian event was associated with the onset of a globally occurring negative carbon isotope excursion, produced by the release of isotopically light CO₂ by the Karroo–Ferrar large igneous province (Palfy and Smith, 2000) and/or methane hydrates (Hesselbo et al. 2000; Kemp et al. 2005), warming of 5–10°C (Bailey et al. 2003; Ruebsam et al. 2019), an enhanced hydrological cycle (Bailey et al. 2003), marine transgression (Hallam and Wignall 1999) and the Toarcian oceanic anoxic event (TOAE; Jenkyns 1988), which was largely a regional phenomenon.

The TOAE was particularly well developed in the epicontinental basins of NW Europe, with the deposition in the Serpentinum Zone of finely laminated, organic-rich (black) shales, such as the Schistes Cartons of France, the Posidonia Shale of Germany and the Mulgrave Shale (formerly the Jet Rock) of northern England (Jenkyns 1988). Within these restricted basins of NW Europe, the TOAE is thought to have been the main cause of the early Toarcian mass extinction event, especially for benthic faunas, because the extinction level is coincident with the beginning of black shale deposition at the end of the Tenuicostatum Zone (Little 1995; Little and Benton 1995; Schmid-Röhl et al. 2002; Wignall et al. 2005; Caswell et al. 2009; Danise et al. 2013, 2015; Caswell and Frid 2017). In other regions, such as the Sub-Boreal and Tethyan realms, these oxygen-restricted facies were weakly or not at all developed during the Toarcian and, instead, a rapid increase in temperature may have been the main driver of extinction (e.g. McArthur et al. 2008; Gómez and Goy 2011; Baeza-Carratalá et al. 2015; Miguez-Salas et al. 2017, Danise et al. 2019). Earlier extinctions occurred at the Pliensbachian–Toarcian boundary in ammonites (Dera et al. 2010; Caruthers et al. 2013) and in Tethyan coral faunas (Vasseur et al. 2021), the cause(s) of which may have been warming from the first major eruptive phase of the Karroo–Ferrar large igneous province (Caruthers et al. 2013; Percival et al. 2015).

Early studies on extinction and recovery spanning the Pliensbachian and Toarcian stages were conducted at low resolution, with the species counts compiled by ammonite zone across much of Europe (Hallam 1976, 1986, 1987, 1996). These studies showed a delayed recovery among the bivalves (with species richness beginning to approach pre-extinction values by the latest Toarcian), whereas the brachiopods displayed little to no recovery within the studied time interval. Similar recovery patterns were found in more recent regional studies from the NW Caucasus (Ruban 2004), the Boreal–Arctic region (Zakharov et al. 2006) and in the Andean Basin of Chile (Aberhan and Fürsich 2000), with the latest Toarcian diversity not regaining its pre-extinction richness. However, the ecological structure of mid Toarcian bivalve assemblages within the Andean Basin attained a level of ecological complexity similar to that of the late Pliensbachian, despite a lower species richness (Aberhan and Fürsich 2000). Recovery occurred at a greater rate within the shallow water succession of northern and central Spain, where brachiopods recovered their diversity and body size by the Variabilis Zone (García Joral and Goy 2009; García Joral et al. 2011, 2018) or, locally, earlier still during the Serpentinum Zone (Danise et al. 2019). This ‘enhanced’ recovery was likely linked to the absence of marine oxygen depletion (García Joral and Goy 2009; García Joral et al. 2011; Danise et al. 2019).

The Cleveland Basin (northern UK) was home to high-diversity benthic fossil assemblages in the late Pliensbachian, with most species being evenly represented (no one species dominated) and featuring a wide range of ecological groups occupying all levels of benthic tiering (Little 1995; Caswell et al. 2009; Danise et al. 2013; Caswell and Frid 2017). Although there was an extinction event in ammonites at the Pliensbachian–Toarcian boundary (the loss of the family Amaltheidae, which defines the end of the Pliensbachian; Howarth 1955), most of the benthic taxa passed through into the early Toarcian unaffected.

The benthic extinction horizon occurred within a few metres of the sedimentary section in the Semicelatum Subzone of the Tenuicostatum Zone at c. 183 Ma (Fig. 2), leading up to the first appearance of thick laminated black shales (Little and Benton 1995; Caswell et al. 2009; Danise et al. 2013; Caswell and Frid 2017). With the black shales of the succeeding Serpentinum Zone, there was a marked shift in the structure of the benthic assemblages, with a collapse of the benthic tiering structure (only low-level epifauna remained) and assemblages dominated by large numbers of the opportunistic, and probable low-oxygen specialist, bivalves Pseudomytiloides and Bositra (Little 1995; Danise et al. 2013; Caswell and Dawn 2019). At the same time, sphaeromorph- and Tasmanites-dominated plankton communities replaced earlier dinoflagellate communities (Slater et al. 2019; Danise et al. 2022).

The benthic faunal diversity remained low well into the Commune Subzone of the Bifrons Zone, c. >2 myr after the extinction horizon (Fig. 2), despite the negative carbon isotope excursion lasting only c. 0.9 myr (Suan et al. 2008). The Commune Subzone records the onset of the recovery, as documented by the first appearance of shallow and deep infaunal tiering. The return of infaunal species marked the loss of well-laminated mudstone and is represented by a bloom of the deposit-feeding bivalve Dacryomya ovum (Little 1995; Caswell et al. 2009; Danise et al. 2013; Caswell and Dawn 2019). Other new species arose from surviving lineages, such as the deep infaunal suspension-feeding bivalve Gresslya donaciformis (Little 1995; Caswell et al. 2009; Danise et al. 2013). Epifaunal tiering increased shortly afterwards, with the first post-extinction occurrence of crinoids (Little 1995). Although plankton communities returned to dinoflagellate-dominated assemblages within the Bifrons Zone (Slater et al. 2019; Danise et al. 2022), neither trace nor body fossil diversity (both taxonomic and functional) had returned to pre-extinction levels at that time (Little 1995; Little and Benton 1995; Martin 2001, 2004; Danise et al. 2013, 2022; Caswell and Frid 2017).

All previous studies of extinction and recovery during the Toarcian in the Cleveland Basin have their upper limits near the top of the Bifrons Zone because the succeeding ammonite zones of the Middle and Upper Toarcian were removed by later Mid-Jurassic erosion at most of the Yorkshire sections (Knox 1984; Milsom and Rawson 1989; Little 1995; Caswell et al. 2009; Powell 2010; Danise et al. 2013), meaning that our knowledge of the recovery pattern from the early Toarcian extinction event is incomplete. However, there is one section on the North Yorkshire coast, at Ravenscar, where a near-complete middle and upper Toarcian section crops out in the cliffs and foreshore (Dean 1954; Knox 1984; Powell 1984). A few studies have sampled this section for palynology (Riding 1984), trace fossil assemblages (Martin 2001, 2004) and geochemical data (Newton et al. 2011), but all at relatively low resolution. Here, we report a high-resolution sampling study of benthic macrofossils through the middle and upper Toarcian from the Ravenscar section using similar methodologies to the previous studies of Little (1995) and Danise et al. (2013). We incorporate new and existing data to produce new species ranges and palaeoecological data for the entire late Pliensbachian to Toarcian time interval in the Cleveland Basin.

During the Jurassic, the Cleveland Basin consisted of a series of small extensional basins within a shallow epicontinental seaway (Powell 2010). The oldest unit within this study area is the Staithes Sandstone Formation of the Pliensbachian Stage (Powell 1984). This fine- to medium-grained sandstone is interpreted as storm wave-influenced, sublittoral and shoreface sands (Howard 1984; Powell 2010). This is followed by the Cleveland Ironstone Formation, which is divided into the Penny Nab and Kettleness members (Powell 1984; Howard 1985). This formation consists of cyclical rhythms of chamosite–ooidal ironstone seams, mudstones and sandstones, which may have been influenced by glacio-eustasy (Howarth 1955; Howard 1985; Korte and Hesselbo 2011). These sediments are thought to have been deposited within a sediment-starved, laterally extensive shallow lagoon or littoral sea (Powell 2010).

Following the Cleveland Ironstone, a marine transgression resulted in the deposition of the Whitby Mudstone Formation, around the Pliensbachian–Toarcian boundary (Powell 1984). This argillaceous unit is represented in most coastal sections by three members: the Grey Shale; the Mulgrave Shale (formerly the Jet Rock Member of Powell 1984); and the Alum Shale (Fig. 2). The Mulgrave Shale Member has the greatest total organic carbon content of the sequence and represents the peak transgression (Hallam 1997; McArthur et al. 2008; Pearce et al. 2008). The water depth has been estimated at 50–100 m, although this was likely to have been towards the shallower estimate (Hallam 1997; Wignall et al. 2005; Powell 2010).

Marine regression began during deposition of the Alum Shale Member (Powell 2010). Across much of the Cleveland Basin, the Alum Shale Member is truncated in the Crassum Subzone of the Bifrons Zone by a period of Mid-Jurassic erosion. However, a more complete section is preserved within the Peak Trough Member (Fig. 3), which was actively subsiding during the Jurassic (Milsom and Rawson 1989), and this upper Toarcian section is exposed in the cliffs and wave-cut platforms near Ravenscar, 6 km south of Robin Hood's Bay (Fig. 3). This succession includes the uppermost members of the Whitby Mudstone Formation: the topmost Alum Shale, the Peak Mudstone and the Fox Cliff Siltstone members (Fig. 2). As the names of these members indicate, they exhibit an increasing grain size and are a continuation of a regressive sequence (Knox 1984).

In the Ravenscar section, the Whitby Mudstone Formation grades upwards into the Blea Wyke Sandstone Formation (Knox 1984). This formation is separated into the Grey and the Yellow Sandstone members (Knox 1984), the contact between the two forming a sharp colour change in the cliff running from Ravenscar to Blea Wyke. The Grey Sandstone Member has a higher mud content and features hard siderite-rich horizons, whereas the Yellow Sandstone Member is a cleaner sandstone. Both members are intensely bioturbated and their deposition occurred within the sublittoral zone, above the wave base (Knox 1984).

During the summers of 2013 and 2017, we collected from 37 sample horizons, spaced mostly at 1 m intervals (range 0.5–4 m due to differences in rock availability), from 45 m of the middle and upper Toarcian stratigraphy, from the topmost Alum Shale Member (Bifrons Zone) to the top of the Yellow Sandstone Member (Aalensis Zone) in the cliff and foreshore exposures between Ravenscar and Blea Wyke Point, North Yorkshire (54° 23′ 58″ N, 00° 28′ 33′ W′) (Fig. 3). The sedimentary logs of Howarth (1962), Knox (1984) and Dean (1954) were used to locate the sample points in the stratigraphic section. At each sample horizon, c. 1 kg of bulk rock was collected and the benthic species richness and abundances were counted from wave-cut platforms (where available) or cliff exposure, for up to 2 h per horizon.

The bulk samples were taken back to the University of Leeds and broken using a hammer and chisel to expose their benthic fossil content. The bulk samples gave species diversity and abundance data from equal rock volumes. Counts were made in the field from wave-cut platforms and cliff sections to increase the size of the available dataset for a total evidence approach. These datasets were plotted separately as a test of sampling bias due to differences in exposure quality, but were combined for the range chart and palaeoecological analysis.

For brachiopods and bivalves, the number of individuals was determined by counting the sum of articulated and disarticulated valves. For gastropods and scaphopods, single shells were counted as one individual. For crustaceans, cephalothorax specimens were counted as one individual and, for crinoids, discrete clusters of ossicles were considered as one individual. Fossils in the Yellow Sandstone Member are largely preserved as moulds, so the surfaces of the bulk-rock samples from this unit were cast using silicone rubber and these casts were used to identify the fossils. Benthic tiering and feeding strategies were assigned using the scheme from Bambach et al. (2007; see Supplementary Material). Trace fossil occurrences for the Ravenscar section were taken from Martin (2001).

Benthic species richness and abundance counts from the Ravenscar section were combined with those of Danise et al. (2013) to form a dataset spanning the late Pliensbachian Staithes Sandstone Formation to the Blea Wyke Sandstone at the top of the Toarcian. It is important to note that Danise et al. (2013) used a larger rock volume (3 kg) for their bulk samples. This combined dataset was used to construct species range charts for that time interval, with range extensions from younger and older units in the basin taken from previously published work (Sowerby 1821; Morris and Lycett 1854; Tate and Blake 1876; Dean 1954; Howarth 1962, 1992; Hallam 1967; Howard 1984, 1985; Johnson 1984; Simms 1989; Ager 1990; Doyle 1990; Caswell et al. 2009; Atkinson and Wignall 2020). Combined species richness data (bulk + foreshore) were subsampled to standardize species richness on the basis of the number of fossil occurrences (Dunhill et al. 2018). Subsampling quotas were set at 50 occurrences per sampled horizon and any of these that did not meet this quota was omitted. In each case, subsampling was carried out across 1000 iterations. Simpson's index of diversity 1-D (SID) was used to assess the evenness of benthic assemblages per sample horizon for the full dataset, following the methods used in Danise et al. (2013). This returns a result between 0 and 1, where 0 is an assemblage dominated by a single species and 1 means all species are equally represented.

To study the ecological composition from sample horizons, we constructed percentage area plots for the 11 identified benthic palaeoecologies (erect epifaunal suspension feeders, low-level epifaunal suspension feeders, epifaunal grazers, epifaunal predators, semi-infaunal suspension feeders, semi-infaunal predators, shallow infaunal suspension feeders, shallow infaunal deposit feeders, shallow infaunal/epifaunal microcarnivores, shallow infaunal species with chemosymbionts and deep infaunal suspension feeders). Proportional diversity was assessed within the orders of bivalves, using raw species counts with taxonomic information from Cox et al. (1969).

An absolute timescale is provided from Gradstein et al. (2020) using dated ammonite zone boundaries (Fig. 2). For the upper Toarcian, the ammonite zone schemes of Knox (1984) and Dean (1954) have been correlated to the current scheme presented in Page (2003), which was used by Gradstein et al. (2020). As many old subzones have since been elevated to zones, there is a good correspondence for most of the zonal biostratigraphy of the Ravenscar succession. However, determining the position of some subzones is more problematic and it is necessary to detail a few points where the correlation between the schemes was not straightforward.

There is no evidence for the Fluitans Subzone on the Yorkshire Coast, so it seems likely that this was removed by Mid-Jurassic erosion below the ‘Terebratula beds’ that mark the base of the Dogger Formation (Dean 1954). The Mactra Subzone is indicated by the presence of a single Dumortieria moorei specimen found by Dean (1954) in the upper half of the Yellow Sandstone Member. The Pseudoradiosa Subzone is unconfirmed, but there is no sedimentary evidence for a break in deposition to suggest that it is not present. It is likely, instead, that the poor state of preservation of ammonites from this part of the section is the reason for this subzone not being recorded. Indeed, we collected several ammonites belonging to the genus Hudlestonia from the upper half of the Yellow Sandstone Member, but none that could be identified to species level and therefore used for biostratigraphy. The Dispansum Subzone of Knox (1984) and Dean (1954) corresponds well with the Dispansum Zone of Page (2003) and is confirmed by the presence of Phlyseogrammoceras dispansum, the index taxa for the Dispansum Zone, from much of the Grey Sandstone Member (Dean 1954). Despite this, the Gruneri and Insigne subzones of the Dispansum Zone were not confirmed. There is no record of the index ammonites from the Vitiosa Subzone, although this may be due to poor accessibility towards the top of the Peak Mudstone Member. The index ammonites for the Illustris Subzone were recorded in the Peak Mudstone by Dean (1954).

In total, 23 345 benthic macrofossil specimens were counted from 37 sample horizons in the Ravenscar section. When combined with data from Danise et al. (2013), this gave a total of 35 615 benthic fossils from 219 sample horizons (Figs 4–6). On the basis of the range charts and palaeoecological analysis, we distinguish four phases related to benthic extinction and recovery in the Cleveland Basin: (1) the pre-extinction interval (Margaritatus Zone to Tenuicostatum Zone); (2) the immediate post-extinction interval (Tenuicostatum Zone, Semicelatum Subzone to Bifrons Zone, Commune Subzone) represented by the communities tolerant to low oxygen levels; (3) the early phase recovery interval (Bifrons Zone, Commune Subzone to Dispansum Zone) defined by the reappearance of ecologies that were exclusively lost from the basin at the extinction horizon; and (4) the late phase recovery interval (Dispansum Zone to Aalensis Zone), when the species richness increased markedly and populated ecological niches.

Prior to the extinction event, the species richness was highly variable, reaching highs of 20 species per sample horizon (Fig. 5a). On average, these sample horizons record six benthic ecological life modes, with low-level epifaunal suspension feeders being the most diverse (up to ten species from one sample point), but also a rich infaunal diversity with up to three deep infaunal species and four shallow infaunal species from a single sample horizon (Fig. 7). Towards the extinction event, there was a gradual decrease in the abundance of deep infaunal suspension feeders (Danise et al. 2013). Immediately prior to the extinction event, shallow infaunal deposit-feeding bivalves, such as Palaeonucula navis, form a significant proportion of these assemblages and, for a short interval straddling the extinction event, the bivalve Nucinella sp. is briefly abundant, forming up to 79% of the benthos in specific horizons (Danise et al 2013; Caswell and Frid 2017).

Following the extinction event, species richness sharply decreased to typically two species and remained low throughout much of this time. During the deposition of laminated black shales, a distinctive series of fossil assemblages were present, formed of low-level epifaunal suspension feeders with an abundance of individuals. Monospecific shell pavements of the ‘paper pecten’ Bositra radiata occur near the top of the Grey Shale Member in the Semicelatum Subzone (Little 1995; Caswell and Coe 2013; Danise et al. 2015). These bivalves have valves up to 50 mm in height (Little 1995). B. radiata was found in only a narrow stratigraphic interval in the Grey Shale Member, although Caswell et al. (2009) report a transient reoccurrence of the species higher in the Mulgrave Shale, from the upper Exaratum Subzone, Serpentinum Zone.

Above this level, in the Falciferum Subzone (Serpentinum Zone), Bositra buchii occurs in thin, monospecific shell beds. This species is smaller than B. radiata, with a maximum shell height of c. 10 mm. Pseudomytiloides dubius occurs sporadically in low numbers throughout the Grey Shale Member, but dominates the benthic fauna in the organic-rich facies of the Mulgrave Shale, forming thin shell beds restricted to discrete horizons. These shell beds occasionally also contain a few individuals of Meleagrinella substriata and B. buchii (Caswell and Coe 2013; Danise et al. 2015). These assemblages gradually disappear in the lower Alum Shale Formation, with P. dubius last recorded from the Commune Subzone (Bifrons Zone).

Following the extinction event, there were no infaunal shelly components for slightly over 2 myr, representing the Sepentinum Zone to the lower part of the Bifrons Zone (Commune Subzone). Levels of ecological tiering increased during the Bifrons Zone, with the first significant infaunal component being the deposit-feeding bivalve  D. ovum (Little 1995; Danise et al. 2013; Caswell and Dawn 2019). At this time, small, simple, pyritized feeding traces appeared, indicating an increase in bioturbation in the sediment, and small (1 mm) Chondrites burrows were present towards the top of the Alum Shale Member (Martin 2004). Deep infaunal organisms returned shortly after the first appearance of D. ovum and were accompanied by a transient reoccurrence of erect epifaunal animals (crinoids). Their presence shows that a tiering structure was briefly restored (albeit following the skeleton crew model) during the Fibulatum Subzone (Bifrons Zone), c. 2.7 myr following the extinction event. Until now, this has been the extent of the recovery studied within the Cleveland Basin; above this stratigraphic level, our new results show a continuation of the increasing benthic species richness in both the subsampled data and in the data from bulk collections (Fig. 5).

For the top of the Alum Shale Member, the fluctuation of SID values results from transient, high abundances of the bivalve M. substriata (Fig. 6a). The ornate and thick-shelled trigoniid bivalve Vaugonia literata is also a common faunal element of the upper Alum Shale Member. The topmost Alum Shale Member marks a temporary highpoint in species richness before a decrease during the subsequent Peak Mudstone and Fox Cliff Siltstone members (Fig. 5a). In much of the Peak Mudstone and Fox Cliff Siltstone members, insufficient specimens were found in both the foreshore and bulk samples to allow for subsampling, which resulted in highly sparse SID values (Fig. 6a, b).

For much of the Peak Mudstone Member (notable for its dark coloration, small rounded phosphatic nodules and pyritized fossil wood), the number of ecologies present and the levels of tiering show a reduction from the underlying Alum Shale, with a temporary loss of both shallow and deep infauna (Fig. 7a). This loss of diversity is recorded in the bulk samples, suggesting that the decline in the Peak Mudstone Member is a true signal and not just an artefact of the generally poor foreshore exposure for most of this unit. This latter fact prevented Martin (2004) from reporting any trace fossils from the Peak Mudstone Member, although small pyritized Planolites were encountered during our field collections. Fossil abundances rose through the Fox Cliff Siltstone Member. This was mirrored in the increasing SID values towards a more evenly represented assemblage and a tiering structure similar to that in the upper Alum Shale, with the return of erect epifauna (Chariocrinus wuerttembergicus) and both deep and shallow infaunal elements (Fig. 7).

Although species richness increases through the Fox Cliff Siltstone Member, as recorded in both the foreshore and bulk collections, the ecological groups are still only represented by a few species each, resulting in a fluctuating ecological diversity (Fig. 7c). For example, shallow infaunal suspension feeders from the Fox Cliff Siltstone Member were represented solely by the brachiopod Lingula beani, with the previously abundant V. literata and D. ovum of the upper Alum Shale Member having disappeared from the basin (Fig. 4). Trace fossil assemblages are represented by Planolites, small Skolithos-type burrows, Chondrites and Rhizocorallium (Martin 2001).

This phase saw a marked increase in the number of species present, which resulted in the filling of ecological guilds, a cessation in the fluctuations of ecological diversity that marked the early recovery phases, and the occurrence of a broadly persistent trend towards increased diversity (Fig. 7c). The late phase recovery interval is also coincident with the return of sandy facies within the Cleveland Basin.

Fossil abundances continued to increase, reaching their acme for this section at the top of the Grey Sandstone Member (14 512 counted individuals). This is due to the extreme numbers of the tube worm Serpula deplexa, which formed ‘nests’. The abundance of this one species resulted in a sharp and dramatic reduction in the SID value at this level (Fig. 6). The species richness fluctuates throughout the Grey Sandstone Member, although with a greater number of species than the underlying two members. This trend is shown most obviously in the foreshore counts, whereas the bulk sampling shows a progressive increase in the number of species (Fig. 5c). Subsampling shows the converse, suggesting there was a progressive decrease in the number of species. This is probably due to the extremely high numbers of S. deplexa swamping the diversity (Fig. 5b). The number of benthic ecologies found within this member exceeded that typical of the pre-extinction interval (Fig. 7c). The return to sandy facies also marks the return of pedically attached epifaunal brachiopods and a greater level of bioturbation, with Chondrites, Rhizocorallium, Taenidium and large (25 mm diameter) Thalassinoides being reported (Martin 2001).

Following the temporary decrease in SID values in the serpulid-rich beds towards the top of the Grey Sandstone Member, there is an increase to stable and high SID values of c. 0.82. These are comparable with the pre-extinction interval and indicate a benthic macrofossil assemblage with an even representation of species (Fig. 6a). These high values of SID occur alongside the greatest values of species richness, which increase rapidly at the base of the Yellow Sandstone Member to surpass the pre-extinction values. Around 30 benthic species occur in samples towards the top of the Yellow Sandstone Member, which is c. 7 myr after the extinction event (Fig. 5). The Yellow Sandstone Member also contains the greatest number of ecologies of the study interval, with 11 recorded. These include all those present prior to the extinction event (Fig. 7). Within the Yellow Sandstone Member, bulk samples consistently reveal a greater benthic diversity than foreshore sampling (Fig. 5c). This trend is surprising and is explained by the mouldic preservation and small size of many of the specimens, especially gastropods, which were only identifiable after silicon rubber casting and inspection under a microscope. The means that much of the diversity was not appreciated during in-field species counts.

Of the 58 benthic species present in the recovery interval (excluding indeterminate architectibranch and cerithioid gastropods and crustacean limbs), only three species are shared between the post- and pre-extinction assemblages – Pleuromya costata, Oxytoma inequivalvis and Meleagrinella substriata – with only P. costata having a prolonged period of absence from the basin. A further 19 taxa are newly occurring species in the basin, but belong to genera that occupied the basin prior to the extinction. Twenty-six species are of genera new to the basin within the study interval.

The post-extinction interval within the Cleveland Basin is here divided into three phases or intervals. The first of these three intervals, the immediate post-extinction interval, spans c. 2 myr and, as this interval is represented by low-diversity, high-abundance assemblages with minimal tiering occupation, it could be viewed as a drawn-out survival phase or recovery phase one, following the nomenclature of Kauffman and Erwin (1995) and Twitchett (2006), respectively. It may instead be better to regard this interval as one dominated by low-oxygen specialist communities rather than the onset of true recovery. B. radiata and P. dubius are often described as opportunistic in their ecology, but an opportunist should be able to colonize a variety of habitats, whereas P. dubius in the Cleveland Basin is largely restricted to the organic-rich facies —that is, the Mulgrave Shale Member and the Sulfur Bands of the Grey Shale Member (Little 1995; Caswell and Coe 2013). This facies-dependent model for this species is additionally supported by the decreased abundance of P. dubius and its eventual disappearance with improving seafloor ventilation in the Alum Shale Member (Caswell and Coe 2013; Danise et al. 2015).

Phases two and three of the Twitchett (2006) recovery model feature an increase in ecological tiering and thereby an increase in ecological diversity. Our early phase recovery interval includes these same features, with the return of shallow and deep infaunal burrowers and erect epifauna. This provides a full suite of the benthic tiering structure c. 2.7 myr after the extinction event. However, throughout this interval, each of these tiers is represented by only a few species or a single species, or even a single specimen, and the absence of any of these from an individual sample horizon means that the number of tiers fluctuates. Shortly after the re-establishment of this structure, the diversity decreased and the tiering structure collapsed again during the deposition of the Peak Mudstone Member. However, we did not see a return to Pseudomytiloides-dominated assemblages. This interval of decreased diversity is also coincident with a decrease in exposure quality, but because the pattern of diversity was seen in both the bulk collections and the foreshore collections, it is believed to be a genuine signal. A similar trend was recorded by Hallam (1987) for European bivalve species richness per ammonite zone: after an initial increase during the Bifrons Zone, the diversity decreases in the Variabilis Zone.

Our late phase recovery interval features some of the characteristics of recovery stage 4 of Twitchett (2006), with a highly diverse and evenly represented assemblage present for both body and trace fossils. Ecological diversity attained pre-extinction values in the Grey Sandstone Member c. 6 myr after the extinction event. This was followed by the filling-out of most guilds, ensuring their presence in each sample horizon and thereby largely removing the fluctuations in ecological diversity seen in the previous phase. This pattern is an important facet of recovery and has been noted during recovery from the Cretaceous–Paleogene, Permian–Triassic and end-Triassic extinction events (Foster and Twitchett 2014; Alvarez et al. 2019; Atkinson and Wignall 2019). During this phase, the species richness exceeded that of the pre-extinction values and we therefore consider this to represent a fully recovered assemblage.

There were, however, several differences in the structure of the recovered assemblage when compared with the pre-extinction interval in the Cleveland Basin. Epifaunal carnivores (e.g. crustaceans and architectibranch gastropods) were represented by three specimens prior to the extinction event, but were of greater significance in the Yellow Sandstone fauna, with up to 8% of fossils from this member belonging to this ecology. For the gastropods, the pre-extinction interval was dominated by vetigastropods with a few caenogastropods, but the upper Toarcian strata were roughly equally dominated by heterobranchs and caenogastropods, with no vetigastropods (Ferrari et al. 2020). Within the study interval, crinoids did not regain their former diversity, whereas epifaunal grazers underwent a burst of diversity in the Yellow Sandstone Member that far exceeded the diversity in pre-extinction times. In addition, semi-infaunal suspension feeders occurred at higher relative abundances in the Yellow Sandstone Member than in the pre-extinction interval. This altered structure is also reflected in the relative diversity of the various bivalve orders.

The greatest difference between the pre-extinction interval and the Blea Wkye Sandstone assemblages are that the latter have a comparable lack of diversity within the Pectinida and an increase in Pholadida, despite deep infaunal suspension-feeding organisms being rarer in the recovered assemblages (Fig. 8). Solemyida is lacking from the recovered assemblage, although this is not considered significant as these bivalves were not a part of the normal pre-extinction communities, only being represented by a single species (Nucinella sp.) straddling the extinction event. Conversely, several species of bivalves belonging to Trigoniida arrived in the basin during the recovery interval. Little (1995) and Danise et al. (2013) did not find any trigonids in the pre-extinction interval, although rare Liotrigonia lingonensis have been reported from the Spinatum Zone (Morris et al. 2019). Trigonids go on to become an important group in the Mid-Jurassic through to the Late Cretaceous (Francis and Hallam 2003). Their paucity prior to the Toarcian extinction event has been noted across much of the British Lower Jurassic (see Supplementary material in Atkinson and Wignall 2020), suggesting that they were still recovering from the end-Triassic extinction event at the onset of the Toarcian event.

The fossil assemblage seen in the upper Toarcian of the Cleveland Basin represents the establishment of typical Mid-Jurassic faunas, with the presence of the bivalves Pteria plana, Goniomya literata, Myophorella striata, and Actinostreon marshi, as well as the first representatives of the gastropod families Aporrhaidae, Procerithiidae and Gordenellidae. These taxa all became major faunal components of the Mid-Jurassic (Ferrari et al. 2020).

Palaeo-redox indicators (Re/Mo and Mo/total organic carbon and $δ$^98/95Mo) suggest that benthic oxygen levels within the Cleveland Basin remained low until the Bifrons Zone, when changes in ocean circulation brought oxygenated Arctic waters south (McArthur et al. 2008; Danise et al. 2015; Caswell and Frid 2017; van de Schootbrugge et al. 2019). This accounts for the >2 myr gap between the extinction and the appearance of D. ovum. This interval is also characterized by a marked change in the primary producers, with calcareous nannoplankton and dinoflagellates being replaced by blooms of prasinophyte green algae (Palliani et al. 2002; Slater et al. 2019). Dinoflagellates begin to increase in relative abundance from the base of the Alum Shale Member, with diversity increasing up through the remaining Whitby Mudstone Formation and into the Blea Wyke Sandstone (Riding 1984; Slater et al. 2019; Danise et al. 2022). Dinoflagellates are regarded as a better food source for suspension-feeding organisms than green algae as a result of their size and nutrient content (Brown et al. 1997; von Elert et al. 2003; Wacker and von Elert 2008). Their demise and protracted absence, coupled with persistently low oxygen levels, contributed to the low diversity of suspension-feeding organisms in the Cleveland Basin during our immediate post-extinction interval (Danise et al. 2022). This relationship can further be evidenced from the comparably rapid benthic recovery seen in Toarcian successions of the Iberian Range (Spain), where there is no evidence for the development of an anoxic water column and benthic recovery occurred within the Serpentinum Zone (Danise et al. 2019).

In the upper portions of the Alum Shale Member, the Th/U ratios indicate increased levels of benthic oxygen (Martin 2001), this being around the time of the first appearance of D. ovum within the basin. Caswell and Dawn (2019) suggested Dacryomya was tolerant of low-oxygen conditions, allowing it to act as an initial colonizer and to take advantage of organic-rich sediments in areas that were previously anoxic. Dacryomya ovum occurs in a narrow stratigraphic interval, perhaps reflecting this requirement for a high organic content within the sediments, a feature needing the correct balance of productivity and low oxygen levels. Caswell and Dawn (2019) suggested that D. ovum was an environmental engineer; the action of their burrowing may have acted to oxygenate the substrate, enabling other organisms to subsequently colonize it. It is unclear how significant Dacryomya was in the role of opening niches, although its first appearance was shortly followed by the deep infaunal bivalve Gresslya and then subsequently crinoids – neither of these would have been excluded by anoxic sediments alone.

A further c. 5 myr had elapsed from the first occurrence of D. ovum to the fullness of recovery. This stands in contrast with Harries (1999), who, for the Cretaceous–Paleogene and Cenomanian–Turonian extinctions, found that it was the duration of the environmental perturbation that controlled the duration of the ‘survival interval’ and, once this had passed, recovery occurred rapidly. Following the delay to recovery caused by frequent bottom water anoxia during the deposition of the Mulgrave Shale Member, the sluggish recovery thereafter in the Cleveland Basin is attributable to the intervals represented by the Peak Mudstone and Fox Cliff Siltstone members. Species richness was low throughout the Peak Mudstone Member, which is notable for its dark coloration, phosphatic nodules and pyritization of both wood and small simple burrows. These features may indicate that the Cleveland Basin was again experiencing low levels of oxygen close to the seafloor during this interval. If so, this was not accompanied by a temporary rise in sea-level (Knox 1984).

A later pulse of anoxia would account for the loss of the shallow infaunal species V. literata and D. ovum and the temporary absence of deep infaunal bivalves until the topmost metre of sediment in this unit. Apart from low-resolution sulfur isotope data, which suggest a continued long-term decrease in pyrite burial in the basin (Newton et al. 2011), there are currently no high-resolution redox data available for this interval. However, by the Dispansum Zone, the Th/U values suggest that reoxygenation had occurred (Martin 2004), this also being a time of significant advance in recovery.

The duration of recovery was probably influenced by changes in sea-level. It is notable that the extinction is bracketed by shallow water facies (the Staithes Sandstone, the Cleveland Ironstone and the Blea Wyke Sandstone formations), which host a higher diversity than the intervening deeper water mudstones of the extinction event and the succeeding low-diversity interval (the immediate post-extinction interval). This is not to say that the extinction and recovery are merely an expression of an ordinary water depth–benthic marine diversity gradient with species tracking their preferred habitats (Holland and Patzkowsky 2015; Holland 2020). Had this been the case, then we would expect a far higher number of returning species in the recovered assemblage, rather than only one.

In the intervening 7 myr, it is unlikely that a great many bivalve species would have become extinct as a result of the background extinction rates alone because Jurassic bivalve species have a typical range of 15 ammonite zones (Hallam 1975). It is, however, more likely that recovery occurred faster in the shallower water regions of the Cleveland Basin (Caswell and Coe 2012), with species then tracking their preferred habitat basinwards during shallowing in the mid- and late Toarcian. This could account for the exceptionally long duration of the local recovery. However, the full extent of the control of sea-level on the recovery can only be realized by studying a series of rock sections that would provide a depth transect through the Cleveland Basin across the Toarcian. If such a study is undertaken, then it may reduce the temporal duration of the recovery (Holland 2020).

This idea is supported by the findings of Danise et al. (2019) from the Iberian Range of Spain. Their study found that the deeper sections were slightly slower to recover (the Falciferum Subzone rather than the Elegantulum Subzone, equivalent to the Exaratum Subzone; Danise et al. 2019). However, it is important to note that these Spanish sections show no evidence of water column anoxia and the interplay of anoxia and sea-level probably caused the overall delay in the Yorkshire successions.

Temperature has been proposed as the main driver of extinction in sections that do not show evidence of the TOAE, but still record the early Toarcian mass extinction event (Danise et al. 2019). In the Cleveland Basin, temperature curves have been generated from oxygen isotopes from (mostly) belemnite calcite (Korte and Hesselbo 2011, their supplementary material; Korte et al. 2015). These records clearly show the abrupt warming associated with the early Toarcian mass extinction event and the temperatures remained elevated throughout the Serpentinum Zone, but the extent by which the subsequent cooling impacted the duration of the recovery is unclear. Temperatures had begun to decrease by the end of the Serpentinum Zone, but there is a large degree of scatter in the oxygen isotope data, which is attributable to seasonality and the habitat effects of the different belemnite species used. Substantial cooling does not take place until the Aalenian (Mid-Jurassic) (Korte et al. 2015, their supplementary material).

Our findings from the Cleveland Basin add to a developing picture that post-extinction recovery is a mosaic that is highly variable between regions and is also dependant on the criteria used to assess recovery. A range of possible recovery durations could be argued for the Cleveland Basin. Based on ecological tiering alone, the recovery spanned <2.7 myr. For ecological diversity to be restabilized, this duration massively increases to c. 6 myr and, using the traditional species counts as the recovery criteria, it is longer still at c. 7 myr. A delayed recovery was also recorded in the trace fossil communities of the Western Tethys Ocean, as evidenced from the Central Western Carpathians of Slovakia (Müller et al. 2020). Like the Cleveland Basin, the Serpentinum Zone is marked by a reduced faunal assemblage with the only trace fossils being small and thin and recovery only beginning in the Bifrons Zone. Full pre-extinction diversity is not reattained within their study section, meaning that full recovery exceeds 3.4 myr (Müller et al. 2020).

These durations are greater than those in sections from the Sub-Boreal Realm, such as the Iberian and Basque–Cantabrian basins (García Joral et al. 2011; Danise et al. 2019). In these regions, elevated temperatures may have been the main driving mechanism behind the extinction in the absence of lingering water column anoxia, as seen in the NW European epicontinental basins. The brachiopods faunas recovered by the Variabilis Zone (3.4–4.9 myr after the extinction event) in the Basque–Cantabrian Basin (García Joral et al. 2011). Repopulation for brachiopods was faster still in the Iberian Basin, occurring in the Serpentinum Zone (<2 myr after the event), and bivalves were little affected by the extinction (Danise et al. 2019).

This pattern of regional differences has similarities with the recovery following the Permian–Triassic mass extinction event, the pace of which was closely related to regional differences in benthic oxygenation. Well-oxygenated Lower Triassic sections commonly feature a level of recovery that would otherwise not have been anticipated until the Mid-Triassic (Twitchett et al. 2004; Beatty et al. 2008; Hofmann et al. 2011; Chen et al. 2012; Dai et al. 2018). Such regions may also have acted as refuges or origination centres from which species could migrate into other regions once they became habitable again, an idea that has been proposed for migration along the Hispanic Corridor (Hallam 1983; Aberhan 2002).

Large, multi-regional studies conducted using a time-binned approach have previously suggested that recovery from the Toarcian extinction was extremely protracted, extending well into the Mid-Jurassic (Hallam 1996; Ruban 2004). As more local/regional studies are conducted, this picture is being changed to a more complex mosaic of recoveries occurring at variable rates. For the Cleveland Basin, our longest estimate for the duration for complete recovery is c. 7 myr, which would place the Toarcian recovery as being slower than the recovery from the end-Trassic and Cretaceous–Paleogene mass extinctions and on a par with some of the estimates for the Permian–Triassic recovery (Table 1). This is testimony that the intensity of the extinction event is by no means an indication of the duration of recovery (Solé et al. 2002). Despite the differences in the magnitudes of the end-Permian and early Toarcian extinction events, there are parallels: both are considered to have been volcanically triggered and involved the spread of oxygen-restricted facies and elevated temperatures. The continuation of harsh environmental conditions has provided one of most plausible theories for why the end-Permian extinction had a protracted interval of low diversity (e.g. Hallam 1991; Foster et al. 2018) – a theory that may also be transferable to the Toarcian, with further modulation resulting from changes in sea-level.

This study documents the full extent and nature of biotic recovery from the early Toarcian mass extinction event within the Cleveland Basin. Benthic oxygen levels remained low following the extinction event, allowing for specialist assemblages tolerant to low oxygen levels to dominate during an interval of high sea-level. These disappeared with improving ventilation and recovery proper was marked by the expansion of infaunal tiers (c. 2 myr post-extinction), followed by epifaunal tiers (c. 2.7 myr post-extinction). Recovery progressed slowly with the temporary absence of infaunal molluscs in the Peak Mudstone Member (straddling the Variablis and Thoarsense zones), which may have been caused by a return to oxygen-restricted environments. With a continued fall in sea-level and a return to silt- and sand-dominated deposition in the Grey Sandstone Member, the ecological richness matched that of the pre-extinction interval, although many of these ecologies contained only a few species. It was not until the upper Yellow Sandstone Member, when a dramatic increase in species richness occurred, that all the previously present ecologies reappeared and increased diversity facilitated a filling-out of these life modes. The Yellow Sandstone benthic assemblage records the establishment of a typical Mid Jurassic-type assemblage.

Recovery within the Cleveland Basin took at most 7 myr, a duration on par with estimates of recovery rates from the largest mass extinction events of the Phanerozoic. However, recovery within the basin was likely to have been strongly influenced by local sea-levels and the continuation of harsh environmental conditions after the extinction event. To fully appreciate the duration and dynamics of the recovery, a depth transect is required within a single basin, something that is currently unavailable for the Cleveland Basin.

The authors thank Brian Atkinson and Karolina Zarzyczny for their assistance during field collections. We thank Kevin Page, Jean Guex, Andreas Schmidt, Paul Taylor, Robin Knight, Mariel Ferrari, Martin Munt, Carrie Schweitzer and Mike Simms for assistance with resolving some of the taxa within this study. We also thank Silvia Danise and Ian Boomer for providing helpful comments on our paper. Thanks also to Rebecca Bennion for her support and encouragement and for providing comments on an earlier version of this paper.

JWA: conceptualization (equal), data curation (lead), formal analysis (lead), investigation (lead), writing – original draft (lead); CTSL: conceptualization (equal), supervision (lead), writing – review & editing (lead); AMD: formal analysis (supporting), software (lead), writing – review & editing (supporting).

This work was funded by a Geological Society of London undergraduate research bursary in 2013, and The University of Leeds (NERC SPHERES DTP NE/L002574/1).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

All data generated or analysed during this study are included in this published article (and if present, its supplementary information files).