Ocean warming affected faunal dynamics of benthic invertebrate assemblages across the Toarcian Oceanic Anoxic Event in the Iberian Basin (Spain)

The Toarcian Oceanic Anoxic Event (TOAE; Early Jurassic, ca. 182 Ma ago) represents one of the major environmental disturbances of the Mesozoic and is associated with global warming, widespread anoxia, and a severe perturbation of the global carbon cycle. Warming-related dysoxia-anoxia has long been considered the main cause of elevated marine extinction rates, although extinctions have been recorded also in environments without evidence for deoxygenation. We addressed the role of warming and disturbance of the carbon cycle in an oxygenated habitat in the Iberian Basin, Spain, by correlating high resolution quantitative faunal occurrences of early Toarcian benthic marine invertebrates with geochemical proxy data (δ18O and δ13C). We find that temperature, as derived from the δ18O record of shells, is significantly correlated with taxonomic and functional diversity and ecological composition, whereas we find no evidence to link carbon cycle variations to the faunal patterns. The local faunal assemblages before and after the TOAE are taxonomically and ecologically distinct. Most ecological change occurred at the onset of the TOAE, synchronous with an increase in water temperatures, and involved declines in multiple diversity metrics, abundance, and biomass. The TOAE interval experienced a complete turnover of brachiopods and a predominance of opportunistic species, which underscores the generality of this pattern recorded elsewhere in the western Tethys Ocean. Ecological instability during the TOAE is indicated by distinct fluctuations in diversity and in the relative abundance of individual modes of life. Local recovery to ecologically stable and diverse post-TOAE faunal assemblages occurred rapidly at the end of the TOAE, synchronous with decreasing water temperatures. Because oxygen-depleted conditions prevailed in many other regions during the TOAE, this study demonstrates that multiple mechanisms can be operating simultaneously with different relative contributions in different parts of the ocean.


