Abstract

Reduction of body size is a common response of organisms to environmental stress. Studying the early Toarcian succession in the Lusitanian Basin of Portugal, we tested whether the shell size of benthic marine communities of bivalves and brachiopods changed at and before the global, warming–related Toarcian oceanic anoxic event (T-OAE). Statistical analyses of shell size over time show that the mean shell size of communities decreased significantly before the T-OAE. This trend is distinct in brachiopods and is caused by larger-sized species becoming less abundant over time, whereas it is not significant in bivalves, suggesting a decoupled response to environmental stress. Reductions in shell size precede the decline in standardized sample-level species richness associated with the early Toarcian extinction event. Such decreases in the shell size of marine invertebrates, well before the onset of biodiversity change, suggest that reductions in body size more generally may be a precursor of a subsequent loss of species and turnover at the community level caused by climate change. Sedimentological evidence is against hypoxia as a driver of extinction and the preceding size decrease in the brachiopod fauna in the studied succession, although low oxygen levels are widely held responsible for elevated early Toarcian extinction rates globally. Reduction of mean shell size in brachiopods but stasis in bivalves is difficult to explain with ocean acidification, because experimental work shows that brachiopods can be resilient to lowered pH, albeit long-term metabolic costs and potential evolutionary adaptations are unknown. Rising early Toarcian temperatures in the Lusitanian Basin seem to be a plausible factor in both diversity decline associated with the T-OAE and the preceding reductions in mean shell size, because thermal tolerances in modern bivalves are among the highest within marine invertebrates.

Introduction

Global warming and stressors associated with climate change—in particular hypoxia, oceanic acidification, and hypercapnia—have raised concerns over the degradation and transformation of marine ecosystems. While ongoing global warming and its related effects are attributed to CO2 release from the combustion of fossil fuels, climate-related crises in the geologic past are commonly associated with the release of greenhouse gases during episodes of intense volcanism (e.g., Wignall et al. 2005 and references therein). Oceanic acidification is a direct consequence of elevated CO2 through the absorption and dissolution of CO2 and SO2 in marine waters, reducing oceanic pH (Jenkyns 2010; Kelly and Hofmann 2013). Hypoxia is related to water-column stratification due to sluggish circulation and/or increased nutrient input and productivity associated with warming (e.g., Bijma et al. 2013; Breitburg et al. 2018). Yet the relative role of these factors is still uncertain in both present and past biotic crises.

One such crisis interval that has received increased attention is the early Toarcian oceanic anoxic event (T-OAE; ~183 Ma). It marks a global, second-order mass extinction that affected both benthic and nektonic marine species, as well as terrestrial species (e.g., Wignall and Bond 2008; Morten and Twitchett 2009; Mattioli et al. 2009; Martindale and Aberhan 2017). Most likely the ultimate cause of the biotic crisis was intense volcanic activity during the placement of the Karoo-Ferrar Large Igneous Province (Pálfy and Smith 2000), possibly coupled with warming- related methane release to the atmosphere from the ocean (Hesselbo et al. 2000) and from terrestrial sources (Them et al. 2017). However, the proximate killing mechanisms are still poorly understood. Extinctions are mostly attributed to widespread anoxia (Jenkyns 1988, 2010; Hallam and Wignall 1997; Aberhan and Baumiller 2003; Wignall and Bond 2008; Danise et al. 2015), but the roles of ocean acidification (Trecalli et al. 2012), rising temperature (e.g., Suan et al. 2010; Gómez and Goy 2011; Danise et al. 2015), and the potentially synergistic and geographically variable effects of multiple drivers are less clear.

In this study we focus on changes in the body size of benthic marine organisms before the T-OAE. There is evidence of environmental perturbations and temperature fluctuations related to a long-term shift from icehouse to greenhouse conditions starting during the late Pliensbachian, and these factors may have begun to affect organisms before the main warming episode of the T-OAE (Suan et al. 2010; see “Discussion” for more details). Body size is a key parameter related to many aspects of an organism’s ecology, behavior, and biology (e.g., growth, metabolism, feeding, and fecundity) and can be strongly affected by environmental stress (e.g., Daufresne et al. 2009). Identifying the direction, magnitude, timing, and biotic selectivity of body-size change can provide insights on the relative importance of different stressors at times of faunal crisis. Size reduction related to climate change has been observed in living and fossil taxa, both terrestrial and marine (e.g., He et al. 2007, 2017; Gardner et al. 2011; Ohlberger 2013; O’Gorman et al. 2017). It has been suggested that body-size change could be a third universal response to warming along with changes in the geographic distribution of species and changes in phenological events (Gienapp et al. 2008; Daufresne et al. 2009; Gardner et al. 2011). Studies on variations in body size related to extinction episodes usually have focused on the aftermath of a crisis (e.g., Lockwood 2005; Morten and Twitchett 2009; Huang et al. 2010). Yet few studies have focused on changes in body size during precrisis times (e.g., He et al. 2007, 2010; Zhang et al. 2016; García Joral et al. 2018; Kiessling et al. 2018).

Here, we briefly review current knowledge of the physiological responses of marine ectotherms to environmental stress, particularly regarding growth and body size. We use this physiological information to state three specific hypotheses about changes in body size as a response to environmental stress before the onset of the T-OAE. Subsequently, we test these hypotheses by analyzing quantitative field-based occurrence data of macrobenthic marine communities in the Lusitanian Basin of western Portugal.

Physiological Response to Environmental Stress and Expected Effects on Body Size

Body-size reduction is one common organismic response to temperature-induced stress and hypoxia and can be explained by the “thermal window” concept. Thermal windows define the temperature range upon which proper functioning of an organism depends, and within such a range, there is an optimum temperature for biological functions such as growth. Thermal windows differ between species and can shift or narrow as an adaptation to stress (Pörtner 2008; Pörtner and Farrell 2008; Day et al. 2018). With warming, oxygen supply cannot match the increased demand and consumption, such that the aerobic scope is reduced until passive tolerance (metabolic depression) ensures temporary survival by reducing biological functions, including growth (Pörtner and Farrell 2008; Breitburg et al. 2018). Warming beyond critical threshold temperatures finally leads to a collapse of physiological functions with lethal effects (Pörtner and Farrell 2008).

Tolerance to warming and hypoxia may also be reduced by increasing CO2 levels in the water, resulting in ocean acidification and hypercapnia leading to metabolic depression and the interruption of growth (Pörtner et al. 2005; Widdicombe and Spicer 2008; Fabry et al. 2008; Bijma et al. 2013). Acidification affects calcification, but the exact response of growth processes, and hence body size, is difficult to predict, because growth is not constant but varies with temperature, food supply, and predation pressure (Gazeau et al. 2013). The effects of lowered pH for calcifying invertebrates also depend on the organisms’ ability to regulate pH at the site of calcification, the extent of organic-layer coverage of the external shell, and biomineral solubility (Ries et al. 2009; Parker et al. 2013).

The physiological effects of warming, hypoxia, hypercapnia, and acidification are related to one another and may be exacerbated if these stressors act in synergy (Pörtner et al. 2005; Byrne and Przeslawski 2013; Kelly and Hofmann 2013), although their combined effects are still poorly understood. Notwithstanding, reductions in body size have been observed with increased temperature (Daufresne et al. 2009), reduced oxygen concentrations (Diaz and Rosenberg 2008), and ocean acidification (Widdicombe and Spicer 2008; Parker et al. 2013). Everything else being equal, large-sized species are expected to be more strongly affected by growth restrictions than smaller-sized species because of their higher metabolic demands, which would be selected against under stressful conditions (He et al. 2007; Daufresne et al. 2009; Genner et al. 2009; O’Gorman et al. 2017). Consequently, we expect a decrease in body size as a response to increasing environmental stress, with more pronounced patterns in large-sized species.

Assuming that physiological principles play out on both the short-term and geologic timescales (see “Discussion”), we test the following hypotheses: (1) The body size of species and communities decreased before the early Toarcian extinction event; (2) large-sized species are more affected than small ones; and (3) reduction in body size occurred earlier than diversity loss and changes in faunal composition. We find support for each of these three hypotheses and discuss to what extent change in shell size before the main phase of faunal loss and extinction may have been caused by the same factor(s).

Studied Sections and Depositional Environments

Fieldwork was performed near Coimbra, Portugal, at the composite sections of Fonte Coberta (40°03′36.5″N, 8°27′33.4″W) and Rabaçal/Maria Pares (40°03′08.0″N, 8°27′30.5″W), located ca. 1 km apart from each other. These sections are representative of the Early Jurassic succession in the Lusitanian Basin, are highly fossiliferous, and have a well-defined biozonation based on ammonites (e.g., Mouterde et al. 1964–1965; Comas-Rengifo et al. 2013), nannofossils (Ferreira et al. 2015), and dinoflagellate cysts (Correia et al. 2018).

The ca. 28-m-thick studied succession (Fig. 1) spans the interval from the Pliensbachian/Toarcian boundary (base of the Polymorphum Zone = Tenuicostatum Zone of the Submediterranean Province) to the middle of the Levisoni Zone (= Serpentinum Zone). Lithostratigraphically, the studied succession is included in the São Gião Formation, which can be divided into three members (e.g., Duarte and Soares 2002):

  1. Marly limestones with Leptaena Fauna (hereafter referred to as the first member): an alternation of decimeter-thick grayish marlstones and marly mud- to wackestones rich in benthic macroinvertebrates representing the Polymorphum Zone, which comprises the Mirabile Subzone and the Semicelatum Subzone (Fig. 1).

  2. Thin nodular limestones (TNL member): gray to brownish, thin-bedded micritic carbonates and subordinate marlstones representing the lower Levisoni Zone. The transition from the first member to the TNL is marked by a color change from grayish to brownish (Miguez-Salas et al. 2017) and represents a facies change in the studied section and elsewhere in the Early Jurassic of the Lusitanian Basin (Duarte 1997; Duarte et al. 2007; Pittet et al. 2014). Body fossils are extremely rare in the lower part of this member, but become more common up-section. Bioturbation, mainly represented by horizontal networks of Thalassinoides and ferruginous nonbranching tubular burrows, is pervasive (Miguez-Salas et al. 2017; Rodríguez-Tovar et al. 2017).

  3. Marls and marly limestones with Hildaites and Hildoceras: a decimeter- to meter-thick alternation of marls and marly limestones with moderately diverse nektonic and benthic fauna. This member corresponds to the interval from the middle Levisoni Zone up to the middle part of the upper Bifrons Zone.

Figure 1.

