Veronica Piazza1, Clemens V Ullmann2, Martin Aberhan1. 1. Museum für Naturkunde, Leibniz Institute for Evolution and Biodiversity Science, Berlin, Germany. 2. University of Exeter, Camborne School of Mines, College of Engineering, Mathematics and Physical Sciences, Penryn, Cornwall, United Kingdom.
Abstract
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.
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.
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 pCO2 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-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-18]. This implies different regional expressions of the TOAE [19,20]. Additional inferred effects of elevated pCO2 are increased continental weathering rates [21], ocean acidification [22,23] and a crisis in carbonate production [24].The patterns and rates of ecosystem recovery from environmental stress during mass extinctions are as yet insufficiently studied [25]. However, understanding community stability when facing disturbances is central to inform present-day conservation management [26]. Biodiversity loss and extinctions in fossil communities might provide a long-term perspective of ecosystem response to global warming. The TOAE coincides with elevated extinction rates for nektonic and benthic organisms, i.e. molluscs, brachiopods, ostracods, foraminifera and dinoflagellates, and with changes in community structure (e.g., [18,24,27-33]. The proximate causes are debated and most studies see widespread dysoxia–anoxia as the main stressor (e.g., [5,24,27-30,34-37], with Them et al. [38] even arguing for global ocean deoxygenation to have already started at the Pliensbachian/Toarcian boundary. Yet, the record of dysoxia–anoxia is far from ubiquitous as is evident from faunal and facies analyses in oxygenated environments from SW Europe [9,11,13-15,17,39]. Even from restricted basins where anoxia developed, oxygenated bottom conditions are recorded before the onset of the TOAE (e.g., [29]). In the absence of low-oxygen conditions, warming was suggested as an important stressor (e.g., [11,40-43]), but rigorous statistical tests of the relationship between temperature change and faunal change are few (e.g., [9,44-46]).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 dysoxia–anoxia 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].
Fig 1
Geographical location and paleogeographical position of Barranco de la Cañada.
Location map of the investigated section in the Iberian Range system (A) and paleogeographical reconstruction of the western Tethyan basins (B). Studied area marked by the black dot. The grey shading in (A) represents the Iberian Chain. Abbreviations in (B): LB = Lusitanian Basin; IM = Iberian Massif; IB = Iberian Basin. (A) and (B) modified respectively after Figs 1 and 6 in ref. [51] under a CC BY license, with permission from Elsevier, original copyright 2015.
Geographical location and paleogeographical position of Barranco de la Cañada.
Location map of the investigated section in the Iberian Range system (A) and paleogeographical reconstruction of the western Tethyan basins (B). Studied area marked by the black dot. The grey shading in (A) represents the Iberian Chain. Abbreviations in (B): LB = Lusitanian Basin; IM = Iberian Massif; IB = Iberian Basin. (A) and (B) modified respectively after Figs 1 and 6 in ref. [51] under a CC BY license, with permission from Elsevier, original copyright 2015.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.
Fig 2
Stratigraphy and ranges of benthic macroinvertebrate species at Barranco de la Cañada.
Taxa are ordered by last occurrences, and subdivided into taxonomical groups. Each circle represents the presence of a taxon at a sampled level, the size reflecting the abundance of the taxon. The star-shaped symbols indicate the sampled levels. The shaded grey area represents the extent of the chemostratigraphically defined Toarcian Oceanic Anoxic Event (TOAE) interval. Chronostratigraphic boundaries modified from ref. [14]. The horizontal dashed line marks the level with the last appearances of the pre-TOAE brachiopod fauna.
Fig 3
Stratigraphy, transgressive-regressive cycles, geochemistry, and paleoecology at Barranco de la Cañada.
Biodiversity metrics (richness, evenness, and Shareholder Quorum Subsampling [SQS] diversity) are based on the taxonomic composition of faunal samples. Non-metric multidimensional scaling (NMDS) curves represent the values of axis 1 and axis 2 of the NMDS ordination of Fig 5A. SQS-diversity calculated after Alroy [56]. Third order transgressive-regressive cycles from refs. [15,53]. Isotope data from ref. [57]; the temperature change ΔT was calculated relative to the lowermost value recorded in the section using Brand’s equation [58], following refs. [45,57]. Horizontal bars in the isotope time series represent the double standard error for each sample. For additional information, see caption of Fig 2.
Stratigraphy and ranges of benthic macroinvertebrate species at Barranco de la Cañada.
Taxa are ordered by last occurrences, and subdivided into taxonomical groups. Each circle represents the presence of a taxon at a sampled level, the size reflecting the abundance of the taxon. The star-shaped symbols indicate the sampled levels. The shaded grey area represents the extent of the chemostratigraphically defined Toarcian Oceanic Anoxic Event (TOAE) interval. Chronostratigraphic boundaries modified from ref. [14]. The horizontal dashed line marks the level with the last appearances of the pre-TOAE brachiopod fauna.
Stratigraphy, transgressive-regressive cycles, geochemistry, and paleoecology at Barranco de la Cañada.
Biodiversity metrics (richness, evenness, and Shareholder Quorum Subsampling [SQS] diversity) are based on the taxonomic composition of faunal samples. Non-metric multidimensional scaling (NMDS) curves represent the values of axis 1 and axis 2 of the NMDS ordination of Fig 5A. SQS-diversity calculated after Alroy [56]. Third order transgressive-regressive cycles from refs. [15,53]. Isotope data from ref. [57]; the temperature change ΔT was calculated relative to the lowermost value recorded in the section using Brand’s equation [58], following refs. [45,57]. Horizontal bars in the isotope time series represent the double standard error for each sample. For additional information, see caption of Fig 2.
Fig 5
Non-metric multidimensional scaling ordination (NMDS) of faunal samples at Barranco de la Cañada.
Results plotted by taxonomical (A) and ecological (B) composition. The symbols and the shading of the convex-hull polygons are coded by interval (pre-TOAE, TOAE and post-TOAE).
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). δ18O 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 two-point calibration using in-house standards CAR (Carrara Marble, δ13C = +2.10 ‰; δ18O = -1.93 ‰) and NCA (Namibian Carbonatite, δ13C = -5.63 ‰; δ18O = -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 δ13C and ± 0.15 ‰ for δ18O 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] with inventory numbers MPZ 2019/415 –MPZ 2019/571, MPZ 2019/792 –MPZ 2019/889, MPZ 2019/920 –MPZ 2019/942 and MPZ 2019/1041. All necessary permits were obtained for the described study, which complied with all relevant regulations. Approval of the conduction of paleontological fieldwork was obtained from the Direccion General de Cultura y Patrimonio of Zaragoza (Gobierno de Aragon).
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.
Table 1
Macrobenthic marine invertebrates recorded from the early Toarcian of Barranco de la Cañada, Spain, and their ecological attributes.
Abbreviations: Motility Level, stationary categories: Ped = Pedically attached; Bys = Byssate; Cem = Cemented; Fr(un) = Unattached free-lying; Motility Level, mobile categories: Fac-Mob(un) = Unattached facultative mobile; Fac-Mob(at) = Attached facultative mobile; Reg-Mov = Regularly moving. Tiering: Epi = Epifaunal; S-inf = Semi-infaunal; Shal-inf = Shallow infaunal; D-inf = Deep infaunal. Feeding Mechanisms: Susp = Suspension feeder; Her = Herbivore/grazing; Car = Carnivore/microcarnivore; Surf-dep = Surface deposit feeder.Numbering of taxa as in Fig 2.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 comparisons—were 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: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;Simpson’s evenness index calculated with the ‘diversity’ function in ‘vegan’ [68], method ‘simpson’. This is based on the formula 1- D, where D = Σp2i and pi the proportional abundance of ith species;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 (SQS-diversity, Simpson’s evenness, richness) and the compositional gradients (NMDS scores for axis 1 and 2) with the δ18O and δ13C 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]):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);We checked for autocorrelation with the Durbin-Watson statistic (‘durbinWatsonTest’ function from ‘car’ package [78]) of the relationship X ~ δ18O + δ13C, where X refers to biodiversity indices or NMDS axis scores (S2 Table);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’;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).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 ~ δ18O + δ13C. 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;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.
Fig 4
Cluster analyses of the faunal samples at Barranco de la Cañada.
Results plotted according to taxonomical composition using species as variables (A) and ecological composition using modes of life as variables (B). The identified associations are shaded in grey and each sample is coded by interval (pre-TOAE, TOAE and post-TOAE). The stars mark the levels in the post-TOAE where the recorded isotope values (both δ18O and δ13C) are within their pre-TOAE ranges.
Cluster analyses of the faunal samples at Barranco de la Cañada.
Results plotted according to taxonomical composition using species as variables (A) and ecological composition using modes of life as variables (B). The identified associations are shaded in grey and each sample is coded by interval (pre-TOAE, TOAE and post-TOAE). The stars mark the levels in the post-TOAE where the recorded isotope values (both δ18O and δ13C) are within their pre-TOAE ranges.
Non-metric multidimensional scaling ordination (NMDS) of faunal samples at Barranco de la Cañada.
Results plotted by taxonomical (A) and ecological (B) composition. The symbols and the shading of the convex-hull polygons are coded by interval (pre-TOAE, TOAE and post-TOAE).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.
Fig 6
Paleoecology of faunal samples at Barranco de la Cañada based on ecological composition.
Functional diversity indices (richness, evenness, and SQS-diversity) and NMDS axis scores of Fig 5B were calculated using modes of life as variables. Third order transgressive-regressive cycles from refs. [15,53]. The abundance of samples is shown as their respective number of individuals. For additional information, see captions of Figs 2 and 3.
Paleoecology of faunal samples at Barranco de la Cañada based on ecological composition.
Functional diversity indices (richness, evenness, and SQS-diversity) and NMDS axis scores of Fig 5B were calculated using modes of life as variables. Third order transgressive-regressive cycles from refs. [15,53]. The abundance of samples is shown as their respective number of individuals. For additional information, see captions of Figs 2 and 3.
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 free-lying (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%.
Fig 7
Relative abundance of each mode of life (MOL).
The shaded area represents the TOAE interval, and the dashed horizontal line marks the level with the last appearances of the pre-TOAE brachiopod fauna. Legend: “n” = number of species; “N” = number of individuals. MOL codings: 1 = Unattached free-lying epifaunal suspension feeders; 2 = Pedically attached epifaunal suspension feeders; 3 = Unattached facultatively mobile shallow infaunal suspension feeders; 4 = Byssate epifaunal suspension feeders; 5 = Unattached facultatively mobile epifaunal suspension feeders; 6 = Unattached facultatively mobile deep infaunal suspension feeders; 7 = Byssate semi-infaunal suspension feeders; 8 = Cemented epifaunal suspension feeders; 9 = Attached facultatively mobile epifaunal suspension feeders; 10 = regularly moving epifaunal herbivore; 11 = Unattached facultatively mobile epifaunal carnivore; 12 = Unattached free-lying epifaunal carnivore, 13 = Regularly moving shallow infaunal surface deposit feeders; 14 = Attached facultatively mobile shallow infaunal suspension feeders.
Relative abundance of each mode of life (MOL).
The shaded area represents the TOAE interval, and the dashed horizontal line marks the level with the last appearances of the pre-TOAE brachiopod fauna. Legend: “n” = number of species; “N” = number of individuals. MOL codings: 1 = Unattached free-lying epifaunal suspension feeders; 2 = Pedically attached epifaunal suspension feeders; 3 = Unattached facultatively mobile shallow infaunal suspension feeders; 4 = Byssate epifaunal suspension feeders; 5 = Unattached facultatively mobile epifaunal suspension feeders; 6 = Unattached facultatively mobile deep infaunal suspension feeders; 7 = Byssate semi-infaunal suspension feeders; 8 = Cemented epifaunal suspension feeders; 9 = Attached facultatively mobile epifaunal suspension feeders; 10 = regularly moving epifaunal herbivore; 11 = Unattached facultatively mobile epifaunal carnivore; 12 = Unattached free-lying epifaunal carnivore, 13 = Regularly moving shallow infaunal surface deposit feeders; 14 = Attached facultatively mobile shallow infaunal suspension feeders.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).
Fig 8
Relative change of each mode of life (MOL) among time intervals.
Comparisons are based on medians. White and dark grey bars represent positive and negative change, respectively. The symbol (*) denotes significant p-values from the Kruskal-Wallis and Dunn’s tests (Table 2). For MOL codings, see caption of Fig 7.
Table 2
Results of the non-parametric Kruskal-Wallis test on medians and post-hoc Dunn’s test.
Kruskal-Wallis test
Dunn’s test
pre-TOAE vs. TOAE
pre-TOAE vs. post-TOAE
TOAE vs. post-TOAE
MOL/Trait
χ2
p-value
Z
p-value
Z
p-value
Z
p-value
1
22.919
<0.001
4.520
<0.001
4.096
<0.001
-0.803
1
2
9.423
0.009
-2.877
0.012
-2.658
0.024
0.453
1
3
2.144
0.342
/
/
/
/
/
/
4
7.118
0.029
2.636
0.025
1.363
0.518
-1.634
0.307
5
0.844
0.656
/
/
/
/
/
/
6
1.177
0.555
/
/
/
/
/
/
7
7.690
0.021
-1.381
0.502
-2.733
0.019
-1.438
0.451
8
4.872
0.088
/
/
/
/
/
/
9
14.656
0.001
-1.416
0.470
-3.643
0.001
-2.430
0.045
10
15.892
<0.001
-0.681
1
-3.460
0.002
-3.110
0.006
11
12
13
14
20.113
<0.001
-3.734
0.001
-0.466
1
3.979
<0.001
Pas
4.994
0.082
/
/
/
/
/
/
Epi
1.830
0.401
/
/
/
/
/
/
Susp
22.185
<0.001
1.870
<0.001
4.523
<0.001
2.882
<0.001
Dunn’s post-hoc test was performed in case of a significant result from the Kruskal-Wallis test to investigate which pairwise-comparison by interval is significant, otherwise the symbol “/” was used. The reported p-values for the post-hoc test are adjusted after applying the Bonferroni correction for multiple comparisons. Analyses were performed for each MOL and single traits. Statistically significant values (p < 0.05) are in bold. χ2 values based on two degrees of freedom. For coding of MOLs (1–14) see caption of Fig 7. Pas = Passive; Epi = Epifauna; Susp = Suspension feeders.
Relative change of each mode of life (MOL) among time intervals.
Comparisons are based on medians. White and dark grey bars represent positive and negative change, respectively. The symbol (*) denotes significant p-values from the Kruskal-Wallis and Dunn’s tests (Table 2). For MOL codings, see caption of Fig 7.Dunn’s post-hoc test was performed in case of a significant result from the Kruskal-Wallis test to investigate which pairwise-comparison by interval is significant, otherwise the symbol “/” was used. The reported p-values for the post-hoc test are adjusted after applying the Bonferroni correction for multiple comparisons. Analyses were performed for each MOL and single traits. Statistically significant values (p < 0.05) are in bold. χ2 values based on two degrees of freedom. For coding of MOLs (1–14) see caption of Fig 7. Pas = Passive; Epi = Epifauna; Susp = Suspension feeders.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).
Fig 9
Relative abundance and relative change of the most important ecological traits.
Each depicted trait represents the most abundant category in each of the three ecological characters that define a MOL (motility, tiering and feeding mechanism). The change among intervals is calculated from the median abundance of all sample means for each interval. The symbol (*) denotes significant p-values from the Kruskal-Wallis and Dunn’s tests (Table 2). White and dark grey bars represent positive and negative change, respectively. “n” = number of species, “N” = number of individuals per interval.
Relative abundance and relative change of the most important ecological traits.
Each depicted trait represents the most abundant category in each of the three ecological characters that define a MOL (motility, tiering and feeding mechanism). The change among intervals is calculated from the median abundance of all sample means for each interval. The symbol (*) denotes significant p-values from the Kruskal-Wallis and Dunn’s tests (Table 2). White and dark grey bars represent positive and negative change, respectively. “n” = number of species, “N” = number of individuals per interval.
Correlation of faunal data with δ13C and δ18O values
All taxonomically and functionally based biodiversity indices are significantly correlated with the δ 18O isotope time series whereas correlations with δ13C values are all not significant (Table 3). Furthermore, the abundance distribution of MOLs is strongly tied to δ 18O 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 δ13C record is much more complex than a comparatively simple conversion of δ18O values to temperature estimates. Therefore, we cannot exclude that carbon cycle-relevant processes did have some impact on communities but their expression in δ13C values is non-unique and only one of the contributing factors in shaping the observed δ13C signal.
Table 3
Generalized Least Squares (GLS) correlations of geochemical proxy data with paleoecological variables.
p-value
Coefficients
Standard error
ARIMA (p,d,q)
δ18O
δ13C
δ18O
δ13C
δ18O
δ13C
A Taxonomic composition
SQS
(0,0,0)
0.047
0.598
2.595
-0.312
1.258
0.585
Simpson’s Evenness
(0,0,0)
0.014
0.535
0.212
0.024
0.082
0.0378
Richness
(0,0,0)
0.036
0.797
5.839
-0.320
2.659
1.236
NMDS1
(1,0,2)
0.915
0.428
-0.006
-0.018
0.054
0.022
NMDS2
(0,0,0)
0.290
0.336
0.085
0.036
0.079
0.0367
B Ecological composition
SQS
(0,0,0)
0.024
0.842
1.599
0.063
0.673
0.313
Simpson Evenness
(0,0,0)
0.012
0.703
0.203
0.014
0.076
0.035
Richness
(0,0,0)
0.032
0.890
1.791
0.052
0.800
0.372
NMDS1
(0,0,0)
0.017
0.791
0.379
-0.019
0.150
0.070
NMDS2
(0,0,0)
0.009
0.869
0.328
0.009
0.118
0.055
Paleoecological variable represented by SQS-diversity, evenness, richness, and NMDS axis scores. The geochemical proxy data used in this study are δ18O and δ13C values obtained from brachiopod and oyster shells. The correlations were performed for faunal assemblages characterized by taxonomic composition (A) and by ecological composition (B) on differenced time series. Statistically significant values (p < 0.05) are in bold. Abbreviations of the parameters of the ARIMA process fitted to the GLS model: p = number of time lags; d = degree of differencing; q = order of moving average.
Paleoecological variable represented by SQS-diversity, evenness, richness, and NMDS axis scores. The geochemical proxy data used in this study are δ18O and δ13C values obtained from brachiopod and oyster shells. The correlations were performed for faunal assemblages characterized by taxonomic composition (A) and by ecological composition (B) on differenced time series. Statistically significant values (p < 0.05) are in bold. Abbreviations of the parameters of the ARIMA process fitted to the GLS model: p = number of time lags; d = degree of differencing; q = order of moving average.
Discussion
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 δ18O 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 δ18O 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 δ18O 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 brachiopod-derived 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 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 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. [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 pre-existing 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 δ18O 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 δ18O 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 shallow-to-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 pre-industrial 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-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.
Cluster analyses of the sampled fauna at Barranco de la Cañada.
Clusters from Fig 4 plotted showing the stratigraphic level and lithology of each sample in meters above the Pliensbachian-Toarcian boundary. Results plotted by taxonomical (A) and ecological (B) composition. The identified associations are numbered and shaded in grey and each sample is coded by interval (pre-TOAE, TOAE and post-TOAE). The stars mark the levels in the post-TOAE where the recorded isotope values (both δ18O and δ13C) are within pre-TOAE ranges.(TIF)Click here for additional data file.
Results of the Ordinary Least Squares (OLS) correlations of faunal variables and of geochemical proxy data against time (sampling level).
Faunal variables are SQS-diversity, evenness, richness, and NMDS axis scores, whereas the geochemical proxy data are δ18O and δ13C values. Results are presented for faunal assemblages characterized by both taxonomic and ecological composition respectively and are based on the original time series. Statistically significant values (p < 0.05) are in bold.(DOCX)Click here for additional data file.
Durbin-Watson (D-W) statistics on Ordinary Least Squares (OLS) correlations of faunal variables with geochemical proxy data.
The OLS correlations are in the form X~ δ18O + δ13C, where X is the respective diversity index or NMDS axis score. Results are presented separately for faunal assemblages characterized by taxonomic and by ecological composition based on the original time series. Statistically significant values (p < 0.05) are in bold.(DOCX)Click here for additional data file.10 Jun 2020PONE-D-20-11804Temperature affects faunal dynamics of benthic invertebrate assemblages across the Toarcian Oceanic Anoxic Event in the Iberian Basin (Spain)PLOS ONEDear Dr. Piazza,Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.Both reviewers suggest minor revisions and provide excellent comments to improve the ms. Please address all reviewer comments in your reply.Please submit your revised manuscript by Jul 25 2020 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.Please include the following items when submitting your revised manuscript:A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocolsWe look forward to receiving your revised manuscript.Kind regards,David P. Gillikin, Ph.D.Academic EditorPLOS ONEJournal Requirements:When submitting your revision, we need you to address these additional requirements.1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found athttps://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf andhttps://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf2. In your Methods section, please provide additional information regarding the permits you obtained for the work. Please ensure you have included the full name of the authority that approved the field site access and, if no permits were required, a brief statement explaining why.3. During our internal checks, the in-house editorial staff noted that you conducted research or obtained samples in another country.Please check the relevant national regulations and laws applying to foreign researchers and state whether you obtained the required permits and approvals.Please address this in your ethics statement in both the manuscript and submission information.4. We note that Figure 1 in your submission contains map images which may be copyrighted.All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For these reasons, we cannot publish previously copyrighted maps or satellite images created using proprietary data, such as Google software (Google Maps, Street View, and Earth). For more information, see our copyright guidelines: http://journals.plos.org/plosone/s/licenses-and-copyright.We require you to either (a) present written permission from the copyright holder to publish this figure specifically under the CC BY 4.0 license, or (b) remove the figure from your submission:a. You may seek permission from the original copyright holder of Figure 1 to publish the content specifically under the CC BY 4.0 license.We recommend that you contact the original copyright holder with the Content Permission Form (http://journals.plos.org/plosone/s/file?id=7c09/content-permission-form.pdf) and the following text:“I request permission for the open-access journal PLOS ONE to publish XXX under the Creative Commons Attribution License (CCAL) CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). Please be aware that this license allows unrestricted use and distribution, even commercially, by third parties. Please reply and provide explicit written permission to publish XXX under a CC BY license and complete the attached form.”Please upload the completed Content Permission Form or other proof of granted permissions as an "Other" file with your submission.In the figure caption of the copyrighted figure, please include the following text: “Reprinted from [ref] under a CC BY license, with permission from [name of publisher], original copyright [original copyright year].”b. If you are unable to obtain permission from the original copyright holder to publish this figure under the CC BY 4.0 license or if the copyright holder’s requirements are incompatible with the CC BY 4.0 license, please either i) remove the figure or ii) supply a replacement figure that complies with the CC BY 4.0 license. Please check copyright information on all replacement figures and update the figure caption with source information. If applicable, please specify in the figure caption text when a figure is similar but not identical to the original image and is therefore for illustrative purposes only.The following resources for replacing copyrighted map figures may be helpful:USGS National Map Viewer (public domain): http://viewer.nationalmap.gov/viewer/The Gateway to Astronaut Photography of Earth (public domain): http://eol.jsc.nasa.gov/sseop/clickmap/Maps at the CIA (public domain): https://www.cia.gov/library/publications/the-world-factbook/index.html and https://www.cia.gov/library/publications/cia-maps-publications/index.htmlNASA Earth Observatory (public domain): http://earthobservatory.nasa.gov/Landsat: http://landsat.visibleearth.nasa.gov/USGS EROS (Earth Resources Observatory and Science (EROS) Center) (public domain): http://eros.usgs.gov/#Natural Earth (public domain): http://www.naturalearthdata.com/5. In your manuscript, please provide additional information regarding the specimens used in your study.Ensure that you have reported specimen numbers and complete repository information, including museum name and geographic location.If permits were required, please ensure that you have provided details for all permits that were obtained, including the full name of the issuing authority, and add the following statement:'All necessary permits were obtained for the described study, which complied with all relevant regulations.'If no permits were required, please include the following statement:'No permits were required for the described study, which complied with all relevant regulations.'For more information on PLOS ONE's requirements for paleontology and archaeology research, see https://journals.plos.org/plosone/s/submission-guidelines#loc-paleontology-and-archaeology-research.6. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.[Note: HTML markup is below. Please do not edit.]Reviewers' comments:Reviewer's Responses to QuestionsComments to the Author1. Is the manuscript technically sound, and do the data support the conclusions?The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.Reviewer #1: YesReviewer #2: Yes**********2. Has the statistical analysis been performed appropriately and rigorously?Reviewer #1: YesReviewer #2: Yes**********3. Have the authors made all data underlying the findings in their manuscript fully available?The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.Reviewer #1: YesReviewer #2: Yes**********4. Is the manuscript presented in an intelligible fashion and written in standard English?PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.Reviewer #1: YesReviewer #2: Yes**********5. Review Comments to the AuthorPlease use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)Reviewer #1: The work of Piazza et al. 2020 is a thorough and comprehensive look at a locality which represents a better case scenario for marine life trying to survive through the Toarcian extinction event. Rather than the organic-rich shales typically deposited during the severe anoxia of the period, the locality described by the authors exhibited subtler patterns of environmental change, yet still overall displayed a dramatic turnover and transition of the benthic macrofauna known from the site.The work represents a useful and rigorous case study for researchers expanding our global understanding of biotic change across the Toarcian, an interval of warming and anoxia which caused global disruptions to the carbonate factory paralleling the present crisis our oceans face. I don't have many comments, but present a few opportunities to expand and search for additional lessons from this obviously fruitful fossil locality. No new data would be needed, but these represent a few additional analyses which could bring to light other lessons from your impressive dataset!68: Here you mention the application of this study to what could be termed conservation paleobiology (using dynamics during past change to better understand our present time), but there is not much follow up. One way I could think of to tie in your results to the present crisis would be to provide a little more detail in the discussion or conclusion on the types of modern biomes the fauna of your site best represent, and how the stressors of the Toarcian are relevant. It could potentially go around line 515 where you start talking about modern community studies. Regarding temperature, it doesn't pass muster for a formal paleotemperature reconstruction but you could use a seawater d18O value of -1 per mil (assuming ice-free world) and as a best-guess taking off point of how the temperature conditions at your site relate to places in the modern day, before and during the extinction event.260: Many benthic ecologists would hesitate to group true epifaunal taxa (recliners, cementers, etc) with semi-infaunal taxa. Semi-infaunal bivalves often have very distinct needs in terms of substrate grain size, food type, environmental energy level from reclining bivalves, at least calling on my experience with glycymerid bivalves, which can highly stenotopic in their substrate preference. Was this grouping made to improve statistical power? Could another analysis with them separated be run, and reported if it said anything different?I'd be very interested in how dynamics in occurrence semi-infaunal bivalves relate to pedically attached brachiopods specifically, and also how reclining bivalves relate to the unattached brachiopods.444: It is perhaps unsurprising that you find that temperature/d18O is the main determinant of faunal change while d13C's influence is messy/absent. I went back to the supplement of the Ullman et al 2020 paper and saw in Supplemental figure S9 that there is an offset between the brachiopod and bivalve results. While combining them is of course the best way to get a large sample size to reconstruct a detailed record of the excursion, perhaps the fact that the excursion itself is poorly correlated to faunal change is not a big shock, as it is a global signal that might not have a huge amount of relevance to the more granular local faunal composition in the way that temperature does.Could you "un-combine" the record and compare correlation of brachiopod d13C to brachiopod occurrence and vice versa for bivalves? As they can have dietary differences and I'd wonder if there were any dynamics regarding tiering over the interval. You could even see if the two records have an interaction, or subset the brachiopod d13C record by ecology, considering it looks like you have a very large sample size for them. Could be an additional NMDS to run which might help elucidate any ecological trade offs or succession between the groups across the boundary.500: The energy may have remained relatively low throughout the site, but were there any changes in faunal composition associated with the packstones/rudstones? Indicating possible ecological disruptions from the intermittent storms which required a mini recovery interval? I understand if the resolution you're working with is not enough to make this characterization but might be relevant to those interested in how subtle substrate changes influence the community composition through time. Just bringing up because I had gone back to Ullmann's recent Scientific Reports "Warm Afterglow" paper which described "rhythmic alternations of marlstones and partly argillaceous limestones. The latter primarily comprise mudstones, wackestones, and floatstones." Any faunal associations with these alternations?I see a couple lines below you say diversity trends are "not consistent" with sea level changes, but could you elaborate on how you made this determination?Regards and good work,Dan KillamPostdoctoral Researcher, Biosphere 2https://dantheclamman.blogReviewer #2: Review of manuscript PONE D-20-11804Entitled “Temperature affects faunal dynamics of benthic invertebrate assemblages across the Toarcian Oceanic Anoxic Event in the Iberian Basin (Spain) by Piazza et al.The manuscript for the most part is well written, and provides a lot of detail about one of the many oxygen minimum zones in the Toarcian. The isotope trends are clear markers of this particular event in the Iberian Basin.I have a few and minor comments that the authors should address before acceptance:Line 146 d13C -> carbon213 i.e. -> such as310 numbers -> …in terms of species and number of specimens.314 On the contrary -> In contrast; whereas not while, sentence is not clear318 , while -> whereas species abundance distributions are more even in the post TOAE326 while ? if at the SAME time, then ok , if not then use ‘whereas’338 IBID (326)389 delete ‘down’427-430 long sentence -please re-write445 …d18O ‘values’; l. 448 …d13C ‘values’ in isotope terminology the d13C is treated as an adjective and thus needs a ‘noun’IBID 486, etc457 …d18O and d13C ‘values’ of what??444, 474 non—significant -> not significant524-527 sentence unclear, 527-528 ibid567 i.e. -> …with cooler temperatures indicating the end..”L 442 if all d18O and d13C values are based on brachiopods, of short mention in Table 3 and text would be nice.**********6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.If you choose “no”, your identity will remain anonymous but your review may still be made public.Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.Reviewer #1: Yes: Daniel KillamReviewer #2: Yes: Uwe Brand[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.6 Oct 2020POINT-BY-POINT REBUTTAL LETTERJournal RequirementsWhen submitting your revision, we need you to address these additional requirements.1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdfWe have checked the templates and made the few adjustments to meet the journal requirements. Moreover, an additional private email address of the corresponding author V.P. on the title page has been added to ensure continuity in communication in cases of future inquiries about the manuscript. The reason is that the given institutional address will be deactivated in the near future.2. In your Methods section, please provide additional information regarding the permits you obtained for the work. Please ensure you have included the full name of the authority that approved the field site access and, if no permits were required, a brief statement explaining why.More details on the authority that approved of our fieldwork and issued the permit have been given at the end of the “Material and methods” section.3. During our internal checks, the in-house editorial staff noted that you conducted research or obtained samples in another country. Please check the relevant national regulations and laws applying to foreign researchers and state whether you obtained the required permits and approvals. Please address this in your ethics statement in both the manuscript and submission information.The original sentence provided in the first version of the manuscript has been modified with the inclusion of the statement provided under Requirement #5.4. We note that Figure 1 in your submission contains map images which may be copyrighted. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For these reasons, we cannot publish previously copyrighted maps or satellite images created using proprietary data, such as Google software (Google Maps, Street View, and Earth). For more information, see our copyright guidelines: http://journals.plos.org/plosone/s/licenses-and-copyright. We require you to either (a) present written permission from the copyright holder to publish this figure specifically under the CC BY 4.0 license, or (b) remove the figure from your submission:a. You may seek permission from the original copyright holder of Figure 1 to publish the content specifically under the CC BY 4.0 license. We recommend that you contact the original copyright holder with the Content Permission Form (http://journals.plos.org/plosone/s/file?id=7c09/content-permission-form.pdf) and the following text:“I request permission for the open-access journal PLOS ONE to publish XXX under the Creative Commons Attribution License (CCAL) CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). Please be aware that this license allows unrestricted use and distribution, even commercially, by third parties. Please reply and provide explicit written permission to publish XXX under a CC BY license and complete the attached form.”Please upload the completed Content Permission Form or other proof of granted permissions as an "Other" file with your submission. In the figure caption of the copyrighted figure, please include the following text: “Reprinted from [ref] under a CC BY license, with permission from [name of publisher], original copyright [original copyright year].”b. If you are unable to obtain permission from the original copyright holder to publish this figure under the CC BY 4.0 license or if the copyright holder’s requirements are incompatible with the CC BY 4.0 license, please either i) remove the figure or ii) supply a replacement figure that complies with the CC BY 4.0 license. Please check copyright information on all replacement figures and update the figure caption with source information. If applicable, please specify in the figure caption text when a figure is similar but not identical to the original image and is therefore for illustrative purposes only. The following resources for replacing copyrighted map figures may be helpful:USGS National Map Viewer (public domain): http://viewer.nationalmap.gov/viewer/The Gateway to Astronaut Photography of Earth (public domain): http://eol.jsc.nasa.gov/sseop/clickmap/Maps at the CIA (public domain): https://www.cia.gov/library/publications/the-world-factbook/index.html and https://www.cia.gov/library/publications/cia-maps-publications/index.htmlNASA Earth Observatory (public domain): http://earthobservatory.nasa.gov/Landsat: http://landsat.visibleearth.nasa.gov/USGS EROS (Earth Resources Observatory and Science (EROS) Center) (public domain): http://eros.usgs.gov/#Natural Earth (public domain): http://www.naturalearthdata.com/The required permission was requested from the original copyright holder through the RightsLink service. The permission is submitted together with the revised manuscript. The caption of the copyrighted figure was adjusted with the addition of the relevant information.5. In your manuscript, please provide additional information regarding the specimens used in your study. Ensure that you have reported specimen numbers and complete repository information, including museum name and geographic location.If permits were required, please ensure that you have provided details for all permits that were obtained, including the full name of the issuing authority, and add the following statement:'All necessary permits were obtained for the described study, which complied with all relevant regulations.'If no permits were required, please include the following statement:'No permits were required for the described study, which complied with all relevant regulations.'For more information on PLOS ONE's requirements for paleontology and archaeology research, see https://journals.plos.org/plosone/s/submission-guidelines#loc-paleontology-and-archaeology-research.Specimen numbers and complete repository information are given at the end of the “Material and methods” section of the manuscript. Details about fieldwork permits are also provided at the end of the same section (see answers to Requirements #2 and #3), with the addition of the statement provided above.6. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.We understand this, and will provide the DOI number as soon as possible upon acceptance of the manuscript.Comments to the AuthorReviewer #1The work of Piazza et al. 2020 is a thorough and comprehensive look at a locality which represents a better case scenario for marine life trying to survive through the Toarcian extinction event. Rather than the organic-rich shales typically deposited during the severe anoxia of the period, the locality described by the authors exhibited subtler patterns of environmental change, yet still overall displayed a dramatic turnover and transition of the benthic macrofauna known from the site.The work represents a useful and rigorous case study for researchers expanding our global understanding of biotic change across the Toarcian, an interval of warming and anoxia which caused global disruptions to the carbonate factory paralleling the present crisis our oceans face. I don't have many comments, but present a few opportunities to expand and search for additional lessons from this obviously fruitful fossil locality. No new data would be needed, but these represent a few additional analyses which could bring to light other lessons from your impressive dataset!We thank the reviewer for his kind comments and useful insights. His suggestions touch indeed interesting aspects!68: Here you mention the application of this study to what could be termed conservation paleobiology (using dynamics during past change to better understand our present time), but there is not much follow up. One way I could think of to tie in your results to the present crisis would be to provide a little more detail in the discussion or conclusion on the types of modern biomes the fauna of your site best represent, and how the stressors of the Toarcian are relevant. It could potentially go around line 515 where you start talking about modern community studies. Regarding temperature, it doesn't pass muster for a formal paleotemperature reconstruction but you could use a seawater d18O value of -1 per mil (assuming ice-free world) and as a best-guess taking off point of how the temperature conditions at your site relate to places in the modern day, before and during the extinction event.We thank the reviewer for this insight. We have added a new paragraph at the end of the discussion comparing early Toarcian and modern faunas, and how the studied event is relevant to modern climate studies and future warming predictions.Regarding temperature, we understand that the reviewer’s comment concerns absolute temperature values. Those values are given in our related paper (Ullmann et al. 2020), which includes a discussion of the underlying assumptions. In Figure 3 we kept the more conservative presentation of relative changes in temperature. However, to accommodate the reviewer’s suggestion, we provide absolute temperature estimates in the newly added paragraph mentioned above with reference to the related paper Ullmann et al. (2020). This decision is based on our previous peer-review experience with using absolute temperatures rather than relative temperature change in Ullmann et al. (2020) and Piazza et al. (2020), where reviewers took an adverse stance against absolute temperature estimates, stressing the uncertainties of making inferences on absolute temperatures. There is ongoing debate which of the growing number of published oxygen isotope thermometer equations is the most suitable. Also, absolute temperature reconstructions require assumptions of seawater composition, which are possible, but somewhat uncertain. Given the above limitations, absolute temperature estimates likely have an uncertainty of at least 4°C. Most of this uncertainty is cancelled out when computing temperature change with respect to a reference datum, e.g., the pre-TOAE average. The sensitivity of oxygen isotope thermometers (permil change per °C) is similar and changes only slightly or not at all as a function of temperature in the range of 10 to 30 °C so that the choice of thermometer and the exact position on the temperature scale loose importance. Therefore, the use of relative temperatures allows to assess changing ambient conditions at one site with high confidence and does prevent the (understandable) temptation of comparing absolute temperatures to possible modern analogues.260: Many benthic ecologists would hesitate to group true epifaunal taxa (recliners, cementers, etc) with semi-infaunal taxa. Semi-infaunal bivalves often have very distinct needs in terms of substrate grain size, food type, environmental energy level from reclining bivalves, at least calling on my experience with glycymerid bivalves, which can highly stenotopic in their substrate preference. Was this grouping made to improve statistical power? Could another analysis with them separated be run, and reported if it said anything different?I'd be very interested in how dynamics in occurrence semi-infaunal bivalves relate to pedically attached brachiopods specifically, and also how reclining bivalves relate to the unattached brachiopods.We understand and agree with this observation, and we thank the reviewer for pointing this out. Indeed, the main reason for lumping semi-infaunal and fully epifaunal taxa, both being exposed to potential predators at the sea floor, was to increase statistical power. Only ca. 10% of the specimens belong to the semi-infauna (see MOL 7 in Fig. 7) – too few to be considered on its own. To see the influence of the semi-infauna we followed the reviewer’s suggestion and analysed the two life-habit categories separately. Overall, the results obtained for epifauna + semi-infauna and for true epifauna are congruent (i.e. no significant changes). To address the concerns of the reviewer, we now present the results for the true epifauna instead of the lumped epifauna + semi-infauna. The text, Table 2 and Fig. 9 have been modified accordingly.As for semi-infaunal bivalves vs. brachiopods, we feel there is not much to add in this respect. The important information was already presented in S2 Fig. To make this more evident we now moved previous S2 Fig to the main manuscript where it became the new Fig. 7. It shows the relative abundance change of all MOLs including semi-infaunal bivalves (MOL 7) and pedically attached brachiopods (MOL 2). The only other realised life habit in brachiopods is represented by the shallow infaunal attached Lingularia sp. (MOL 14) which is very rare. No unattached brachiopods are present in the studied fauna as the reviewer seems to imply.444: It is perhaps unsurprising that you find that temperature/ δ18O is the main determinant of faunal change while δ13C's influence is messy/absent. I went back to the supplement of the Ullman et al 2020 paper and saw in Supplemental figure S9 that there is an offset between the brachiopod and bivalve results. While combining them is of course the best way to get a large sample size to reconstruct a detailed record of the excursion, perhaps the fact that the excursion itself is poorly correlated to faunal change is not a big shock, as it is a global signal that might not have a huge amount of relevance to the more granular local faunal composition in the way that temperature does.Could you "un-combine" the record and compare correlation of brachiopod δ13C to brachiopod occurrence and vice versa for bivalves? As they can have dietary differences and I'd wonder if there were any dynamics regarding tiering over the interval. You could even see if the two records have an interaction, or subset the brachiopod δ13C record by ecology, considering it looks like you have a very large sample size for them. Could be an additional NMDS to run which might help elucidate any ecological trade offs or succession between the groups across the boundary.The reviewer has made a series of quite interesting points. "Uncombining" the carbon isotope records of bivalves and brachiopods can easily be done as all the raw data are available and are listed in the supplements of Ullmann et al. (2020), but we don't think that it would yield significant further insights. Bivalves are invariably enriched in 13C versus coeval brachiopods and most enrichment factors agree with the average within uncertainty. Only in the range of 5 to 7 m height in the section the values are outside 95 % confidence, but this seems to be a part of the section where nothing really noteworthy changes regarding assemblages. Even though overall the dataset is large, taking the dataset apart might lead to over-interpreting some of the results because we are dealing with signals much smaller than 0.5 permil, if they are present at all. We feel that the safest approach - given the heterogeneity of macrofossil datasets - is stating that brachiopod and bivalve carbon isotope records run parallel with an approximately 0.6 permil offset towards higher δ13C in bivalves.It may be useful to revisit the observation of bivalve 13C enrichment (but missing 18O offset) and to try to find a satisfactory explanation rather than keeping it at the pragmatic solution we employed then, but so far it seems to us that we (the community in general that is) don't really know enough about exactly how these different organisms form shell material and what this may mean for isotopic signatures.Regarding inter-species differences, there are no clear differences between any of the rhynchonellid genera (they are certainly not larger than between individuals of the same species and this point has been explored in detail in the supplements of Ullmann et al., 2020). We would hesitate making environmental inferences from comparison to terebratulid data even if we had sufficient data to do so because terebratulid vital effects make this too prone to biased outcomes. Bivalve data are Gryphaea-only, so it would be impossible to look into finer detail there.500: The energy may have remained relatively low throughout the site, but were there any changes in faunal composition associated with the packstones/rudstones? Indicating possible ecological disruptions from the intermittent storms which required a mini recovery interval? I understand if the resolution you're working with is not enough to make this characterization but might be relevant to those interested in how subtle substrate changes influence the community composition through time. Just bringing up because I had gone back to Ullmann's recent Scientific Reports "Warm Afterglow" paper which described "rhythmic alternations of marlstones and partly argillaceous limestones. The latter primarily comprise mudstones, wackestones, and floatstones." Any faunal associations with these alternations? I see a couple lines below you say diversity trends are "not consistent" with sea level changes, but could you elaborate on how you made this determination?This is a valid point. We are aware of the influence that substrate changes can have on faunal composition. Yet, our samples are derived from fairly homogenous lithologies (wackestones and floatstones with just one sample frpm a rudstone. We now made this clear in the manuscript text. We also performed again the cluster and NMDS analyses on both taxonomical and functional composition of the studied faunas, including lithology as a grouping factor to investigate whether faunal associations related to specific lithologies could be observed. It is quite apparent that lithological differences did not have any larger effect on faunal composition. These new results are briefly mentioned in the text and the clusters in Supplementary Figure 1 were adjusted to include lithological differences.Regarding diversity trends in relation to sea level change, we have made adjustments to make clearer the basis of our statement.Reviewer #2The manuscript for the most part is well written, and provides a lot of detail about one of the many oxygen minimum zones in the Toarcian. The isotope trends are clear markers of this particular event in the Iberian Basin. I have a few and minor comments that the authors should address before acceptance:We thank the reviewer for his positive assessment. As illustrated in detail below, we have addressed most of the comments, and justified the few ones we decided not to incorporate in the revised manuscript.Line 146 δ13C -> carbonDone.213 i.e. -> such asThe Reviewer suggested substitution of “i.e.” with “such as”. We feel that it would imply the explanation for MOL given here as an example, while we provide a definition. We use then “–“ instead.310 numbers -> …in terms of species and number of specimens.Done.314 On the contrary -> In contrast; whereas not while, sentence is not clear.Done.318 , while -> whereas species abundance distributions are more even in the post TOAEDone.326 while ? if at the SAME time, then ok , if not then use ‘whereas’As in that sentence we are still discussing the same time interval (the pre-TOAE), we keep “while” instead of “whereas”.338 IBID (326)In this instance we changed “while” into “whereas”.389 delete ‘down’Done.427-430 long sentence -please re-writeWe have rearranged the sentence.445 … δ18O ‘values’; l. 448 … δ13C ‘values’ in isotope terminology the d13C is treated as an adjective and thus needs a ‘noun’We have adjusted it here, and also when needed elsewhere in the text.IBID 486, etcWe did not see any mention of isotopes at line 486 of the originally submitted manuscript. This comment probably referred to line 468 and we adjusted the text here and elsewhere in the manuscript.457 … δ18O and d δ13C ‘values’ of what??Adjusted.444, 474 non—significant -> not significantDone.524-527 sentence unclear, 527-528 ibidSentences rearranged to make the content clear.567 i.e. -> …with cooler temperatures indicating the end..”Done.L 442 if all δ18O and δ13C values are based on brachiopods, of short mention in Table 3 and text would be nice.As stated in the methods, isotope data are derived from both oysters and brachiopods. We added this information again in the caption of Table 3, and feel therefore that no further change in the text is required.Submitted filename: Response to Reviewers.pdfClick here for additional data file.2 Nov 2020Ocean warming affected faunal dynamics of benthic invertebrate assemblages across the Toarcian Oceanic Anoxic Event in the Iberian Basin (Spain)PONE-D-20-11804R1Dear Dr. Piazza,We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.Kind regards,David P. Gillikin, Ph.D.Academic EditorPLOS ONE17 Nov 2020PONE-D-20-11804R1Ocean warming affected faunal dynamics of benthic invertebrate assemblages across the Toarcian Oceanic Anoxic Event in the Iberian Basin (Spain)Dear Dr. Piazza:I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.If we can help with anything else, please email us at plosone@plos.org.Thank you for submitting your work to PLOS ONE and supporting open access.Kind regards,PLOS ONE Editorial Office Staffon behalf ofDr David P. GillikinAcademic EditorPLOS ONE
Fig 1
Fig 2
Fig 3
Fig 5
Fig 4
Fig 6
Fig 7
Fig 8
Fig 9