Introduction
The early Toarcian Oceanic Anoxic Event (TOAE, ca. 182 Ma) [1], is one of the major carbon cycle perturbations of the Mesozoic [2], marked globally in the sedimentary, geochemical and paleontological records (e.g., [3,4]). Its ultimate cause is most likely the rapid and extensive increase in pCO 2 triggered by large-scale eruptions of the Karoo-Ferrar igneous province [5] and methane release from both oceanic [6] and terrestrial sources [7]. Consequences for marine ecosystems included increased water temperatures, widespread anoxia evidenced by black shale deposition, and possibly ocean acidification. The TOAE thus provides an opportunity to investigate the effects of climate warming on marine communities as recorded by fossils and may serve as a potential deep-time analog for the consequences of current and projected climate change on extant communities. The magnitude of seawater warming is estimated to have ranged between +3˚C and +7˚C, depending on latitude [8][9][10][11][12]. While dysoxic-anoxic conditions and black shales are typical for the TOAE in many regions, they can be very limited in extent and duration or even absent in other areas such as SW Europe and northern Africa [9,11, [13][14][15][16][17][18]. This implies different regional expressions of the TOAE [19,20]. Additional inferred effects of elevated pCO 2 are increased continental weathering rates [21], ocean acidification [22,23] and a crisis in carbonate production [24].
Experiments have shown that raised temperatures can lead to physiological stress for the marine ectotherms that are above their thermal optimum, through increased metabolic oxygen demand and reduced oxygen transport (e.g., [47,48]). Thus, temperature-driven stress may decrease population performance, leading to extirpations and ultimately to extinction [49]. What seems clear from the debate on the TOAE is that there is no universal trigger of the biotic crisis, and that regional differences in environmental conditions are important. As dysoxiaanoxia is not a plausible main stressor in oxygenated settings, the role of other factors should be tested statistically.
Here, we investigated the taxonomic and ecological responses of marine benthic macroinvertebrate assemblages across the TOAE in an oxygenated setting of the Iberian Basin, Spain, with three objectives: (i) Establish the patterns, timing, and magnitude of change in taxonomic composition and ecological structure of fossil communities from pre-TOAE times into the TOAE and post-TOAE intervals; (ii) test for negative faunal effects of elevated ocean temperatures and perturbations of the carbon cycle as recorded by geochemical proxy data; and (iii) interpret the results in terms of ecosystem stability, faunal recovery, and the relationship between temperature and faunal change.

Geological framework
Faunal and geochemical data were collected from the same samples along the sedimentary succession of Barranco de la Cañada (40˚23'53.4"N 1˚30'07.4"W), located near Albarracín, in the Iberian Basin, Spain (Fig 1). This section was chosen for its well-preserved and abundant benthic macroinvertebrate fauna, its relative stratigraphic completeness and its biostratigraphic control based on brachiopods and ammonites [14,50].
Ca. 28 m of the succession were sampled, spanning from the Pliensbachian/Toarcian boundary to the lower Bifrons Zone (middle Toarcian) (see sedimentary log in Fig 2). This time span corresponds to the Turmiel Formation [52], which is characterized by alternations of marlstones and argillaceous limestones. The depositional environment is interpreted as a low-energy, mid-ramp setting at 40-70 m depth, below storm wave base [14,15,53]. Rare packstones and rudstones indicate episodes of increased water energy. The absence of black shales, the relatively abundant fauna with widespread bioturbation across the TOAE interval [20,54] and the generally low TOC (<1.2 wt. %) measured for the formation [40] are evidence of oxygenated bottom conditions during the event. The sequence stratigraphic architecture of the basin has been subdivided into second and third order cycles [14,15,53,55]. The studied interval can be assigned to the second-order transgressive-regressive cycle termed LJ-3 in ref. [53] (Fig 3), which spans from the early Pliensbachian to the middle Toarcian.
The succession at Barranco de la Cañada shows the well-known TOAE negative carbon isotopic excursion which ranges from the uppermost Tenuicostatum Zone up to the lower Serpentinum Zone ( [57], Fig 3). This negative carbon isotope excursion is recorded globally from bulk rock carbonates, organic matter, wood and biogenic calcite (e.g., [12,22,59,60]). Its duration is estimated as *300 to 500 kyr according to recent studies [61,62], but uncertainty around those numbers remains (see ref. [63] for a summary). These changes in carbon isotope ratios reflect changes in the global carbon cycle, but the causes of the variation are complex, including the release of isotopically light volcanic, thermogenic, and/or biogenic carbon into the global ocean-atmosphere system, changes in primary productivity, and changes in the rates of organic matter burial [4].
Oxygen isotope ratios are available from the same samples at Barranco de la Cañada ( [57]; see 'Material and methods' below). δ 18 O values in shells are controlled by ambient water temperatures and other factors [58] but, in the studied material, they reliably represent temperature signals ( [57]; see 'Discussion'). Accordingly, early Toarcian temperature increased rapidly at the onset of the TOAE and stayed high throughout the interval (Fig 3). Local bottom water temperatures were elevated by ca. 3.5˚C through the entire TOAE [57]. Although somewhat smaller compared to some previous studies (see above), this temperature change is considered reliable as it represents averages for pre-TOAE, TOAE and post-TOAE intervals from a large quantity (>400) of measurements with high stratigraphic resolution. Towards the end of the TOAE temperatures decreased but during the post-TOAE remained on average higher than during the pre-TOAE (Fig 3).

Material and methods
We performed quantitative bed-by-bed sampling from carbonate beds, standardized by consistently collecting~13 kg of bulk rock for each sample. We investigated the complete coverage of benthic macroinvertebrates, regardless of preservation quality. A total number of 3668 individuals belonging to 93 taxa was recorded, comprising 2078 bivalves (62 taxa), 1400 brachiopods (21 taxa), 152 gastropods (nine taxa) and 38 individuals of one coral species (Fig 2). Taxa were identified preferentially at species level, wherever preservation allowed. Specimens that could not reliably be identified at genus level were disregarded. This set of specimens has already been described and investigated for trends in body size in ref. [45].
We obtained oxygen and carbon isotopes from calcite shells of the best-preserved rhynchonellid brachiopods and oysters from each sample (see ref. [57] for distribution of samples). Sample extraction was performed with preparation needle or scalpel after sediment matrix, altered crusts and, in the case of brachiopods, the primary shell layer was removed. Sampling areas where the shell calcite visually appeared well preserved were targeted and analytical results further scrutinized by comparison to ultrastructure information for each sampled specimen and measurement of element/Ca ratios for each sample. The methods applied for the geochemical analyses are described in detail in ref. [57]. In short, samples of ca. 500 μg size were analysed using a Sercon 20-22 Gas Source isotope ratio mass spectrometer in continuous flow setup, with the reaction of carbonate with nominally anhydrous phosphoric acid under a He atmosphere taking place at 70˚C. Instrumental drift and biases were corrected with a twopoint calibration using in-house standards CAR (Carrara Marble, δ 13 C = +2.10 ‰; δ 18 O = -1.93 ‰) and NCA (Namibian Carbonatite, δ 13 C = -5.63 ‰; δ 18 O = -21.90 ‰). These standards were previously calibrated against international standards, ensuring accuracy of the data. Precision of the measurements (2 s.d.) was found to be ± 0.07 ‰ for δ 13 C and ± 0.15 ‰ for δ 18 O based on 125 repeat analyses of CAR measured together with the fossil calcite samples.
Fossil preservation differs by taxonomical group: shells are well preserved in brachiopods, while in bivalves preservation is either as shell or moulds (both internal or external). Gastropods and corals occur generally as internal moulds. Because originally aragonitic taxa are preserved as moulds, the species composition and abundance distribution of the shelly macrobenthos apparently were not strongly influenced by diagenetic processes. Moreover, the absence of preferential shell orientations and lack of size sorting [45] suggest that shell transport was minimal, shells were buried within their original habitat, and faunal composition was not heavily biased by taphonomic processes.
The specimens collected in the field are deposited at the Museo de Ciencias Naturales de Universidad de Zaragoza, Spain [64]

Faunal analyses
Faunal analyses are based on a database of raw abundances of species standardized by rock volume and were performed in the R environment (v. 3.5.3, [65]). They address potential changes in taxonomic composition as well as in the ecological structure of benthic communities. Accordingly, we assigned each species to a specific mode of life (MOL)-a unique combination of feeding strategy, tiering and motility level reflecting the realized ecospace of a community [66,67]. MOL assignations (Table 1) are based on the published literature or were inferred from living relatives. Counts of species and MOLs were used to calculate the diversity indices (see below). Relative abundances (percentages) were used when applying ordination and clustering methods. In all analyses, levels with less than 17 specimens were excluded.
Our aim was to visualize to which degree pre-TOAE, TOAE and post-TOAE intervals differ in terms of taxonomical and ecological composition. Q-mode hierarchical clustering ('hclust' function of 'stats' package [65]) was based on a Bray-Curtis distance matrix ('vegdist' function in 'vegan' [68]), and Ward's clustering criterion [69] was implemented. The similarity profile test (SIMPROF) from the 'clustsig' package [70] was applied to determine the number of significant clusters [71]. Non-metric Multidimensional Scaling (NMDS) on a Bray-Curtis distance matrix was performed with the 'metaMDS' function from 'vegan' [68], with three dimensions and 100 attempts to search for a stable ordination solution. The fit of the NMDS was assessed by the stress value: the larger the stress value, the less efficient is the NMDS transformation [72]. Transformations with stress values lower than 0.2 were considered acceptable [72]. This ordination method is efficient in detecting distinct groupings in multivariate space and makes no assumptions about the distribution of the variables (e.g., [73]). The non-parametric Analysis of Similarities test (ANOSIM; 'anosim' function in 'vegan' [68], with 999 permutations) was used to test whether there is a significant difference between the three intervals represented in the NMDS.
We estimated the ecological importance of each MOL by calculating the average percentages per sample to visualize the relative abundance changes. The non-parametric Kruskal-Wallis test ( [74]; 'kruskal.test' function in 'stats' [65]) and the associated post hoc Dunn's test [75] ('dunnTest' function in 'FSA' [76])-with Bonferroni correction for multiple comparisonswere used to investigate differences in the median relative abundance between intervals. The post hoc test was performed only when the output of the Kruskal-Wallis test was significant. To avoid decreasing the statistical power of the tests, we refrained from performing these analyses in the cases a MOL accounted for insufficient (< 10) numbers of individuals per interval. The median values by interval were then calculated and used to estimate the difference in percentage of each MOL between intervals. This allowed us to illustrate the magnitude and direction of change of single MOLs. These same procedures were performed for the most abundant trait belonging to each of the three ecological variables (feeding, motility and tiering). For simplification, motility was categorized as either active or passive (i.e. all taxa with no ability to move).
We calculated diversity indices, which gave us different aspects of biodiversity change, to investigate the taxonomical and functional biodiversity trends in time: 1. Alroy's Shareholder Quorum Subsampling (SQS [56]; R code available at 'http://bio.mq. edu.au/~jalroy/SQS-3-3.R'). This method was developed to calculate how many species (and MOLs in the case of functional diversity) can be found in each sample given a certain "quorum" (i.e. the desired frequency coverage level) of the abundance distribution. This method corrects for changing abundance distributions, thus balancing the samples. To include as many samples as possible, a quorum of 0.8 (= 80% coverage) was used for the ecological analyses, and of 0.6 to calculate taxonomical diversity; 2. Simpson's evenness index calculated with the 'diversity' function in 'vegan' [68], method 'simpson'. This is based on the formula 1-D, where D = Sp 2 i and p i the proportional abundance of ith species; 3. Raw richness, i.e. the total count of species per sampled horizon, was calculated with the function 'specnumber' in 'vegan' [68]. This measure is already standardized by amount of bulk rock sampled in the field.

Correlation of faunal and geochemical data
To establish links of faunal variables with temperature change and/or processes determining the carbon cycle, we performed correlation tests for each of the biodiversity time series (SQSdiversity, Simpson's evenness, richness) and the compositional gradients (NMDS scores for axis 1 and 2) with the δ 18 O and δ 13 C isotope time series using Generalized Least Squares regression (GLS). The advantage of this method is that it can include terms to account for temporal autocorrelation among samples. Rows with missing values in any of the time series to be correlated were deleted prior to correlation. We performed several steps before the GLS proper ('gls' function in 'nlme' package [77]): i. The presence of any trend in time was investigated for each time series through a simple Ordinary Least Squares (OLS) regression (X~sampling level, where X is each faunal, ecological and geochemical time series) (results reported in S1 Table); ii. We checked for autocorrelation with the Durbin-Watson statistic ('durbinWatsonTest' function from 'car' package [78]) of the relationship X~δ 18 O + δ 13 C, where X refers to biodiversity indices or NMDS axis scores (S2 Table); iii. As some of the time series indicate non-stationarity (step i) and/or autocorrelation (step ii), we performed generalized differencing for the purpose of detrending [79] with the R script available at 'http://www.graemetlloyd.com/pubdata/functions_2.r'; iv. We repeated step ii) to check that the model assumptions of stationarity and non-autocorrelation were now satisfied (results of the Durbin-Watson statistic not reported).
v. We fitted an AutoRegressive Integrated Moving Average (ARIMA) process with the auto. arima function in the "forecast" R package [80,81] to the residuals of the differenced OLS regressions X~δ 18 O + δ 13 C. Because it is possible that unit roots are present as a different facet of an ARIMA process, the residual series of the fitted ARIMA process was investigated with the 'ndiffs' function from the "forecast" package ( [80,81], with 'kpss' and 'adf' stationarity tests and maximum fifth order differencing applied). Given the negative outcome, we proceeded with the last step; vi. The ARIMA fit was finally applied in the GLS regression to specify the expected correlation structure of the residuals.

Results
Bivalves dominate the studied fauna in terms of both species and number of specimens, with brachiopods being the second most abundant group (Fig 2). Gastropods and corals are rare faunal components (Table 1). Bivalves are most common and diverse in the pre-TOAE, but their diversity decreases approaching the TOAE. Many taxa, despite being rare or absent in the TOAE, range through and are again found in the post-TOAE. In contrast, brachiopods experience the extinction of pre-TOAE taxa in the lowermost Serpentinum Zone and a complete turnover leading to a new brachiopod fauna in the post-TOAE (Fig 2). The pre-TOAE and TOAE intervals are characterized by high abundances in just a few or only one species, whereas species abundance distributions are more even in the post-TOAE. Ecologically, two MOLs are dominant: free-lying epifaunal suspension feeders in the pre-TOAE, and pedically attached epifaunal suspension feeders in the TOAE and in the post-TOAE.

Change in taxonomic composition across the TOAE
Overall, the raw number of species increases from interval to interval, with 49 species sampled in the pre-TOAE, 55 in the TOAE and 68 in the post-TOAE. Moderately high and relatively stable diversity characterizes the pre-TOAE (Fig 3). Evenness and SQS-diversity vary little during the pre-TOAE, while raw richness increased during the Semicelatum Subzone until the onset of the TOAE. The TOAE itself is characterized by generally low values of raw species richness and contains the samples with the lowest evenness values. Fluctuations in all biodiversity indices are more pronounced compared to the other intervals (Fig 3). A shift to relatively high values in all three diversity metrics took place during uppermost TOAE and earliest post-TOAE times. Diversity values remain high and, in the case of evenness, stable throughout the studied post-TOAE interval.
Concerning faunal composition in more detail, cluster analysis on the 46 sampled horizons yielded 12 faunal associations, each consisting of faunal samples with similar taxonomic composition and abundance distribution of species (Fig 4A; version with stratigraphic levels provided in the Supplement as S1A Fig). The pre-TOAE, TOAE, and post-TOAE intervals are generally well separated, which is confirmed by the NMDS ordination ( Fig 5A) and NMDS axes scores (Fig 3) (ANOSIM: R = 0.701, p-value = 0.001). Whereas some overlap in NMDS space exists between pre-TOAE and TOAE intervals, the post-TOAE samples are distinct from those of both other intervals.
The pre-TOAE is taxonomically characterized by the dominance of the oyster Gryphaea sublobata, and secondarily by terebratulid brachiopods belonging to the genus Lobothyris ( Fig  2). These taxa are still an important part of the fauna during the early stages of the TOAE, hence the similarity of some TOAE assemblages with those from the pre-TOAE (Figs 2 and  5A). After the typical pre-TOAE brachiopods disappeared during the lower part of the TOAE (Fig 2), the rhynchonellid brachiopod Soaresirhynchia bouchardi becomes the dominant taxon in the upper part of the TOAE. The levels strongly dominated by S. bouchardi correspond to the evenness minima in Fig 3. The most common bivalves of the TOAE interval are Entolium corneolum and Pinna cf. folium. Brachiopods are still an important faunal component in the post-TOAE interval (Fig 2) and are represented mostly by terebratulid species of the genus Telothyris. As for the bivalves, Grammatodon concinnus, E. corneolum and Protocardia striatula are the most common species.
In terms of absolute abundances, numbers of individuals are high during the pre-TOAE and decline with the onset of the TOAE, reaching lowest values in the upper part of the TOAE (Fig 6). As an exception, two TOAE samples exhibit particularly high abundances, owing to high abundances of the rhynchonellid brachiopod Soaresirhynchia bouchardi (see below). Abundances again increase in the post-TOAE but mostly stayed below pre-TOAE levels.

Changes in ecological composition
The trajectories of functional diversity (Fig 6) resemble those based on taxonomic composition. Specifically, the TOAE is distinct from the other two intervals by exhibiting marked fluctuations in all curves. As was the case for taxonomic composition, pre-and post-TOAE intervals also differ in ecological composition (Figs 4B and 5B) (ANOSIM: R = 0.488, p-value = 0.001). The ecologically fairly homogeneous pre-TOAE (Fig 4B) is dominated by freelying (ca. 50%) and pedically attached (ca. 20%) epifaunal suspension feeders (Fig 7), represented respectively by ostreid bivalves and terebratulid brachiopods. The TOAE marks a major decline of epifaunal free-lying suspension feeders and an increase in pedically attached suspension feeders (Fig 7). The latter, represented by two species of the brachiopod genus Lobothyris in the early part of the TOAE and S. bouchardi after the brachiopod turnover, are dominant in the TOAE (> 40%), while free-lying epifaunal individuals are reduced to < 10%.
In ecological composition, TOAE samples often overlap with post-TOAE samples (Figs 4B and 5B) and the dominant MOL in the TOAE is still the most important in the post-TOAE (on average ca 40% of the post-TOAE fauna are pedically attached brachiopods). No MOL present in the pre-TOAE is lost after the TOAE.
By calculating the change in the percentage of each MOL among intervals (Fig 8), we confirm that the most substantial changes take place within two MOLs, i.e. the free-lying and the pedically attached suspension feeders. In particular, changes in these two MOLs are statistically significant for the pre-TOAE interval when compared to both other intervals ( Table 2). Most of the change occurred in the early phase of the TOAE, whereas subsequent changes during the post-TOAE are relatively minor in extent (Figs 7 and 8). When we considered motility, tiering, and feeding mechanism separately to pinpoint which traits within a MOL are most affected (Fig 9), we find that most change is located within taxa with passive motility. Changes in the relative abundance of epifaunal taxa and of suspension feeders are fairly small and significant only for suspension feeders ( Table 2).

Correlation of faunal data with δ 13 C and δ 18 O values
All taxonomically and functionally based biodiversity indices are significantly correlated with the δ 18 O isotope time series whereas correlations with δ 13 C values are all not significant (Table 3). Furthermore, the abundance distribution of MOLs is strongly tied to δ 18 O values, as is evident from significant correlations with both NMDS axes scores (Table 3). These results suggest an important role of temperature change in the observed faunal patterns. It should be kept in mind that reconstructing the carbon cycle from the δ 13 C record is much more complex than a comparatively simple conversion of δ 18 O values to temperature estimates. Therefore, we cannot exclude that carbon cycle-relevant processes did have some impact on communities but their expression in δ 13 C values is non-unique and only one of the contributing factors in shaping the observed δ 13 C signal.

The role of warming versus other environmental stressors
We document substantial faunal changes across the studied early Toarcian interval that are correlated with changes in water temperature. For another locality of the Iberian Range (Castrovido), Danise et al.
[9] found significant correlations of δ 18 O values with the taxonomic composition of Pliensbachian-Toarcian macrobenthic communities and bivalve subcommunities. Although community composition in our analysis was the only paleoecological variable that did not significantly correlate with the δ 18 O time series (Table 3), this finding is an independent indicator that temperature variations have influenced the composition of early Toarcian communities. Also, unlike our study, correlations between δ 18 O values and diversity indices (richness, evenness) were not significant for Castrovido and another Spanish locality (Sierra Palomera; [9]). This lack of correlations might be influenced by their relatively low number of observations for the TOAE interval and the stratigraphic clustering of brachiopodderived isotope values.
Apart from heat stress during the TOAE hyperthermal, the observed faunal patterns could also be caused by other or additional stressors related to climate warming. These include seawater dysoxia, reduced salinity, ocean acidification, changes in primary productivity, and

PLOS ONE
Responses of marine benthic macroinvertebrates to ocean warming during the Toarcian Oceanic Anoxic Event changes in habitat via sea level change. For the studied site, we discussed these factors in detail in two recent articles [45,57] and concluded that changes in temperature are the best   supported and most plausible factor. In brief, the well bioturbated and shell-rich TOAE sediments preclude oxygen-depleted bottom waters at Barranco de la Cañada. This is quite unlike more northern areas where black shales are common, and faunal change is related to ocean deoxygenation and changes in primary productivity (e.g., [28,29,82]). Freshening of seawater  (Table 2). White and dark grey bars represent positive and negative change, respectively. "n" = number of species, "N" = number of individuals per interval.
https://doi.org/10.1371/journal.pone.0242331.g009 Table 3. Generalized Least Squares (GLS) correlations of geochemical proxy data with paleoecological variables.  is incompatible with paleo-oceanographic modelling and with faunal composition (stenohaline ammonites, crinoids, and rhynchonelliform brachiopods are continuously present), whereas evidence of freshening exists for central and northern Europe [19,83], making salinity a stressor of regional importance. This geographic variability in salinity and oxygen regimes shows that, on larger spatial scales, multiple mechanisms can be operating simultaneously with different relative contributions in different parts of the ocean. Ocean acidification is more difficult to assess but is not supported by faunal patterns at Barranco de la Cañada, where brachiopods, which might be more tolerant against lowered pH than bivalves [84,85], are more strongly affected ( [45], this study). Likewise, the selective extinction of brachiopods may argue against low productivity because their typically low metabolic rates and their capability to feed on both particulate and dissolved organic matter should buffer them against reduced food supply [86,87]. Finally, patterns of faunal change may also arise from facies change owing to fluctuations in sea level, and Danise et al.

p-value
[9] showed the clustering of last occurrences of early Toarcian species in the Iberian Basin at flooding surfaces. At Barranco de la Cañada, facies changes are fairly subtle and we sampled largely homogeneous lithologies (marly wackestones to floatstones) with just one sample being from a rudstone (S1 Fig). Overall, the depositional environment stayed in a low-energy, mid-ramp position. Also, diversity trends (both taxonomical and ecological) are not consistent in equivalent stages of transgressive-regressive cycles. This is for example evident in Fig 3 where we would expect congruent diversity trends in each transgressive-regressive cycle if diversity fluctuations were primarily driven by fluctuations in sea level, which is not the case.

Timing of faunal change
Changes in multiple community attributes (Figs 3 and 6) took place at the onset of the TOAE, synchronous with an increase in water temperatures. Within the first samples of the TOAE, taxonomical and functional diversity indices decrease and enter into a mode of strongly fluctuating values. Similarly, the two most important MOLs show drastic abundance changes right at the onset of the event (Fig 7). This suggests that in terms of ecological organisation the preexisting ecosystem was not resistant to the geologically rapid temperature rise. By contrast, taxonomic composition shifted more gradually (see NMDS scores in Fig 3), and the main change happened around the transition between the early and the late phase of the TOAE, coinciding with the final disappearance of the pre-existing brachiopod fauna and the appearance of Soaresirhynchia assemblages. Similar patterns at the local scale are also known from modern ecosystems where-on much shorter time scales-anthropogenic disturbances predominantly decreased abundance and species richness [88]. In contrast to our findings, however, the meta-analysis of Murphy & Romanuk [88] showed that the effects of temperature increase on modern aquatic ectotherms were insignificant. This difference could be due to the magnitude and prolonged duration of raised temperatures in our study system or because increases in temperature need not be deleterious as long as organisms shift towards or stay close to their thermal optima. Our observation of functional shifts at the onset of the TOAE is consistent with the timing of functional changes in dysoxic to anoxic settings [28,29,31,89]. However, apart from functional change, these dysoxic-anoxic environments also exhibit compositional shifts at the onset of the TOAE negative carbon isotopic excursion. This contrasts with our present findings from an oxygenated environment, where compositional changes shifted into the early phase of the TOAE. This discrepancy could suggest that the combination of warming and dysoxia in oxygen-depleted environments led to earlier, more abrupt shifts comprising both community composition and structure.

Ecological state shift during hyperthermal conditions
The TOAE interval clearly differs from the intervals before and after in several faunal characteristics that are indicative of sustained stressful and ecologically unstable conditions in the benthic marine ecosystem. First, the species composition changed with the complete loss of the pre-TOAE brachiopod fauna in the first half of the TOAE. At a wider geographic scale, none of the pre-TOAE articulate brachiopod species of the Iberian platform system survived this episode [40]. In the second half of the TOAE, Soaresirhynchia-dominated assemblages appeared and prevailed, whereas most bivalve species disappeared or became rare. The dominance of Soaresirhynchia in the late phase of the TOAE is known from other basins in the Mediterranean [9, [15][16]40,44,90,91]. S. bouchardi is a small-sized species with high morphological variability, geographically widespread (ranging from North Africa to England [40]), and occurs in high numbers. This taxon was successful especially during the rising limb of the TOAE negative carbon isotope excursion where it thrived despite persisting high water temperatures. Its presumably slow growth and low metabolic rates, as suggested by the geochemical composition and low organic matter content of the shell [57], probably put it at an advantage in a stressed environment over other species with higher energy demands [9,57]. Also, the associated bivalve Parvamussium pumilum and the lingulid brachiopod Lingularia sp. have been interpreted as opportunistic taxa capable of surviving in stressed environments [92,93], making the second half of the TOAE an interval dominated by opportunists.
Second, the abundance of specimens in samples (which is standardized by rock volume) is on average lower during the TOAE than in the pre-TOAE. Combined with the observation from the same samples that larger-sized species selectively become less abundant during the TOAE [45], this suggests an overall decline of macrobenthic biomass in the TOAE. Notable exceptions are two samples with proliferations of the opportunistic S. bouchardi.
Third, the shift from diverse pre-TOAE assemblages to the paucispecific assemblages of the TOAE is paralleled by distinct oscillations of both taxonomic and functional biodiversity indices during the TOAE (Figs 3 and 6). Along with pronounced fluctuations in MOL-based NMDS scores, in the proportions of passive and epifaunal individuals, and with high-amplitude excursions in the abundance of specimens (Figs 6 and 9), this points to ecological instability during the TOAE.

Recovery from the TOAE hyperthermal
The latest phase of the TOAE is characterized by a swift recovery of both taxonomical and ecological biodiversity. Such a geologically quick recovery contrasts with other areas where the TOAE is characterized by black shales (e.g., [28,29,31,89]) and has been related to the more favourable conditions of Tethyan areas where anoxia and dysoxia did not develop [9]. The return to high diversity values in the latest TOAE is synchronous with the shift to more positive δ 18 O values, with cooler temperatures indicating the end of the hyperthermal interval. We interpret this synchronicity of trends in temperature and biodiversity at the beginning and the end of the TOAE hyperthermal as evidence for an important role of warming in diversity loss and of cooling in the subsequent recovery of diversity. The speed of diversity recovery as soon as temperatures decline suggests that benthic diversity was resilient, albeit taxonomic composition had changed. Subsequent stable diversity in the post-TOAE interval indicates a return to ecosystem stability (Figs 3 and 6).
Comparing pre-and post-TOAE assemblages, we found that no MOL present in the pre-TOAE was lost after the TOAE. Yet, ordinations and cluster analyses (Figs 4 and 5) show that pre-and post-TOAE intervals differ taxonomically and ecologically from each other. Taxonomic disparity can readily be explained by the strong turnover in brachiopods and marked reductions in the abundance of the oyster Gryphaea sublobata that had dominated the pre-TOAE interval. In the face of extinction, full recovery of taxonomic composition remains impossible. Yet, functional composition can be restored in the absence of taxonomic compositional recovery because different species can be functionally the same [26]. Post-TOAE temperatures were mostly higher compared to the pre-TOAE ( [57] , Fig 3). To test whether this difference in temperature explains the difference in ecological composition, we compared only those post-TOAE samples with pre-TOAE assemblages where δ 18 O values stayed within their pre-TOAE range (see selected samples in Fig 4 and S1 Fig). Even these selected post-TOAE assemblages have a different ecological composition than before which shows that communities did not return to their previous state despite similar inferred temperatures.

Comparison with modern ecosystems and projected warming scenarios
The early Toarcian sedimentary rocks of the Barranco de la Cañada section were deposited in a shallow to mid shelf/ramp environment at low, subtropical paleolatitudes. The dominance of bivalves and brachiopods in marine, shelly, macrobenthic associations in these level bottom ecosystems is common during the early Jurassic whereas the higher level taxonomic composition has changed to bivalve-gastropod associations in the late Cenozoic [94]. Although locally common in shallow waters, the distribution of living brachiopods reflects the continuation of a post-Permian trend in brachiopod retreat to offshore habitats, and today they typically are species of the deeper continental shelf and upper slope [86]. Because brachiopods in our study were more severely hit by temperature rise than bivalves one might expect that, everything else being equal, modern shallowto-mid shelf assemblages are possibly more resistant to warming than those of the early Jurassic.
In terms of absolute temperatures, using the brachiopod oxygen isotope thermometer of Brand et al. [58], we estimated average pre-TOAE water temperatures at the seafloor of ca. 21˚C followed by a geologically rapid mean temperature increase during the TOAE of ca. 3.5˚C [57]. Climate change predictions until the end of the century range from a low greenhouse gas emission scenario with a projected mean of 0.73˚C warming of global mean sea surface temperature to a high emission scenario with mean warming of 2.58˚C relative to preindustrial levels [95]. In this respect, the amount of early Toarcian warming inferred in our analysis [57] is beyond the range of ocean warming scenarios for the near future. Yet, this does not imply that the local ecological effects of future warming are less severe than those established for the early Toarcian. These temperature estimates are not equivalent because they represent local point estimates (Barranco de la Cañada) versus global mean sea surface temperature (in ref. [95]). Furthermore, the rate of projected warming is likely much higher than during the TOAE, making adaptations to warming oceans more difficult. Several other factors add to the uncertainties in assessing the consequences of current and projected climate change for marine ecosystems from deep-time warming analogs as recorded by fossils. These include the time-averaged nature of fossil assemblages, the complex interactions with other warming-related stressors which may vary geographically, differences in the thermal structure of communities, and differences in the starting state of the climate system [96][97][98]. Despite these complexities and uncertainties we tentatively consider the identified faunal changes across the TOAE hyperthermal-taxonomic and ecological reorganisations; declines in diversity, abundance, and biomass; long-term ecological instability with predominance of opportunistic species; and reductions in community size structure (see ref. [45])-to also be plausible threats to present-day shallow-water benthic communities.

Conclusions
High-resolution early Toarcian faunal and geochemical data from marine mid-ramp environments at Barranco de la Cañada provide statistical support for a relationship between bottom water temperatures and taxonomical and ecological changes in benthic macroinvertebrate assemblages. This result provides strong support of the regional importance of warming as a proximate cause of the TOAE-related biotic crisis. Because dysoxia-anoxia prevailed in many other regions during the TOAE, our study alongside others suggests that multiple mechanisms can be operating simultaneously with different relative contributions in different parts of the ocean.
The local faunal assemblages were reorganized taxonomically and ecologically across the TOAE. This suggests that the local ecosystem could not withstand the disturbance, heat stress exceeded the communities' resistance capacity, and an environmental threshold was crossed leading to a non-transient change with a new distribution of species-specific traits in the TOAE aftermath. Brachiopods were most severely hit, as they experienced species extinctions in the lower part of the TOAE followed by the proliferation of an opportunistic species and a complete species turnover in the post-TOAE interval. Molluscs, on the contrary, mostly ranged through the event, hinting at rather different responses to environmental stress among taxonomical groups. Most of the ecological changes were focused at the onset of the TOAE hyperthermal and comprise declines in multiple diversity metrics, abundance, and biomass. During the TOAE, the studied marine communities were ecologically unstable, as evidenced by oscillatory dynamics in most ecological variables, whereas stable and diverse post-TOAE faunal assemblages were established simultaneously with decreasing water temperatures at the end of the TOAE. This suggests geologically rapid functional recovery to an undisturbed state, albeit with different frequency distributions of modes of life than before the event.
This confirms the need for a better understanding of ecosystem and organism responses under heat stress. Investigations of deep-time episodes of climate change like the TOAE can provide insights how marine communities might ultimately respond in the face of the present global warming and to changes in future climatic regimes.