Stratigraphic log of the early Toarcian succession of the composite sections at Fonte Coberta and Rabaçal near Coimbra, Portugal, with the lithostratigraphic and biostratigraphic zonations of Mouterde et al. (1964–1965), Duarte and Soares (2002), and Comas-Rengifo et al. (2013). Stratigraphic ranges of species are based on recorded occurrences (black dots) ordered by last appearances and separately for bivalves, brachiopods, and gastropods. The shaded area marks the extent of the Toarcian oceanic anoxic event (T-OAE) as defined by carbon isotope data. The dashed horizontal line represents the level where faunal loss is severe in terms of both species richness and fossil abundance. Dashed lines within the stratigraphic ranges are used when the extent of the range is uncertain. The star-shaped symbols indicate sampling levels. Abbreviations for lithology: M, marlstone; CM, calcareous marlstone; ML, marly limestone; L, limestone.

Figure 1.

Stratigraphic log of the early Toarcian succession of the composite sections at Fonte Coberta and Rabaçal near Coimbra, Portugal, with the lithostratigraphic and biostratigraphic zonations of Mouterde et al. (1964–1965), Duarte and Soares (2002), and Comas-Rengifo et al. (2013). Stratigraphic ranges of species are based on recorded occurrences (black dots) ordered by last appearances and separately for bivalves, brachiopods, and gastropods. The shaded area marks the extent of the Toarcian oceanic anoxic event (T-OAE) as defined by carbon isotope data. The dashed horizontal line represents the level where faunal loss is severe in terms of both species richness and fossil abundance. Dashed lines within the stratigraphic ranges are used when the extent of the range is uncertain. The star-shaped symbols indicate sampling levels. Abbreviations for lithology: M, marlstone; CM, calcareous marlstone; ML, marly limestone; L, limestone.

This hemipelagic Pliensbachian–Toarcian succession was deposited in a low-energy environment on a middle to distal homoclinal ramp dipping toward the northeast (Duarte 1997, 2007; Duarte et al. 2007). In particular, the monotonous succession of micritic limestones and marlstones of the interval studied herein in detail (see below) suggests that the depositional environment remained below storm-wave base without any obvious changes in water depth. Thus, while a decrease in shell size with depth has been reported, for example, in modern terebratulid brachiopods (Peck and Harper 2010), any such size changes in our material could not be explained by variations in water depth. Moreover, the shells are evenly distributed within the sedimentary rocks. Lack of size sorting, absence of preferential orientations of shells, and lack of faunal amalgamation suggest that shell transport was insignificant and that the distribution of shell size within samples was not biased by physical agents such as currents and waves.

Unlike well-studied early Toarcian successions elsewhere, such as those in England (e.g., Wignall and Bond 2008; Danise et al. 2015) and Germany (e.g., Röhl et al. 2001), the black shales typical for the T-OAE are not developed. The interval corresponding to the T-OAE is mostly represented by the micritic carbonates and marlstones of the TNL member (Fig. 1). The position of the T-OAE in our section is defined by the globally recorded negative carbon isotope excursion (e.g., Hesselbo et al. 2007; Jenkyns 2010), because black shales are absent and the total organic carbon content is low (Duarte et al. 2005). The T-OAE starts when the isotopic values begin to decrease and terminates when the values have returned to background level and thus, in our section, spans the uppermost Polymorphum Zone and the lower part of the Levisoni Zone (Fig. 1). The T-OAE has been estimated to last between 300 kyr to 900 kyr using geochronology, astronomical calibrations, and biostratigraphy (Suan et al. 2008, 2015; Boulila et al. 2014; Sell et al. 2014). Using the timescale provided in Suan et al. (2015: Fig. 5), we estimated the Polymorphum Zone, which is the time interval focused on in our study, to last ~900 kyr.

Materials and Methods

This study uses two types of data: (1) quantitative field data with abundance counts of specimens are used for analyzing trends in shell size and diversity; and (2) occurrence data from the Paleobiology Database (www.paleobiodb.org) and the literature are applied to reconstruct the paleolatitudinal distribution of species.

In the field, we sampled the limestone beds of the succession at the sampling levels indicated in Figure 1. Sampling intensity was standardized by collecting the same amount of bulk rock, ca. 15 kg per sample, from one spot. Through quantitative bed-by-bed sampling of the three stratigraphic members, we recorded a total of 1321 specimens belonging to 58 taxa of brachiopods, bivalves, and gastropods (Fig. 1). While our study focuses on benthic macroinvertebrates, we also recorded occurrences of ammonoids, belemnites, ichnofossils, and wood fragments. As far as permitted by preservation, the benthic fauna was identified to the species level using the relevant literature (e.g., Alméras 1994; Gahr 2002; Baeza-Carratalá 2013; Comas-Rengifo et al. 2013).

Overall, preservation quality of the benthic fauna is moderate. Bivalves are commonly preserved as internal molds, more rarely as external molds, and occasionally with parts of the shell. As an exception, the calcitic shells of the plicatulid Harpax spinosa and those of oysters are preserved as complete single valves or articulated specimens. Different groups of brachiopods are variably preserved: terebratulids and rhynchonellids exhibit recrystallized shells, while spiriferinids generally preserve a thin shell layer and are filled with sediment. Preservation is poorest in gastropods, which consist of incomplete internal molds. Because gastropods are rare and fragmented, we excluded them from all quantitative analyses, which thus examine only the bivalve–brachiopod community.

Shell size was measured with calipers to the nearest 0.1 mm. Measurements were inferred when only a small fraction of the shell was missing, while incomplete specimens were disregarded. The size of a specimen was calculated as the log2 of the geometric mean of shell length (L) and height (H) in bivalves and shell width (W) and length (L) in brachiopods (LogGeoMean), which is a good proxy for body size (Kosnik et al. 2006). For calculation of the geometric mean (in millimeters), the following equations were applied:  
$${\rm GeoMea}{\rm n}_{{\rm Bivalves}} = \sqrt {L{\rm {^\ast}}H} $$
1
 
$${\rm GeoMea}{\rm n}_{{\rm Brachiopods}} = \sqrt {L{\rm {^\ast}}W} $$
2

To test our hypotheses, we focused on the precrisis interval (see “Results”), which covers the lowermost 7.5 m of the section. Three samples are from the Mirabile Subzone, and 12 samples are from the Semicelatum Subzone (Fig. 1). For this precrisis interval we recorded 942 occurrences of 35 taxa of bivalves and brachiopods, of which 588 were measured (400 brachiopods and 188 bivalves; see Supplementary Material). This data set was used to construct size–frequency histograms and to analyze size trends in individual species/genera (discussed later). All other analyses of shell size were performed on an extended database in which a shell size value was assigned to each occurrence. For individuals lacking measurements, the mean shell size of the respective species in the sample was allocated. In cases in which no size measurements were available for a species in a sample, we used the mean size of this species in the sample directly below and above or, when this option was not feasible, the mean size of the species in all samples of the first member. After deletion from the data set of a few taxa that had insufficient size data, the extended data set consists of 830 records of shell size (571 brachiopods and 259 bivalves; see Supplementary Material), that is, 71% of the size data are from actual measurements of specimens and 29% are inferred following the procedure described above.

Shell-size analyses were performed at the community level (i.e., all individuals of brachiopods and bivalves of a faunal sample), separately for all brachiopods pooled and for all bivalves pooled, and for individual taxa. Changes in community-level shell size through time were investigated using the mean of all individuals of a sample irrespective of taxonomic identity. We consider this measure as the average shell size of the time-averaged relics of ancient bivalve–brachiopod communities. We also calculated this measure separately for brachiopods and bivalves. Ideally, tests for temporal changes in shell size were performed for each species separately. In the majority of cases this was hampered by low numbers of samples and/or specimens per species. Therefore, we illustrate the time series of mean size and maximum size of five taxa for which sample sizes were deemed sufficient, that is, a taxon was present in at least 10 samples with at least 30 specimens.

Bivalves and brachiopods are double-valved organisms. To determine the number of individuals from counts of specimens, we initially recorded left, right, and articulated valves for bivalves, and dorsal and ventral valves and double-valved specimens in brachiopods separately. For the final analyses, however, each specimen was counted as one individual. This approach is justified, because single valves of the same species in any of the samples did not match in shell size and therefore must represent different individuals.

To address our second hypothesis, that is, large species are more strongly affected by environmental stress than small species, each species was categorized as “smaller-sized” or “larger-sized.” To this end, we determined the mean shell size of each species in the pooled samples of the precrisis interval and used the median of these values to separate larger-sized species from smaller-sized species. For each of the two such defined subsets, the change in shell size over time was analyzed in the same way as described earlier. In addition, the abundance of individuals of larger-sized species per sample, expressed as their percentage relative to all individuals of a sample, was tracked over time. Both steps were performed for the bivalve–brachiopod communities as a whole and separately for brachiopods and bivalves.

To illustrate the direction of a trend in shell size, if any, we applied weighted LOESS smoothing. To statistically test for the existence of a trend in shell size, we fit several models of trait evolution, using the R package ‘paleoTS’ (Hunt 2006, 2015), to the time series of mean geometric mean at each level for the whole bivalve–brachiopod community, for brachiopods and bivalves separately, and for larger-sized and smaller-sized brachiopods and bivalves. The three models tested were stasis, random walk (URW), and directional trend (GRW). For each model and each time series, we report the corrected Akaike information criteria (AICc), the difference between the AICc and the minimum AICc (ΔAICc), and the Akaike weights. Following Burnham and Anderson (2003), we used a rule of thumb of ΔAICc < 2 to consider models that are not significantly less plausible than the best-fit model. Finally, to assess the adequacy of each model, package ‘adePEM’ (Voje 2018) was used to test the models (for autocorrelation, length of runs, fixed variance over time, and in the case of the stasis model, net evolution) by running a large number (here 10,000) of simulated time series using the parameters of the fitted models and checking whether they are likely to belong to the same distribution (here, the null hypothesis, with p > 0.05 being the qualifying significance). In addition, we applied Spearman’s rank correlation, as all these time series showed no autocorrelation, having passed the Box-Pierce test of autocorrelation (Box and Pierce 1970) as implemented in base R by function Box.test, and thus their values can be considered to be independent. All the results of the statistical tests are shown in Tables 13.

Table 1.

Results of the statistical test for autocorrelation (Box-Pierce test) and of the relative Hunt’s (2015) model fit estimates of the analyzed time series. Significant p-values from the Box-Pierce test and estimates of the best model(s) are marked in bold. For the autocorrelation test, p-values < 0.05 would mean that samples are autocorrelated. AICc, Akaike information criteria; ΔAICc, the difference between the AICc and the minimum AICc.

