Reductions in body size of benthic macroinvertebrates as a precursor of the early Toarcian (Early Jurassic) extinction event in the Lusitanian Basin, Portugal

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 CO₂ 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 CO₂ through the absorption and dissolution of CO₂ and SO₂ 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 CO₂ 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. 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. 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. 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.

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 log₂ 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:

(1)$${\rm GeoMea}{\rm n}_{{\rm Bivalves}} = \sqrt {L{\rm {^\ast}}H} $$

(2)$${\rm GeoMea}{\rm n}_{{\rm Brachiopods}} = \sqrt {L{\rm {^\ast}}W} $$

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 1–3.

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.

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).

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.

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 log₂ 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 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.

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 log₂ 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.

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.

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 CO₂ 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 CO₂ 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 pCO₂ 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.