Analyzed group
(Figure in this study)
Directional changeRandom walkStasis
AICcΔAICcWeightAICcΔAICcWeightAICcΔAICcWeight
Community (Fig. 2A)3.3111.7410.2761.5700.666.2524.6820.064
Box-Pierce test0.199
Bivalves (Fig. 2B)18.8039.5540.00713.284.0310.1179.24900.876
Box-Pierce test0.596
Brachiopods (Fig. 2B)9.2981.7860.2817.51200.68713.6356.1230.032
Box-Pierce test0.073
Smaller-sized taxa (/)−12.8033.3140.087−16.11200.456−16.11200.456
Box-Pierce test0.309
Larger-sized taxa (/)6.7916.9950.0263.9744.1790.107−0.20400.867
Box-Pierce test0.892
Smaller-sized brachiopods (Fig. 3A)−8.8242.8030.110−11.62600.445−11.62600.445
Box-Pierce test0.490
Larger-sized brachiopods (Fig. 3A)11.5076.2840.0215.22200.4895.22200.489
Box-Pierce test0.970
Smaller-sized bivalves (Fig. 3B)19.2777.0900.02515.9933.8060.12712.18700.849
Box-Pierce test0.543
Larger-sized bivalves (Fig. 3B)29.1189.8820.00725.5206.2540.04119.23600.952
Box-Pierce test0.961
Analyzed group
(Figure in this study)
Directional changeRandom walkStasis
AICcΔAICcWeightAICcΔAICcWeightAICcΔAICcWeight
Community (Fig. 2A)3.3111.7410.2761.5700.666.2524.6820.064
Box-Pierce test0.199
Bivalves (Fig. 2B)18.8039.5540.00713.284.0310.1179.24900.876
Box-Pierce test0.596
Brachiopods (Fig. 2B)9.2981.7860.2817.51200.68713.6356.1230.032
Box-Pierce test0.073
Smaller-sized taxa (/)−12.8033.3140.087−16.11200.456−16.11200.456
Box-Pierce test0.309
Larger-sized taxa (/)6.7916.9950.0263.9744.1790.107−0.20400.867
Box-Pierce test0.892
Smaller-sized brachiopods (Fig. 3A)−8.8242.8030.110−11.62600.445−11.62600.445
Box-Pierce test0.490
Larger-sized brachiopods (Fig. 3A)11.5076.2840.0215.22200.4895.22200.489
Box-Pierce test0.970
Smaller-sized bivalves (Fig. 3B)19.2777.0900.02515.9933.8060.12712.18700.849
Box-Pierce test0.543
Larger-sized bivalves (Fig. 3B)29.1189.8820.00725.5206.2540.04119.23600.952
Box-Pierce test0.961
Table 2.

Results of the adequacy tests for the three analyzed models of change in shell size (directional change, random walk, and stasis). The p-values are provided for each test; p-values marked in bold indicate a test was passed, otherwise a test has failed. The p-values of the adequacy tests represent the portion, divided by 0.5, of the simulated test statistics that is larger/smaller than the test statistics calculated on the actual data. A test is passed if the value of the test statistic falls within the distribution range provided by simulated test statistics (Voje 2018).

Analyzed group (Figure in this study)Adequacy tests (p)
AutocorrelationRunsFixed varianceNet change
Community (Fig. 2A)
Directional change0.3670.7300.695NA
Random walk0.9780.9300.009NA
Stasis0.0070.9890.4580.007
Bivalves (Fig. 2B)
Directional change0.0170.0030.908NA
Random walk0.2890.0220.011NA
Stasis0.5780.1640.0580.039
Brachiopods (Fig. 2B)
Directional change0.1460.1900.824NA
Random walk0.4950.7920.201NA
Stasis0.0050.3470.6050.017
Smaller-sized taxa (NA)
Directional change0.3360.7800.499NA
Random walk0.4260.5940.248NA
Stasis0.4260.8600.4450.350
Larger-sized taxa (NA)
Directional change0.8620.5880.754NA
Random walk0.5610.7400.455NA
Stasis0.6330.8710.9980.250
Smaller-sized brachiopods (Fig. 3A)
Directional change0.2790.7750.839NA
Random walk0.6800.5280.600NA
Stasis0.2230.5710.6710.769
Larger-sized brachiopods (Fig. 3A)
Directional change0.8820.8370.956NA
Random walk0.4510.8960.805NA
Stasis0.5730.8730.5680.640
Smaller-sized bivalves (Fig. 3B)
Directional change0.3320.7660.122NA
Random walk0.2880.4560.005NA
Stasis0.3470.7380.0010.554
Larger-sized bivalves (Fig. 3B)
Directional change0.7380.3890.786NA
Random walk0.6440.4580.156NA
Stasis0.7110.5240.4160.180
Analyzed group (Figure in this study)Adequacy tests (p)
AutocorrelationRunsFixed varianceNet change
Community (Fig. 2A)
Directional change0.3670.7300.695NA
Random walk0.9780.9300.009NA
Stasis0.0070.9890.4580.007
Bivalves (Fig. 2B)
Directional change0.0170.0030.908NA
Random walk0.2890.0220.011NA
Stasis0.5780.1640.0580.039
Brachiopods (Fig. 2B)
Directional change0.1460.1900.824NA
Random walk0.4950.7920.201NA
Stasis0.0050.3470.6050.017
Smaller-sized taxa (NA)
Directional change0.3360.7800.499NA
Random walk0.4260.5940.248NA
Stasis0.4260.8600.4450.350
Larger-sized taxa (NA)
Directional change0.8620.5880.754NA
Random walk0.5610.7400.455NA
Stasis0.6330.8710.9980.250
Smaller-sized brachiopods (Fig. 3A)
Directional change0.2790.7750.839NA
Random walk0.6800.5280.600NA
Stasis0.2230.5710.6710.769
Larger-sized brachiopods (Fig. 3A)
Directional change0.8820.8370.956NA
Random walk0.4510.8960.805NA
Stasis0.5730.8730.5680.640
Smaller-sized bivalves (Fig. 3B)
Directional change0.3320.7660.122NA
Random walk0.2880.4560.005NA
Stasis0.3470.7380.0010.554
Larger-sized bivalves (Fig. 3B)
Directional change0.7380.3890.786NA
Random walk0.6440.4580.156NA
Stasis0.7110.5240.4160.180
Table 3.

Results of statistical tests of the correlation of shell size with time (sampling level) for different faunal groupings. Significant p-values are marked in bold. Tests on the larger- and smaller-sized groups of taxa were performed on the mean of the species means.

Analyzed groupFigureSpearman’s rank correlationMann-Whitney rank test
p-valueρp-value
Community2A<0.001−0.796NA
Bivalves2B0.123−0.418NA
Brachiopods2B<0.001−0.868NA
Smaller-sized taxaNA0.259−0.311NA
Larger-sized taxaNA0.267−0.319NA
Smaller-sized brachiopods3A0.630−0.136NA
Larger-sized brachiopods3A0.614−0.173NA
Smaller-sized bivalves3B0.8440.059NA
Larger-sized bivalves3B0.409−0.240NA
Lower half vs. upper half (community)5NANA<0.001
Lower half vs. upper half (bivalves)5NANA0.772
Lower half vs. upper half (brachiopods)5NANA<0.001
Analyzed groupFigureSpearman’s rank correlationMann-Whitney rank test
p-valueρp-value
Community2A<0.001−0.796NA
Bivalves2B0.123−0.418NA
Brachiopods2B<0.001−0.868NA
Smaller-sized taxaNA0.259−0.311NA
Larger-sized taxaNA0.267−0.319NA
Smaller-sized brachiopods3A0.630−0.136NA
Larger-sized brachiopods3A0.614−0.173NA
Smaller-sized bivalves3B0.8440.059NA
Larger-sized bivalves3B0.409−0.240NA
Lower half vs. upper half (community)5NANA<0.001
Lower half vs. upper half (bivalves)5NANA0.772
Lower half vs. upper half (brachiopods)5NANA<0.001

As a final step of our shell-size analyses, we constructed separate size–frequency histograms for the lower and upper parts of the precrisis interval, spurred by the observation that shell size decreased over this time interval. The boundary between the two parts was set at the level where the abundance of larger-sized species started to decrease. We applied the nonparametric Mann-Whitney rank test (Mann and Whitney 1947) to assess whether shells in the upper part are significantly smaller than in the lower part.

To compare the timing of potential changes in shell size relative with changes in biodiversity, we applied Alroy’s shareholder quorum subsampling (SQS), which calculates the number of species using a defined “coverage” or quorum (Alroy 2010). The quorum for the SQS of bivalve–brachiopod communities was fixed at 0.8. SQS-based diversity was also obtained for bivalves and brachiopods separately using a quorum of 0.6. SQS diversity was calculated for samples that had at least 35 occurrences of brachiopods and bivalves. A few successive samples in the crisis and postcrisis interval that individually did not reach this number were pooled for this purpose.

The second type of data, occurrence data from the Paleobiology Database and the literature, was used to assess the role of temperature as a stressor and potential driver of change in shell size. Using these data sources, we reconstructed the paleolatitudinal ranges of each bivalve and brachiopod species recorded from our study site during the Pliensbachian–early Toarcian interval. Assuming that the latitudinal range of a species reflects its realized thermal niche (Day et al. 2018), we expect that heat stress will most strongly affect growth in narrowly distributed, stenothermal species or in those taxa for which Portugal represents the warm edge of their latitudinal ranges. When just the modern location of an occurrence was known, paleocoordinates were obtained using A Paleolatitude Calculator for Paleoclimate Studies (van Hinsbergen et al. 2015). We focused on records from the European epicontinental seas and the western Tethys and excluded data from distant regions (e.g., Japan, South America, western Canada) to avoid mixing of regions with possibly differing latitudinal temperature gradients. We assigned each species to one of four categories—eurythermal, warm adapted, cool adapted, or stenothermal—that were then pooled into two categories—warming tolerant (eurythermal and warm adapted) and warming sensitive (cool adapted and stenothermal). The latitudinal ranges of species were thus interpreted in terms of thermal affinities relative to the geographic position of our study site, and we tested whether the taxa that are sensitive to warming are those that as a group show a significant decrease in shell size.

All analyses were conducted with the R software (v. 3.5.1; R Core Team 2017).

Results

Size data of specimens along with the position of the respective samples in the studied section are given in the Supplementary Material for all measured brachiopod and bivalve specimens and for the extended data set.

Stratigraphic Distribution of Species

Occurrence-based stratigraphic ranges provide an overview of the communities before, across, and after the T-OAE (Fig. 1). Based on faunal changes we subdivided the section into precrisis, crisis, and postcrisis intervals. The precrisis interval corresponds to the entire Polymorphum Zone. The benthic fauna is generally small and moderately diverse, with brachiopods being most abundant. The crisis interval, spanning the lower part of the Levisoni Zone, sets in at the level where faunal loss of bivalves, brachiopods, and gastropods is severe. Brachiopods in particular experienced an almost complete faunal turnover across the T-OAE. The lower part of the crisis interval is devoid of benthic shelly fossils. The upper part is characterized by high abundances of a single species, the brachiopod Soaresirhynchia bouchardi, which has been interpreted as a disaster taxon (Gahr 2002; Baeza-Carratalá 2013). A few other species associated with S. bouchardi remain very rare. The last occurrence of the athyridid species Koninckella liasina in this part of the section apparently represents a late survivor before its final demise. These faunal assemblages, although low in taxonomic richness and in evenness, indicate the early phase of the recovery of the benthic fauna. Herein, we consider the crisis interval to terminate with the last occurrence of the brachiopod S. bouchardi. The beginning of the postcrisis interval is defined by the first appearances of the so-called Spanish brachiopod fauna, which is characterized by species of the brachiopod genera Telothyris, Lobothyris, and Homoeorhynchia (Alméras 1994; Comas-Rengifo et al. 2013), in the middle Levisoni Zone. Benthic taxa in the postcrisis interval are generally less abundant but larger than in the precrisis interval. Bivalve species present in the precrisis interval commonly reappear in the postcrisis interval (Fig. 1).

Patterns and Trends in Shell Size

Analyzing for temporal trends in shell size at the community level we find that the best-fitting model is the random walk (URW). Still, the directional model (GRW) cannot be rejected because of a low (<2) ΔAICc (Table 1). When the adequacy of both models is tested (Table 2), the GRW passes all tests, while the URW fails one adequacy test for the heteroscedasticity of the residuals, making the GRW the better model of the two. This is consistent with the result of the Spearman’s rank correlation test, which indicates a significant decrease in the mean size of bivalve–brachiopod communities during the precrisis interval (Fig. 2A, Table 3). Similarly, a GRW cannot be excluded in brachiopods (for which all adequacy tests were fulfilled for both GRW and URW), while for bivalves the stasis model is clearly the best fit, despite failing one of its adequacy tests (p-value of ~0.04 when it needed to be >0.05) (Fig. 2B, Table 2). A significant decrease in the shell size of the brachiopod fauna is also confirmed by Spearman’s rank correlation test, while the trend in bivalves is not significant (Table 3).

Figure 2.

Trends in shell size in the precrisis interval. Shell size is expressed as the mean of the log2 geometric mean of shell length and height in bivalves and shell width and length in brachiopods. A, Mean shell size of the whole bivalve–brachiopod community per sample. B, Mean shell size per sample shown separately for bivalves and for brachiopods. Trend lines are based on weighted LOESS smoothing. The onset of the T-OAE (shaded area) is marked by the vertical dashed line on the right. The vertical dashed line on the left marks the Pliensbachian/Toarcian boundary.

Figure 2.

Trends in shell size in the precrisis interval. Shell size is expressed as the mean of the log2 geometric mean of shell length and height in bivalves and shell width and length in brachiopods. A, Mean shell size of the whole bivalve–brachiopod community per sample. B, Mean shell size per sample shown separately for bivalves and for brachiopods. Trend lines are based on weighted LOESS smoothing. The onset of the T-OAE (shaded area) is marked by the vertical dashed line on the right. The vertical dashed line on the left marks the Pliensbachian/Toarcian boundary.

Analysis of the shell size of smaller- and larger-sized species separately, without differentiating between bivalves and brachiopods, shows that stasis is the best fit for the group of larger-sized species (Table 1). For the group of smaller-sized species, stasis and random walk receive equal support. Both models fulfill all the adequacy tests (Table 2), so neither of the two can be preferred on this basis, but a directional trend can be rejected. Also, the related statistical tests using Spearman’s rank correlation show no significant trend in either group (Table 3). When bivalves and brachiopods are considered separately (Fig. 3A,B), stasis is the best fit in both smaller- and larger-sized bivalves (Table 1). In smaller-sized and larger-sized brachiopods, both random walk and stasis can be applied, and both models pass all their adequacy tests, while a directional trend can be rejected (Table 2). Spearman’s rank correlation tests again show that a significant trend is absent in all these subgroupings (Table 3). Yet, larger-sized species become relatively less abundant over time (Fig. 4A), a pattern prominent in brachiopods (Fig. 4B) but not in bivalves (Fig. 4C). Thus, the significant decrease in shell size across the precrisis interval can primarily be related to the larger-sized brachiopod taxa becoming less abundant with time.

Figure 3.

Per sample shell size of larger- and smaller-sized species of brachiopods (A), and bivalves (B). Each symbol in a given sample represents a different species. For further explanations, see Fig. 2. The weighted LOESS smoothing was applied to the mean of the species means.

Figure 3.

Per sample shell size of larger- and smaller-sized species of brachiopods (A), and bivalves (B). Each symbol in a given sample represents a different species. For further explanations, see Fig. 2. The weighted LOESS smoothing was applied to the mean of the species means.

Figure 4.

Relative abundance of individuals of larger-sized species in each sample expressed as the percentage of all individuals for the whole bivalve–brachiopod community (A) and for brachiopods (B) and bivalves (C). For further explanations, see Fig. 2.

Figure 4.

Relative abundance of individuals of larger-sized species in each sample expressed as the percentage of all individuals for the whole bivalve–brachiopod community (A) and for brachiopods (B) and bivalves (C). For further explanations, see Fig. 2.

Comparing the size–frequency distributions of the lower and the upper part of the precrisis interval reveals the size classes in which the overall reduction in shell size occurs (Fig. 5). The histograms are right-skewed, with the most common size classes ranging from 3.0 to 7.0 mm. The size classes larger than 21 mm are lost up-section. In brachiopods, larger-sized shells (7.0 to 19 mm) become less common, and shells larger than 19 mm disappear altogether. In bivalves, the few specimens larger than 21 mm disappear, but the other size classes are equally well represented throughout. Statistical comparison of shell sizes in the lower part with those in the upper part of the precrisis interval confirms that brachiopods are smaller in the upper part, whereas no significant difference is evident in bivalves (Table 3).

Figure 5.

Size–frequency distribution histograms (expressed as percentage) for the lower (A) and upper parts (B) of the precrisis interval. The proportion of bivalves and brachiopods are shown separately as stacked histograms. N = number of measured bivalve and brachiopod specimens.

Figure 5.

Size–frequency distribution histograms (expressed as percentage) for the lower (A) and upper parts (B) of the precrisis interval. The proportion of bivalves and brachiopods are shown separately as stacked histograms. N = number of measured bivalve and brachiopod specimens.

The size patterns of those few species/genera that are quantitatively best represented in our data are illustrated in Figure 6. In bivalves, H. spinosa shows a fairly uniform trend line of mean and maximum size, with somewhat lower size values in the uppermost part, but this may be caused by low numbers of specimens in some samples. Little net change is also evident in the brachiopods K. liasina and Nannirhynchia pygmaea, with the latter having a slight increase in maximum size across the studied interval. The shells of the terebratulid Zeilleria culeiformis tend to get smaller in the lower part of the section but are again close to their initial sizes at the end of the precrisis interval. The mean size of specimens of Liospiriferina becomes smaller up-section, but maximum size first increases until the trend is reversed before the taxon disappears from our samples.

Figure 6.

Trends in shell size in the precrisis interval for selected species. Shell size is expressed as maximum and mean log2 geometric mean (as defined in Fig. 2). For the two species Koninckella liasina and Nannirhynchia pygmaea, the maximum shell length from García Joral et al. (2018) is plotted as a comparison. Trend lines fit to our data are based on weighted LOESS smoothing. N = number of measured bivalve and brachiopod specimens.

Figure 6.

Trends in shell size in the precrisis interval for selected species. Shell size is expressed as maximum and mean log2 geometric mean (as defined in Fig. 2). For the two species Koninckella liasina and Nannirhynchia pygmaea, the maximum shell length from García Joral et al. (2018) is plotted as a comparison. Trend lines fit to our data are based on weighted LOESS smoothing. N = number of measured bivalve and brachiopod specimens.

Diversity

Bivalve–brachiopod communities are moderately diverse throughout the precrisis interval (Fig. 7A). Biodiversity fluctuated in the lower part of the interval and reached its lowest values approximately 1.3 m below the onset of the crisis interval. This diversity minimum interval corresponds to samples strongly dominated by N. pygmaea. Comparing the SQS-based diversity of bivalves and brachiopods (Fig. 7B,C), no distinct trend is observed in bivalves, while brachiopods experience a decline in the upper half of the precrisis interval. Precrisis diversity values are again reached after the crisis in the uppermost 1.8 m of the studied section.

Figure 7.

Standardized species richness of faunal samples using the shareholder quorum subsampling (SQS) of Alroy (2010). Time series of the SQS metric are shown for the whole bivalve–brachiopod community (A) and separately for bivalves (B) and brachiopods (C). For further explanations, see Fig. 1.

Figure 7.

Standardized species richness of faunal samples using the shareholder quorum subsampling (SQS) of Alroy (2010). Time series of the SQS metric are shown for the whole bivalve–brachiopod community (A) and separately for bivalves (B) and brachiopods (C). For further explanations, see Fig. 1.

Paleolatitudinal Ranges and Thermal Niches

When the thermal affinities of bivalve and brachiopod species are compared (Fig. 8), bivalves present a larger variety of thermal affinities than brachiopods (all four categories are represented) and can be grouped into seven warming-tolerant and nine warming-sensitive species. Brachiopods were assigned to just two categories, being either eurythermal (six species) or stenothermal (three species). If heat stress were the cause of the decline in the proportion of larger-sized brachiopods, the species assigned to this group should be sensitive to temperature rise. However, only one of the four larger-sized brachiopod species is categorized as being sensitive to warming. These results do not corroborate warming as the main driver toward smaller community-level shell size.

Figure 8.

Raw paleolatitudinal distribution of each species recorded from the precrisis interval. Data represent the Pliensbachian–early Toarcian time interval, apart from the bivalve species Homomya gibbosa, for which all Jurassic occurrences were used to circumvent lack of data. The larger-sized species in both bivalves and brachiopods are marked with an asterisk (*). Assuming that latitudinal ranges reflect thermal affinities, species were grouped into four categories relative to the geographic position of our study site. The taxon Liospiriferina spp. includes all species belonging to this genus as recorded in the studied section. For the purpose of this analysis, species identified with reservation (i.e., with the identifier “cf.”) in the Paleobiology Database and the literature were considered as true representatives of the respective species. The vertical dashed line marks the paleolatitude of the Fonte Coberta section.

Figure 8.

Raw paleolatitudinal distribution of each species recorded from the precrisis interval. Data represent the Pliensbachian–early Toarcian time interval, apart from the bivalve species Homomya gibbosa, for which all Jurassic occurrences were used to circumvent lack of data. The larger-sized species in both bivalves and brachiopods are marked with an asterisk (*). Assuming that latitudinal ranges reflect thermal affinities, species were grouped into four categories relative to the geographic position of our study site. The taxon Liospiriferina spp. includes all species belonging to this genus as recorded in the studied section. For the purpose of this analysis, species identified with reservation (i.e., with the identifier “cf.”) in the Paleobiology Database and the literature were considered as true representatives of the respective species. The vertical dashed line marks the paleolatitude of the Fonte Coberta section.

Discussion

Selective Faunal Response of Brachiopods

We find that brachiopods are more strongly affected before and at the T-OAE than bivalves. Most species belonging to the brachiopod fauna of the Polymorphum Zone went extinct in their distribution area, that is, the northwest European basins and the Mediterranean Province (García Joral et al. 2011, 2018; Baeza-Carratalá 2013; Comas-Rengifo et al. 2013, 2015). This selectivity against brachiopods is evident in a drop in brachiopod diversity (Fig. 7) followed by an almost complete faunal turnover across the T-OAE (Fig. 1), as well as significant reductions in mean shell size of brachiopods. Thus, our first hypothesis—that body size declined before the T-OAE—is supported by our results, albeit only selectively for brachiopods. Abundance declines and extirpations of larger-sized taxa before the event are the underlying mechanism of this size reduction in the brachiopod fauna. By contrast, single species and subgroupings of species, such as the group of larger-sized brachiopod species, do not show significant declines in shell size. Thus, our second hypothesis—that larger-sized species were more affected than smaller-sized ones—is only supported in the sense that larger-sized brachiopod taxa became less common over time. It is also apparent that shell size decreased fairly continuously during the precrisis interval (Fig. 2), while diversity only declined in the uppermost part of the precrisis interval (Fig. 7).

Assuming constant sedimentation rates and a duration of the precrisis interval of ~900 kyr, the offset between the abundance decline of larger-sized brachiopod species between 3 and 4 m above the Pliensbachian/Toarcian boundary (Fig. 4) and the drop in species richness at 6 m above the boundary (Fig. 7) would correspond to ~300 kyr. However, this is likely to be an overestimation, because the reduced thickness of the Mirabile Subzone at Fonte Coberta and comparison with the Polymorphum Zone at Peniche (Rita et al. 2016) suggest that the lower part of the Polymorphum Zone is condensed. Notwithstanding this, a temporal decoupling remains, supporting the third hypothesis—that reductions in body size occurred earlier than diversity loss. So-called early warning signals have been hypothesized to predict sudden changes in species composition and ecological structure in present-day ecosystems (Biggs et al. 2012; Gsell et al. 2016). In analogy, a decline in body size may serve as a precursor of imminent turnover at the community level at geologic timescales.

A recent comparative study of brachiopod shell size during the late Pliensbachian–early Toarcian interval among basins surrounding the Iberian Massif found decreases in shell size as the depositional environments generally became deeper, more turbid, and less oxygenated from the Iberian basins in Spain to the Lusitanian Basin in Portugal (García Joral et al. 2018). Rather than being attributable to miniaturization of species, this trend was caused by a change in taxonomic composition, with small-sized species being more common in the Lusitanian Basin. García Joral et al. (2018) also investigated within-basin and within-section size changes through time. For the Lusitanian Basin as a whole, they reported smaller maximum and mean shell sizes of N. pygmaea in the Mirabile Subzone than in the Semicelatum Subzone (García Joral et al. 2018: Figs. 3, 4), and thus inferred an increase in shell size over time. Similarly, when plotting maximum shell size of N. pygmaea and of K. liasina in several beds of the Fonte Coberta section (García Joral et al. 2018: Fig. 5), that is, the same section studied here, they found the largest specimens of these small-sized species in the upper part of the Semicelatum Subzone. We contrast the Fonte Coberta size data of García Joral et al. (2018) with our own data for these two brachiopod species (Fig. 6). The maximum sizes of shells measured by García Joral et al. (2018) are distinctly larger than those of our specimens. This is likely caused by differing sampling procedures, as we confined our analyses to specimens from the standardized sample size of 15 kg of rock described earlier. Even though our time series of measurements are more complete, the slight increase in the maximum size of N. pygmaea is inferred in both studies. In contrast to our study, García Joral et al. (2018) also report an increase in maximum shell length for K. liasina, but this interpretation relies on just one data point and should be considered with caution (see Fig. 6). In any case, our main conclusion—that size decrease in the brachiopod fauna was caused by larger-sized species becoming less common or disappearing from the record entirely rather than individual species becoming smaller—is not affected by the interpretations of García Joral et al. (2018), who did not examine this aspect.

Causes of Biotic Responses

The early Toarcian mass extinction was global in extent and involved the worldwide demise of the brachiopod orders Spiriferinida and Athyridida (Vörös 2002; García Joral et al. 2011; Vörös et al. 2016). The ultimate cause(s) of this event must therefore be global in nature, while the proximate killing mechanisms might still vary depending on the environmental context. Although the factors causing faunal turnover and extinction need not necessarily be the same factors that caused previous change in body size, the shared selectivity against brachiopods makes a common-cause scenario plausible.

A particular threat for marine organisms today is the stress that arises from the coupling of global warming, ocean acidification, and ocean deoxygenation, collectively termed the “deadly trio” (Bijma et al. 2013). Episodes of mass extinction in the geologic past commonly involve these three stressors (Jenkyns 2010), and they also have been considered as a cause for the early Toarcian extinction event (see “Introduction”). Changes in nutrient cycling and productivity may be added as a fourth warming-related stressor, thus adding up to a “deadly quartet.” Because sustained warming promotes ocean stratification, nutrients may be transferred to and trapped in deeper waters, thus reducing global productivity (Moore et al. 2018). Alternatively, an accelerated hydrological cycle could increase regional continental weathering rates and nutrient input into the sea, resulting in eutrophication, expansion of oxygen minimum zones onto the shelf, and toxic algal blooms. Below, we evaluate each of these four main stressors with regard to our own findings.

Identifying the abiotic causes of faunal change in the geologic past often involves facies analysis and evaluation of geochemical proxy data. Integrating multiple biological disciplines by using fossil ecological data, modern ecological data, and physiological data is a relatively new approach in Earth system science to improve understanding of past mass extinctions and current biosphere change and to predict ecological consequences of climate change in the near future (e.g., Knoll et al. 2007; Calosi et al. 2019). Biological concepts such as the concept of oxygen- and capacity-limited thermal tolerance can successfully bridge between physiology and ecology (Pörtner et al. 2017), albeit quantitative links between physiological processes and ecosystem-level processes are still limited. It is less clear to what extent results from physiological experiments can be applied to explain paleoecological patterns. Although physiological experiments can investigate the short-term metabolic costs of environmental stress and thus go beyond merely reporting survival or failure, temporal scaling represents an important difference. While the time-averaged nature of paleoecological data hampers forecasting the near future, physiological experiments with animals last much shorter than even geologically abrupt events of warming or the time necessary for evolutionary adaptation. Physiological limits will still be important at longer timescales, but other processes not captured by experiments will come into play (Peck et al. 2009). The following discussion of potential causes of the identified paleoecological patterns includes physiological facets and assumes that the physiological tolerances of species explored in organismic-scale lab experiments can provide information about the sensitivity of taxa to climate-related stressors in the geologic past.

Hypoxia

There is no evidence of low-oxygen conditions in the benthic environments of the Rabaçal section. Total organic carbon levels are generally low (Duarte et al. 2005), black shales are lacking (see also Duarte 1997; Duarte et al. 2007), and bioturbation during both the precrisis and the crisis intervals indicates oxygenated bottom conditions throughout. In particular, the interval corresponding to the T-OAE elsewhere is intensely bioturbated, as indicated by pervasive networks of Thalassinoides burrows, which were produced by crustaceans (see also Miguez-Salas et al. 2017; Rodríguez-Tovar et al. 2017). Crustaceans are sensitive to hypoxia (Vaquer-Sunyer and Duarte 2008), and their apparent former abundance is therefore additional evidence against low-oxygen conditions on and within the seafloor. Limited oxygen availability may have played a role in the deeper parts of the Lusitanian Basin, as represented by the sedimentary record at Peniche (e.g., Hesselbo et al. 2007), but apparently not in the mid-ramp setting investigated here.

Ocean Acidification

Bivalves appear to be generally negatively affected by acidification, although some species exhibit neutral or even positive effects (Parker et al. 2013). In the early Toarcian fauna studied herein, bivalves as a group are less affected by size reduction and diversity decline than brachiopods. Their low-magnesium calcite shell makes rhynchonelliform brachiopods more resistant to acidification than organisms with aragonitic or high-magnesium calcite shells (Ries et al. 2009). Because organisms with aragonitic shells are represented in the precrisis interval (e.g., nuculoid bivalves), we expect that we would have observed any selective pattern against aragonitic species, if present. In previous studies, rhynchonelliform brachiopods had been considered to have minimal physiological buffering capacities of the calcification process, making them sensitive to CO2 stress (e.g., Knoll et al. 2007). Conversely, recent work on living brachiopods and on material from museum collections from the last 120 years shows little effect of lowered pH conditions on shell growth (Cross et al. 2015, 2016, 2018). Selectivity against brachiopods in the reduction of shell size and loss of species at Rabaçal do not conform with the inferred physiological pH tolerance of bivalves and brachiopods; thus, any clear support for acidification from faunal patterns is missing.

Productivity Decline

Unlike mollusks, the brachiopods, and rhynchonellids in particular, can fare relatively well in low-nutrient environments, as they are capable of feeding on both dissolved and particulate organic matter (Steele-Petrović 1976) and have very low metabolic rates (Peck and Harper 2010). Carbon isotope and nannofossil abundance data indicate mesotrophic to eutrophic, albeit unstable, conditions throughout the early Toarcian in the Lusitanian Basin (Mattioli et al. 2009; Ferreira et al. 2015), while high productivity and upwelling are suggested for the anoxic settings in northern and central Europe (e.g., Jenkyns 1988, 2010). Adding our finding that brachiopods are preferentially affected by environmental perturbations at Rabaçal, a low-productivity scenario is unlikely.

Warming

The temperature increase at the T-OAE was not globally homogeneous. Estimates range from +2°C to +3.5°C in subtropical areas (Dera and Donnadieu 2012), and between +6°C and +8°C at higher latitudes (Dera et al. 2009; Gómez and Goy 2011). Reliable oxygen isotope data from the Polymorphum Zone of the studied section are lacking due to the poor preservation of brachiopod shells. However, oxygen isotopes from brachiopods elsewhere in the Lusitanian Basin reveal a first short-lived warming episode in the lowermost Polymorphum Zone, followed by a cooling phase, before a marked warming phase started in the mid-Polymorphum Zone that culminated during the T-OAE in the lowermost Levisoni Zone (Suan et al. 2010). Warming has been suggested as the main cause of brachiopod faunal turnover and diversity decrease in the oxygenated environments of western Europe (e.g., Vörös 2002; García Joral et al. 2011, 2018; Gómez and Goy 2011; Miguez-Salas et al. 2017). Indeed, the loss of fauna is recorded shortly after the onset of the T-OAE at Rabaçal, suggesting a role of warming in brachiopod extinctions. Evidence that warming was also a main driver of the decrease in mean community-level shell size is equivocal. The thermal affinities inferred for brachiopods from their latitudinal distributions do not hint at temperature stress as the main cause of the size patterns recorded from Rabaçal. On the other hand, modern bivalves are among the organisms with highest upper thermal limits and are therefore one of the most thermally tolerant groups (Song et al. 2014), which would match the pattern of stasis in bivalves. Predictions using physiological responses of brachiopods are made difficult by the scarcity of data regarding their thermal tolerances. An exception is the Antarctic rhynchonellid Liothyrella uva, which exhibits lower thermal tolerance to rapid warming than two simultaneously studied Antarctic bivalve species (Clark et al. 2017), but results from a single species need not be universally valid for brachiopods as a whole. Clearly, more experimental work on the physiological tolerances of modern brachiopods in the face of multiple stressors is needed to improve our understanding of selective biotic responses to environmental stress.

Conclusions

Long-term decrease on the order of a few hundred thousand years in the mean body size of early Toarcian marine invertebrate communities from the Lusitanian Basin, well before the main phase of local extirpations and biodiversity decrease occurred, suggests that reductions in body size may be one of the first ecological responses to abiotic stressors. Environmental stress acted selectively against specific size classes and taxonomic groups, that is, larger-sized brachiopod species, which became less abundant over time. Future research has to show to which degree the patterns of decreasing shell sizes observed at Rabaçal can be transferred to other ecosystems and thus are general predictors of forthcoming climate-related community turnover in ancient ecosystems.

In contrast to many other regions, geologic and paleontological evidence from Rabaçal is against a role of hypoxic conditions as a driver of both body-size decrease and faunal loss. Similarly, a reduced oceanic pH alone seems incompatible with the observed faunal trends, although a definitive conclusion is hampered by the scarcity of comparative experimental data from bivalves and brachiopods under elevated CO2 conditions. Heat stress is a plausible main cause of diversity decline and elevated early Toarcian extinction intensity in aerated environments. The rising paleo-temperatures from the mid-Polymorphum Zone to the earliest Levisoni Zone, as inferred from oxygen isotopes, are also compatible with heat stress as a driver of the precrisis decrease in shell size. However, biotic evidence is equivocal at present, because the selective size decline of early Toarcian brachiopods is not matched by their thermal affinities as inferred from their latitudinal distributions, and the physiological response of brachiopods to warming has hardly been studied experimentally. Rising temperature and increasing pCO2 may have operated additively or synergistically in our study system with different combined effects than when acting individually. Examination of more specific physiology-based hypotheses on biotic change, combined with analysis of multiple proxy data on changing environmental conditions, will likely improve our understanding of the interactions between abiotic stressors and biotic responses in the studied early Toarcian succession and beyond. Yet, identifying the exact role of individual players in the “deadly quartet” of warming, ocean acidification, ocean deoxygenation, and productivity change will remain a big challenge in the analysis of ancient ecosystems.

Acknowledgments

This work was funded by the Deutsche Forschungsgemeinschaft grant DFG AB 09/10-1 and is part of the Research Unit TERSANE (FOR 2332: Temperature-related Stressors as a Unifying Principle in Ancient Extinctions). This research is also a contribution to the IGCP-655 (IUGS-UNESCO: Toarcian Oceanic Anoxic Event: Impact on Marine Carbon Cycle and Ecosystems). We are grateful to W. Foster (Museum für Naturkunde) and S. Schneider (Cambridge Arctic Shelf Programme) for commenting on earlier versions of the article. T. Klein and F. Lucassen (Center for Marine Environmental Sciences, University of Bremen) and S. Schneider (Cambridge Arctic Shelf Programme) are warmly thanked for the joint fieldwork in Portugal. Comments from L. Harper (University of Cambridge) and an anonymous reviewer substantially improved an earlier version of this article.

Literature Cited

Aberhan
,
M.
, and
T. K.
Baumiller
.
2003
.
Selective extinction among Early Jurassic bivalves: a consequence of anoxia
.
Geology
 
31
:
1077
1080
.
Alméras
,
Y.
1994
.
Le genre Soaresirhynchia nov. (Brachiopoda, Rhynchonellacea, Wellerellidae) dans le Toarcien du sous-bassin nord-lusitanien (Portugal)
.
Documents Laboratoire Géologie Lyon
 
130
:
1
136
.
Alroy
,
J.
2010
.
Fair sampling of taxonomic richness and unbiased estimation of origination and extinction rates
.
Paleontological Society Papers
 
16
:
55
80
.
Baeza-Carratalá
,
J. F.
2013
.
Diversity patterns of Early Jurassic brachiopod assemblages from the westernmost Tethys (Eastern Subbetic)
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
381–382
:
76
91
.
Biggs
,
R.
,
T.
Blenckner
,
C.
Folke
,
L.
Gordon
,
A.
Norström
,
M.
Nyström
, and
G. D.
Peterson
.
2012
.
Regime shifts
. Pp.
609
617
in
A.
Hastings
and
L.
Gross
, eds.
Encyclopedia of theoretical ecology
 .
University of California Press
,
Berkeley, Calif
.
Bijma
,
J.
,
H.-O.
Pörtner
,
C.
Yesson
, and
A. D.
Rogers
.
2013
.
Climate change and the oceans—What does the future hold?
Marine Pollution Bulletin
 
74
:
495
505
.
Boulila
,
S.
,
B.
Galbrun
,
E.
Huret
,
L. A.
Hinnov
,
I.
Rouget
,
S.
Gardin
, and
A.
Bartolini
.
2014
.
Astronomical calibration of the Toarcian Stage: implications for sequence stratigraphy and duration of the early Toarcian OAE
.
Earth and Planetary Science Letters
 
386
:
98
111
.
Box
,
G. E. P.
, and
D. A.
Pierce
.
1970
.
Distribution of residual correlations in autoregressive-integrated moving average time series models
.
Journal of the American Statistical Association
 
65
:
1509
1526
.
Breitburg
,
D.
,
L. A.
Levine
,
A.
Oschlies
,
M.
Grégoire
,
F. P.
Chavez
,
D. J.
Conley
,
V.
Garçon
,
D.
Gilbert
,
D.
Gutiérrez
,
K.
Isensee
,
G. S.
Jacinto
,
K. E.
Limburg
,
I.
Montes
,
S. W. A.
Naqvi
,
G. C.
Pitcher
,
N. N.
Rabalais
,
M. R.
Roman
,
K. A.
Rose
,
B. A.
Seibel
,
M.
Telszewski
,
M.
Yasuhara
, and
J.
Zhang
.
2018
.
Declining oxygen in the global ocean and coastal waters
.
Science
 
359
:
eaam7240
.
Burnham
,
K. P.
, and
D. R.
Anderson
.
2003
.
Model selection and multimodel inference: a practical information-theoretic approach
 .
Springer
,
Berlin
.
Byrne
,
M.
, and
R.
Przeslawski
.
2013
.
Multistressor impacts of warming and acidification of the ocean on marine invertebrates’ life histories
.
Integrative and Comparative Biology
 
53
:
582
596
.
Calosi
,
P.
,
H. M.
Putnam
,
R. J.
Twitchett
, and
F.
Vermandele
.
2019
.
Marine metazoan modern mass extinction: improving predictions by integrating fossil, modern, and physiological data
.
Annual Review of Marine Science
 
11
:
20.1
20.22
.
Clark
,
M. S.
,
U.
Sommer
,
J. K.
Sihra
,
M. A. S.
Thorne
,
S. A.
Morley
,
M.
King
,
M. R.
Viant
, and
L. S.
Peck
.
2017
.
Biodiversity in marine invertebrate responses to acute warming revealed by a comparative multi-omics approach
.
Global Change Biology
 
23
:
318
330
.
Comas-Rengifo
,
M. J.
,
L. V.
Duarte
,
F.
García Joral
, and
A.
Goy
.
2013
.
The brachiopod record in the Lower Toarcian (Jurassic) of the Rabaçal–Condeixa region (Portugal): stratigraphic distribution and palaeobiogeography
.
Comunicações Geológicas
 
100
:
37
42
.
Comas-Rengifo
,
M. J.
,
L. V.
Duarte
,
F. F.
Félix
,
F.
García Joral
,
A.
Goy
, and
R. B.
Rocha
.
2015
.
Latest Pliensbachian–Early Toarcian brachiopod assemblages from the Peniche section (Portugal) and their correlation
.
Episodes
 
38
:
2
7
.
Correia
,
V. F.
,
J. B.
Riding
,
L. V.
Duarte
,
P.
Fernandes
, and
Z.
Pereira
.
2018
.
The Early Jurassic palynostratigraphy of the Lusitanian Basin, western Portugal
.
Geobios
 
51
:
537
557
Cross
,
E. L.
,
L. S.
Peck
, and
E. M.
Harper
.
2015
.
Ocean acidification does not impact shell growth or repair of the Antarctic brachiopod Liothyrella uva (Broderip, 1833)
.
Journal of Experimental Marine Biology and Ecology
 
462
:
29
35
.
Cross
,
E. L.
,
L. S.
Peck
,
M. D.
Lamare
, and
E. M.
Harper
.
2016
.
No ocean acidification effects on shell growth and repair in the New Zealand brachiopod Calloria inconspicua (Sowerby, 1846)
.
ICES Journal of Marine Science
 
73
:
920
926
.
Cross
,
E. L.
,
E. M.
Harper
, and
L. S.
Peck
.
2018
.
A 120-year record of resilience to environmental change in brachiopods
.
Global Change Biology
 
24
:
2262
2271
.
Danise
,
S.
,
R. J.
Twitchett
, and
C. T. S.
Little
.
2015
.
Environmental controls on Jurassic marine ecosystems during global warming
.
Geology
 
43
:
263
266
.
Daufresne
,
M.
,
K.
Lengfellner
, and
U.
Sommer
.
2009
.
Global warming benefits the small in aquatic ecosystems
.
Proceedings of the National Academy of Sciences USA
 
106
:
12788
12793
.
Day
,
P. B.
,
R. D.
Stuart-Smith
,
G. J.
Edgar
, and
A. E.
Bates
.
2018
.
Species’ thermal ranges predict changes in reef fish community structure during 8 years of extreme temperature variation
.
Diversity and Distributions
 
24
:
1036
1046
.
Dera
,
G.
, and
Y.
Donnadieu
.
2012
.
Modeling evidence for global warming, Arctic seawater freshening, and sluggish oceanic circulation during the Early Toarcian anoxic event
.
Paleoceanography
 
27
:
PA2211
.
Dera
,
G.
,
E.
Pucéat
,
P.
Pellenard
,
P.
Neige
,
D.
Delsate
,
M. M.
Joachimski
,
L.
Reisberg
, and
M.
Martinez
.
2009
.
Water mass exchange and variations in seawater temperature in the NW Tethys during the Early Jurassic: evidence from neodymium and oxygen isotopes of fish teeth and belemnites
.
Earth and Planetary Science Letters
 
286
:
198
207
.
Diaz
,
R. J.
, and
R.
Rosenberg
.
2008
.
Spreading dead zones and consequences for marine ecosystems
.
Science
 
321
:
926
929
.
Duarte
,
L. V.
1997
.
Facies analysis and sequential evolution of the Toarcian–Lower Aalenian series in the Lusitanian Basin (Portugal)
.
Communicações do Instituto Geológico e Mineiro
 
83
:
65
94
.
Duarte
,
L. V.
2007
.
Lithostratigraphy, sequence stratigraphy and depositional setting of the Pliensbachian and Toarcian series in the Lusitanian Basin, Portugal
. Pp.
17
23
in
R. B.
Rocha
, ed.
The Peniche section (Portugal)
 .
Contributions to the definition of the Toarcian GSSP. International Subcommission on Jurassic Stratigraphy
,
Lisbon
.
Duarte
,
L. V.
, and
A. F.
Soares
.
2002
.
Litostratigrafia das séries margo-calcárias do Jurássico inferior da Bacia Lusitânica (Portugal)
.
Comunicacões do Instituto Geoloógico e Mineiro
 
89
:
135
154
.
Duarte
,
L. V.
,
R.
Rodrigues
,
L. C.
Oliveira
, and
F.
Silva
.
2005
.
Avaliação preliminar das variações do carbono orgânico total nos sedimentos margosos do Jurássico inferior da Bacia Lusitânica (Portugal)
. In
XIV Semana de Geoquímica and VIII Congresso de Geoquímica dos Países de Lingua Portuguesa
 
1
:
39
43
.
Aveiro, Portugal.
Duarte
,
L. V.
,
L. C.
Oliveira
, and
R.
Rodrigues
.
2007
.
Carbon isotopes as a sequence stratigraphic tool: examples from the Lower to Middle Toarcian marly limestones of Portugal
.
Boletín Geológico y Minero
 
118
:
3
18
.
Fabry
,
V. J.
,
B. A.
Seibel
,
R. A.
Feely
, and
J. C.
Orr
.
2008
.
Impacts of ocean acidification on marine fauna and ecosystem Processes
.
ICES Journal of Marine Science
 
65
:
414
432
.
Ferreira
,
J.
,
E.
Mattioli
,
B.
Pittet
,
M.
Cachão
, and
J. E.
Spanbenberg
.
2015
.
Palaeoecological insights on Toarcian and lower Aalenian calcareous nannofossils from the Lusitanian Basin (Portugal)
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
436
:
246
262
.
Gahr
,
M. E.
2002
.
Palökologie des Makrobenthos aus dem Unter-Toarc SW-Europas
.
Beringeria
 
31
:
3
204
.
García Joral
,
F.
,
J. J.
Gómez
, and
A.
Goy
.
2011
.
Mass extinction and recovery of the Early Toarcian (Early Jurassic) brachiopods linked to climate change in Northern and Central Spain
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
302
:
367
380
.
García Joral
,
F.
,
J. F.
Baeza-Carratalá
, and
A.
Goy
.
2018
.
Changes in brachiopod body size prior to the Early Toarcian (Jurassic) mass extinction
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
506
:
242
249
.
Gardner
,
J. L.
,
A.
Peters
,
M. R.
Kearney
,
L.
Joseph
, and
R.
Heinsohn
.
2011
.
Declining body size: a third universal response to warming?
Trends in Ecology and Evolution
 
26
:
285
291
.
Gazeau
,
F.
,
L. M.
Parker
,
S.
Comeau
,
J.-P.
Gattuso
,
W. A.
O’Connor
,
S.
Martin
,
H.- O.
Pörtner
, and
P. M.
Ross
.
2013
.
Impacts of ocean acidification on marine shelled molluscs
.
Marine Biology
 
160
:
2207
2245
.
Genner
,
M. J.
,
D. W.
Sims
,
A. J.
Southward
,
G. C.
Budd
,
P.
Masterson
,
M.
McHugh
,
P.
Rendle
,
E. J.
Southall
,
V. J.
Wearmouth
, and
S. J.
Hawkins
.
2009
.
Body size-dependent responses of a marine fish assemblage to climate change and fishing over a century-long scale
.
Global Change Biology
 
16
:
517
527
.
Gienapp
,
P.
,
C.
Teplitsky
,
J. S.
Alho
,
J. A.
Mills
, and
J.
Merilä
.
2008
.
Climate change and evolution: disentangling environmental and genetic responses
.
Molecular Ecology
 
17
:
167
178
.
Gómez
,
J. J.
, and
A.
Goy
.
2011
.
Warming-driven mass extinction in the Early Toarcian (Early Jurassic) of northern and central Spain. Correlation with other time-equivalent European sections
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
306
:
176
195
.
Gsell
,
A. S.
,
U.
Scharfenberger
,
D.
Özkundakci
,
A.
Walters
,
L.-A.
Hansson
,
A. B. G.
Janssen
,
P.
Nõges
,
P. C.
Reid
,
D. E.
Schindler
,
E.
Van Donk
,
V.
Dakos
, and
R.
Adrian
.
2016
.
Evaluating early-warning indicators of critical transitions in natural aquatic ecosystems
.
Proceedings of the National Academy of Sciences USA
 
113
:
E8089
E8095
.
Hallam
,
A.
, and
P. B.
Wignall
.
1997
.
Mass extinctions and their aftermath
 .
Oxford University Press
,
Oxford
.
He
,
W.
,
G. R.
Shi
,
Q.
Feng
,
M. J.
Campi
,
S.
Gu
,
J.
Bu
,
Y.
Peng
, and
Y.
Meng
.
2007
.
Brachiopod miniaturization and its possible causes during the Permian–Triassic crisis in deep water environments, South China
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
252
:
145
163
.
He
,
W.
,
G. R.
Shi
,
Y.
Xiao
,
K.
Zhang
,
T.
Yang
,
H.
Wu
,
Y.
Zhang
,
B.
Chen
,
M.
Yue
,
J.
Shen
,
Y.
Wang
,
H.
Yang
, and
S.
Wu
.
2017
.
Body-size changes of latest Permian brachiopods in varied palaeogeographic settings in South China and implications for controls on animal miniaturization in a highly stressed marine ecosystem
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
486
:
33
45
.
He
,
W.-H.
,
R. J.
Twitchett
,
Y.
Zhang
,
G. R.
Shi
,
Q.-L.
Feng
,
J.-X.
Yu
,
S.-B.
Wu
, and
X.-F.
Peng
.
2010
.
Controls on body size during the Late Permian mass extinction event
.
Geobiology
 
8
:
391
402
.
Hesselbo
,
S. P.
,
D. R.
Gröcke
,
H. C.
Jenkyns
,
C. J.
Bjerrum
,
P.
Farrimond
,
H. S.
Morgans Bell
, and
O. R.
Green
.
2000
.
Massive dissociation of gas hydrate during a Jurassic oceanic anoxic event
.
Nature
 
406
:
392
395
.
Hesselbo
,
S. P.
,
H. C.
Jenkyns
,
L. V.
Duarte
, and
L. C. V.
Oliveira
.
2007
.
Carbon-isotope record of the Early Jurassic (Toarcian) Oceanic Anoxic Event from fossil wood and marine carbonate (Lusitanian Basin, Portugal)
.
Earth and Planetary Science Letters
 
253
:
455
470
.
Huang
,
B.
,
D. A. T.
Harper
,
R.
Zhan
, and
J.
Rong
.
2010
.
Can the Lilliput effect be detected in the brachiopod faunas of South China following the terminal Ordovician mass extinction?
Palaeogeography, Palaeoclimatology, Palaeoecology
 
285
:
277
286
.
Hunt
,
G.
2006
.
Fitting and comparing models of phyletic evolution: random walks and beyond
.
Paleobiology
 
32
:
578
601
.
Hunt
,
G.
2015
.
paleoTS: analyze paleontological time-series. R package, Version 0.5-1
. https://CRAN.R-project.org/package=paleoTS, accessed 22 October 2018.
Jenkyns
,
H. C.
1988
.
The early Toarcian (Jurassic) anoxic event: stratigraphic, sedimentary, and geochemical evidence
.
American Journal of Science
 
288
:
101
151
.
Jenkyns
,
H. C.
2010
.
Geochemistry of oceanic anoxic events
.
Geochemistry, Geophysics, Geosystems
 
11
:
Q03004
.
Kelly
,
M. W.
, and
G. E.
Hofmann
.
2013
.
Adaptation and the physiology of ocean acidification
.
Functional Ecology
 
27
:
980
990
.
Kiessling
,
W.
,
M.
Schobben
,
A.
Ghaderi
,
V.
Hairapetian
,
L.
Leda
, and
D.
Korn
.
2018
.
Pre–mass extinction decline of latest Permian ammonoids
.
Geology
 
46
:
283
286
.
Knoll
,
A. H
,
R. K.
Bambach
,
J. L.
Payne
,
S.
Pruss
, and
W. W.
Fischer
.
2007
.
Paleophysiology and end-Permian mass extinction
.
Earth and Planetary Science Letters
 
256
:
295
313
.
Kosnik
,
M. A.
,
D.
Jablonski
,
R.
Lockwood
, and
P. M.
Novack-Gottshall
.
2006
.
Quantifying molluscan body size in evolutionary and ecological analyses: maximizing the return on data-collection efforts
.
Palaios
 
21
:
588
597
.
Lockwood
,
R.
2005
.
Body size, extinction events, and the early Cenozoic record of veneroid bivalves: a new role for recoveries?
Paleobiology
 
31
:
578
590
.
Mann
,
H. B.
, and
D. R.
Whitney
.
1947
.
On a test of whether one of two random variables is stochastically larger than the other
.
Annals of Mathematical Statistics
 
18
:
50
60
.
Martindale
,
R. C.
, and
M.
Aberhan
.
2017
.
Response of macrobenthic communities to the Toarcian Oceanic Anoxic Event in northeastern Panthalassa (Ya Ha Tinda, Alberta, Canada)
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
478
:
103
120
.
Mattioli
,
E.
,
B.
Pittet
,
L.
Petitpierre
, and
S.
Mailliot
.
2009
.
Dramatic decrease of pelagic carbonate production by nannoplankton across the Early Toarcian anoxic event (T-OAE)
.
Global and Planetary Change
 
65
:
134
145
.
Miguez-Salas
,
O.
,
F. J.
Rodríguez-Tovar
, and
L. V.
Duarte
.
2017
.
Selective incidence of the toarcian oceanic anoxic event on macroinvertebrate marine communities: a case from the Lusitanian basin, Portugal
.
Lethaia
 
50
:
548
560
.
Moore
,
J. K.
,
W.
Fu
,
F.
Primeau
,
G. L.
Britten
,
K.
Lindsay
,
M.
Long
,
S. C.
Doney
,
N.
Mahowald
,
F.
Hoffman
, and
J. T.
Randerson
.
2018
.
Sustained climate warming drives declining marine biological productivity
.
Science
 
359
:
1139
1143
.
Morten
,
S. D.
, and
R. J.
Twitchett
.
2009
.
Fluctuations in the body size of marine invertebrates through the Pliensbachian–Toarcian extinction event
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
284
:
29
38
.
Mouterde
,
R.
,
C.
Ruget
, and
F.
Moitinho de Almeida
.
1964–1965
.
Coupe du Lias au Sud de Condeixa
.
Comunicações dos Serviços Geológicas de Portugal
 
48
:
61
91
.
O’Gorman
,
E. J.
,
L.
Zhao
,
D. E.
Pichler
,
G.
Adams
,
N.
Friberg
,
B. C.
Rall
,
A.
Seeney
,
H.
Zhang
,
D. C.
Reuman
, and
G.
Woodward
.
2017
.
Unexpected changes in community size structure in a natural warming experiment
.
Nature Climate Change
 
7
:
659
666
.
Ohlberger
,
J.
2013
.
Climate warming and ectotherm body size—from individual physiology to community ecology
.
Functional Ecology
 
27
:
991
1001
.
Parker
,
L. M.
,
P. M.
Ross
,
W. H.
O’Connor
,
H.-O.
Pörtner
,
E.
Scanes
, and
J. M.
Wright
.
2013
.
Predicting the response of molluscs to the impact of ocean acidification
.
Biology
 :
651
692
.
Pálfy
,
J.
, and
P. L.
Smith
.
2000
.
Synchrony between Early Jurassic extinction, oceanic anoxic event, and the Karoo-Ferrar flood basalt volcanism
.
Geology
 
28
:
747
750
.
Peck
,
L. S.
, and
E. M.
Harper
.
2010
.
Variation in size of living articulated brachiopods with latitude and depth
.
Marine Biology
 
157
:
2205
2213
.
Peck
,
L. S.
,
M. S.
Clark
,
S. A.
Morley
,
A.
Massey
, and
H.
Rossetti
.
2009
.
Animal temperature limits and ecological relevance: effects of size, activity and rates of change
.
Functional Ecology
 
23
:
248
256
.
Pittet
,
B.
,
G.
Suan
,
F.
Lenoir
,
L.V.
Duarte
, and
E.
Mattioli
.
2014
.
Carbon isotope evidence for sedimentary discontinuities in the Lower Toarcian of the Lusitanian Basin (Portugal): sea level change at the onset of the Oceanic Anoxic Event
.
Sedimentary Geology
 
303
:
1
14
.
Pörtner
,
H.-O.
2008
.
Ecosystem effects of ocean acidification in times of ocean warming: a physiologist’s view
.
Marine Ecology Progress Series
 
373
:
203
217
.
Pörtner
,
H.-O.
, and
A. P.
Farrell
.
2008
.
Physiology and climate change
.
Science
 
322
:
690
692
.
Pörtner
,
H.-O.
,
M.
Langenbuch
, and
B.
Michaelidis
.
2005
.
Synergistic effects of temperature extremes, hypoxia, and increases in CO2 on marine animals: from Earth history to global change
.
Journal of Geophysical Research
 
110
:
C09S10
.
Pörtner
,
H.-O.
,
C.
Bock
, and
F. C.
Mark
.
2017
.
Oxygen- and capacity-limited thermal tolerance: bridging ecology and physiology
.
Journal of Experimental Biology
 
220
:
2685
2696
.
R Core Team
.
2017
.
R: a language and environment for statistical computing
 .
R Foundation for Statistical Computing
,
Vienna, Austria
.
Ries
,
J. B.
,
A. L.
Cohen
, and
D.
McCorkle
.
2009
.
Marine calcifiers exhibit mixed responses to CO2-induced ocean acidification
.
Geology
 
37
:
1131
1134
.
Rita
,
P.
,
M.
Reolid
, and
L. V.
Duarte
.
2016
.
Benthic foraminiferal assemblages record major environmental perturbations during the Late Pliensbachian–Early Toarcian interval in the Peniche GSSP, Portugal
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
454
:
267
281
.
Rodríguez-Tovar
,
F. J.
,
O.
Miguez-Salas
, and
L. V.
Duarte
.
2017
.
Toarcian Oceanic Anoxic Event induced unusual behaviour and palaeobiological changes in Thalassinoides tracemakers
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
485
:
46
56
.
Röhl
,
H.-J.
,
A.
Schmidt-Röhl
,
W.
Oschmann
,
A.
Frimmel
, and
L.
Schwark
.
2001
.
The Posidonia Shale (Lower Toarcian) of SW-Germany: an oxygen-depleted ecosystem controlled by sea level and palaeoclimate
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
165
:
27
52
.
Sell
,
B.
,
M.
Ovtcharova
,
J.
Guex
,
A.
Bartolini
,
F.
Jourdan
,
J. E.
Spangenberg
,
J.-C.
Vicente
, and
U.
Schaltegger
.
2014
.
Evaluating the temporal link between the Karoo LIP and climatic–biologic events of the Toarcian Stage with high-precision U–Pb geochronology
.
Earth and Planetary Science Letters
 
408
:
48
56
.
Song
,
H.
,
P. B.
Wignall
,
D.
Chu
,
J.
Tong
,
Y.
Sun
,
H.
Song
,
W.
He
, and
Li
Tian
.
2014
.
Anoxia/high temperature double whammy during the Permian–Triassic marine crisis and its aftermath
.
Scientific Reports
 
4
:
srep04132
.
Steele-Petrović
,
H. M.
1976
.
Brachiopod food and feeding processes
.
Palaeontology
 
19
:
417
436
.
Suan
,
G.
,
B.
Pittet
,
I.
Bour
,
E.
Mattioli
,
L. V.
Duarte
, and
S.
Mailliot
.
2008
.
Duration of the Early Toarcian carbon isotope excursion deduced from spectral analysis: consequence for its possible causes
.
Earth and Planetary Science Letters
 
267
:
666
679
.
Suan
,
G.
,
E.
Mattioli
,
B.
Pittet
,
C.
Lécuyer
,
B.
Suchéras-Marx
,
L. V.
Duarte
,
M.
Philippe
,
L.
Reggiani
, and
F.
Martineau
.
2010
.
Secular environmental precursors to Early Toarcian (Jurassic) extreme climate changes
.
Earth and Planetary Science Letters
 
290
:
448
458
.
Suan
,
G. B.
,
van de Schootbrugge
,
T.
Adatte
,
J.
Fiebig
, and
W.
Oschmann
.
2015
.
Calibrating the magnitude of the Toarcian carbon cycle perturbation
.
Paleoceanography
 
30
:
495
509
.
Them
,
T. R.
, II
,
B. C.
Gill
,
A. H.
Caruthers
,
D. R.
Gröcke
,
E. T.
Tulsky
,
R. C.
Martindale
,
T. P.
Poulton
, and
P. L.
Smith
.
2017
.
High-resolution carbon isotope records of the Toarcian Oceanic Anoxic Event (Early Jurassic) from North America and implications for the global drivers of the Toarcian carbon cycle
.
Earth and Planetary Science Letters
 
459
:
118
126
.
Trecalli
,
A.
,
J.
Spangenberg
,
T.
Adatte
,
K. B.
Föllmi
, and
M.
Parente
.
2012
.
Carbonate platform evidence of ocean acidification at the onset of the early Toarcian oceanic anoxic event
.
Earth and Planetary Science Letters
 
357–358
:
214
225
.
van Hinsbergen
,
D. J. J.
,
L. V.
de Groot
,
S. J.
van Schaik
,
W.
Spakman
,
P. K.
Bijl
,
A.
Sluijs
,
C. G.
Langereis
, and
H.
Brinkhuis
.
2015
.
A paleolatitude calculator for paleoclimate studies
.
PLoS ONE
 
10
:
e0126946
.
Vaquer-Sunyer
,
R.
, and
C. M.
Duarte
.
2008
.
Thresholds of hypoxia for marine biodiversity
.
Proceedings of the National Academy of Sciences USA
 
105
:
15452
15457
.
Voje
,
K. L.
2018
.
Assessing adequacy of models of phyletic evolution in the fossil record
.
Methods in Ecology and Evolution
 
9
:
2402
2413
.
Vörös
,
A.
2002
.
Victims of the Early Toarcian anoxic event: the radiation and extinction of Jurassic Koninckinidae (Brachiopoda)
.
Lethaia
 
35
:
345
357
.
Vörös
,
A.
,
Á. T.
Kocsis
, and
J.
Pálfy
.
2016
.
Demise of the last two spire-bearing brachiopod orders (Spiriferinida and Athyridida) at the Toarcian (Early Jurassic) extinction event
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
457
:
233
241
.
Widdicombe
,
S.
, and
J. I.
Spicer
.
2008
.
Predicting the impact of ocean acidification on benthic biodiversity: what can animal physiology tell us?
Journal of Experimental Marine Biology and Ecology
 
366
:
187
197
.
Wignall
,
P. B.
, and
D. P. G.
Bond
.
2008
.
The end-Triassic and Early Jurassic mass extinction records in the British Isles
.
Proceedings of the Geologists’ Association
 
119
:
73
84
.
Wignall
,
P. B.
,
R. J.
Newton
, and
C. T. S.
Little
.
2005
.
The timing of paleoenvironmental change and cause-and-effect relationships during the early Jurassic mass extinction in Europe
.
American Journal of Science
 
305
:
1014
1032
.
Zhang
,
Y.
,
G. R.
Shi
,
W.-H.
He
,
H.-T.
Wu
,
Y.
Lei
,
K.-X.
Zhang
,
C.-C.
Du
,
T.-L.
Yang
,
M.-L.
Yue
, and
Y.-F.
Xiao
.
2016
.
Significant pre-mass extinction animal body-size changes: evidences from the Permian–Triassic boundary brachiopod faunas of South China
.
Palaeogeography, Palaeoclimatology, Palaeoecology
 
448
:
85
96
.
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.