The environmental factors limiting the distribution of shallow-water terebratulid brachiopods

Brachiopods were the most successful benthic marine animals during the Paleozoic (Thayer 1986), with about 400 species remaining today (Emig et al. 2013). In Recent seas, terebratulids are by far the most successful of all brachiopod clades (Lee 2008). Like other representatives of the subfamily Terebratulinae, such as Pliothyrina and Maltaia, some species of “Terebratula” appear to stand as an exception to the general progressive trend of brachiopod retreat to deep and/or cryptic habitats after the Permian–Triassic mass extinction (Tomašových 2006). This is because Terebratulinae were numerically abundant above the storm-weather wave base (Kroh et al. 2003; Gramigna et al. 2008; Pervesler et al. 2011), much like other Terebratulida can dominate in present-day, shallow-water hard-bottom environments (Logan and Noble 1971; Richardson 1981; Försterra et al. 2008; Tomašových 2008). European Terebratulinae, however, went extinct in the Calabrian (approximately during or shortly after the Jaramillo Subchron), not surviving beyond the Sicilian (D’Alessandro et al. 2004; Taddei Ruggiero and Taddei 2006; La Perna and Vazzana 2016; Crippa et al. 2019), with the last species being Terebratula terebratula and Terebratula scillae. Although Neogene and Pleistocene Terebratula had to cope with specialized predators, niche competitors (bivalves), and bulldozing and grazing marine organisms (Thayer 1983), it thrived mainly in warm-temperate shallow-water detritic-bottom habitats in the Paratethys and the Mediterranean realms (Pedley 1976; Bitner and Pisera 2000; Reolid et al. 2012). There are also reports from subtropical environments, where species of Terebratula infilled cavities in Porites–Tarbellastraea reefs during the Tortonian (Barbera et al. 1995) or were present on the distal slope of Porites fringing reefs during the lower Messinian (Llompart and Calzada 1982; Obrador et al. 1983; Videt 2003). Little is known, however, about the paleoecology of Terebratula (Benigni and Robba 1990; Pavia and Zunino 2008). Narrowing down the factors that led European Terebratulinae to extinction calls for improving our knowledge about the habitats and the ecological niche occupied by these terebratulids in their prime. This study is designed to help fill this gap by evaluating the environmental conditions in the subenvironments where the Terebratula populations were at their optimum in the Pliocene outcrops of the Águilas Basin (SE Spain), as well as the paleoecological boundary conditions onshore and offshore of this optimum. More specifically, we estimate environmental tolerances of this brachiopod with respect to water depth, bottom oxygenation, bottom temperature, productivity, and organic enrichment. This study area is among the most suitable for this purpose, because abundant and almost monospecific shell beds of Terebratula occur here cyclically in a mixed temperate carbonate–siliciclastic system (García-Ramos and Zuschin 2019).

The Águilas Basin (SE Spain) is located in the southwestern inner sector of the tectonic Águilas Arc (Fig. 1A,B), in the Internal Betic Zone (Coppier et al. 1989). The Águilas Basin was a small embayment (about 14 km wide) during the early Pliocene (Fig. 1C) (García-Ramos and Zuschin 2019).

The studied brachiopods belong to a sequence of late Zanclean age (MPl3 planktic foraminiferal biozone of the Mediterranean Pliocene; for the biozonation, see Iaccarino et al. [2007] and Corbí and Soria [2016]). The present study also yielded scarce specimens of the nannoplankton species Reticulofenestra cf. cisnerosi at the base of the sequence. This occurrence, together with the absence of Discoaster tamalis and the association with Globorotalia puncticulata and Globorotalia margaritae, constrains the age of the base of the sequence to the older part of Subchron C3n.1r (between 4.52 and 4.42 Ma), according to the biostratigraphic scheme of Lancis et al. (2015).

The studied strata define originally inclined units that were deposited on a slope (termed “clinobeds”), which were interpreted to reflect high-frequency and low-amplitude relative sea-level changes (García-Ramos and Zuschin 2019). Sedimentologically, these clinobedded units are consistent with sand-prone subaqueous delta-scale clinoforms (sensu Patruno et al. 2015), which developed entirely below sea level. The most proximal deposits are bioturbated coarse sands, followed distally by rhodolith-rich finer-grained sediments (rhodalgal facies). Both facies correspond to the topset of the clinoforms (Fig. 2A). Following a biological benthic zonation (i.e., Gili et al. 2014), these environments were interpreted as mediolittoral to infralittoral (foreshore to upper shoreface) and lower infralittoral to upper circalittoral (lower shoreface) (García-Ramos and Zuschin 2019). The rhodolith debris disappears gradually toward the basin, giving way to fine-grained sands rich in the polychaete Ditrupa arietina along with benthic and planktic foraminifera. The disappearance of the rhodoliths distally can be interpreted as a transition zone from the upper circalittoral to the lower circalittoral (Basso 1998; Cameron and Askew 2011; Gili et al. 2014). This facies transition between rhodalgal and Ditrupa-rich deposits was also interpreted as the offshore transition zone (OTZ) based on sedimentological evidence; it coincides with the topset–foreset transition (i.e., upper rollover) of the clinoforms (García-Ramos and Zuschin 2019). The Ditrupa facies in this transition zone to the lower circalittoral corresponds to the foreset of the clinoforms (Fig. 2A). The most distal facies in our area is fine-grained muddy sands with the characteristic pectinid Costellamussiopecten cristatum occurring as dispersed, disarticulated shells (Fig. 2A, bottomset of the clinoforms). These deposits can be interpreted as formed in upper offshore environments, also within the transition zone to the lower circalittoral.

In the above facies, 5- to 20-cm-thick monospecific pavements of Terebratula calabra (Figs. 2B and 3A,B) are interspersed cyclically (Fig. 2A). Locally the pavements also yield rare specimens of the rhynchonellid Aphelesia bipartita. As a rule, the greatest Terebratula concentrations occur in the transition between the foreset and the bottomset (the lower rollover zone) (Fig. 2A,B), where the Terebratula pavements show loosely to densely packed biofabrics (Figs. 2B and 3B). The density of Terebratula decreases toward more proximal (foreset) and more distal positions (bottomset) (Figs. 2A and 3C) in the depositional profile. The rhodalgal facies (lower shoreface) also contains rare, isolated specimens of Terebratula (Fig. 3D,E). The shells in the pavements are mostly complete, articulated, and minimally encrusted by epizoans or bioeroded (recording rare traces of Podichnus obliquus produced by conspecifics). These concentrations were interpreted as parautochthonous assemblages in relatively shallow-water but offshore (circalittoral) environments during stages of reduced sedimentation rates (for details, see García-Ramos and Zuschin 2019). A recognizable 1.5- to 2-m-thick bed dominated by Terebratula in the study area shows loose to densely packed biofabrics (Fig. 3F–J). The taxonomic richness in this bed is distinctly higher than in the pavements. In this bed, the “Terebratula biostrome,” highly altered shells occur admixed with clusters of well-preserved terebratulids, some of which can be interpreted as biological clumps (sensu Kidwell et al. 1986). This suggests that the terebratulids were autochthonous (García-Ramos and Zuschin 2019). Many shells in the biostrome are bioeroded by clionaid sponges (Fig. 3K), although other traces also occur in low abundance (Fig. 3L–R). The rasping trace Gnathichnus pentax (Fig. 3L) is represented, whereas the rasping trace Radulichnus isp. is absent (García-Ramos and Zuschin 2019). Molinu et al. (2013), who studied microbioerosion traces affecting Terebratula specimens from the biostrome at the Cañada Brusca E-2 section, reported that the traces produced by fungi were dominant, although the microendolith trace Rhopalia clavigera (the product of chlorophytes) is also represented. The Terebratula biostrome, compared with the pavements, yields additional brachiopod taxa: the craniid Novocrania anomala (Fig. 3S), often encrusting disarticulated shells of Terebratula, is relatively common. In contrast, Megathiris detruncata, Megerlia truncata, Terebratulina retusa, A. bipartita, and Maltaia moysae are rare or very rare.

From a morpho-sedimentary viewpoint, these clinoform systems have also been referred to as “infralittoral prograding wedges” (Hernández-Molina et al. 2000; Pomar and Tropeano 2001). Assuming the latter genetic model, the transition between the topset and the foreset of the clinoforms (i.e., the upper rollover) is coincident with the storm-weather wave base (lower shoreface–offshore transition) (Fig. 2A). A minimum water depth for this environment was probably about 25–30 m, in line with data from Recent examples from the Western Mediterranean (Hernández-Molina et al. 2000; Betzler et al. 2011) and the coralline algal assemblage (García-Ramos and Zuschin 2019). The vertical distance between the lower and upper rollover in the clinoform between clinothems 5 and 6 (García-Ramos and Zuschin 2019) is 29 m (Fig. 2B). This suggests a water depth of 54–59 m for the foreset–bottomset transition of the clinoforms (i.e., the lower rollover), without considering lithostatic compaction. During stages of low-amplitude (~15–20 m) relative sea-level rise pulses (García-Ramos and Zuschin 2019), the water depth of the lower rollover and bottomset was an estimated 70–80 m (Fig. 2A).

Brachiopod assemblages are well represented, and crop out cyclically, in the Águilas Basin (SE Spain) (Fig. 1A,B). To evaluate the range of paleoenvironmental conditions, we studied the micropaleontological content of 26 bulk samples of friable sediment from the Cabezo Alto (CA) section (Fig. 4). The CA section was chosen because it records a vertical succession continuously exposing the main facies in the study area (including two Terebratula pavements: samples CA10 and CA15) (Fig. 4). The CA section is a suitable template for comparison with other Terebratula outcrops from the same stratigraphic sequence (three additional samples from Terebratula pavements [TP1, TP2, and TP3] and one from the Terebratula biostrome [sample TB]). Comparison between CA and other outcrops helps assess the variability of the environmental conditions associated with the presence and absence of Terebratula in the studied localities.

The samples were screened to study benthic and planktic foraminifera and calcareous nannoplankton. Each of the three taxonomic groups was studied and analyzed as separate data sets. For each sample, >200 specimens of benthic foraminifera (>100 in the case of planktics) in the 125–500 μm fraction (Weinkauf and Milker 2018) were identified to species level whenever possible and counted. Only tests >50% complete, which included diagnostic features, were considered. For calcareous nannoplankton, smear slides were prepared using standard procedures and examined under the light microscope (cross and parallel nicols) at 1000× magnification. The occurrence of ascidian spicules was noted. Quantitative data were obtained by counting at least 300 specimens from each smear slide. A further 100 view squares were checked for important species to interpret the biostratigraphy and paleoecology of the calcareous nannoplankton. Among reticulofenestrids, we followed the general distinction based on size (Young 1998): Reticulofenestra minuta (<3 μm), Reticulofenestra haqii (3–5 μm), Reticulofenestra pseudoumbilica (5–7 μm), and R. pseudoumbilica (>7 μm).

Before subsequent analyses, differences in sample size were accounted for by rarefying the Q × R matrix with count abundance data of benthic foraminifera so that each sample contained 200 specimens (100 specimens in the planktic foraminiferal data set). This was accomplished using the function rrarefy of the package vegan (Oksanen et al. 2018). All subsequent analyses were conducted in the R statistical environment, v. 3.5 (R Development Core Team 2018). The resampling process was repeated 100 times, and the mean data set was used for further analysis. For benthic foraminifera, species with relative abundances below 3% were discarded. The variability of the environmental parameters is shown in box plots by comparing samples where Terebratula was abundant, rare, or formed the biostrome. Samples lacking Terebratula were segregated between bottomset and foreset subenvironments (both in the transition zone from upper to lower circalittoral) (Fig. 2A) based on facies analysis conducted by García-Ramos and Zuschin (2019).

We used a Q-mode nonmetrical multidimensional scaling (NMDS) as an ordination method to visualize in the Bray-Curtis multivariate space the position of the samples containing Terebratula along the environmental gradient encompassing bottomset and foreset subenvironments at the CA section. Differences in the composition of benthic foraminiferal assemblages among bottomset, foreset, and Terebratula samples were additionally evaluated with a test of permutational multivariate analysis of variance (PERMANOVA, function adonis in package vegan). The variation in composition of such assemblages was examined with a permutational analysis of multivariate dispersions (PERMDISP, function betadisper in package vegan). Facies occurring in shallower environments (e.g., the rhodalgal facies) were excluded from the analysis, because this facies consists of hardened rock that hampers the extraction of foraminiferal tests. The composition of benthic foraminifera in the samples containing Terebratula is shown in bar plots. The composition of the remaining samples can be checked in a two-way cluster analysis (Supplementary Fig. 1).

Additionally, we provide estimates for bathymetry using the transfer function proposed by Báldi and Hohenegger (2008). This approach was applied using benthic foraminiferal species with both non-overlapping and overlapping depth ranges. The depth ranges of benthic foraminifera were mainly compiled from Sgarrella and Moncharmont Zei (1993), Altenbach et al. (2003), Hohenegger (2005), Rasmussen (2005), Spezzaferri and Tamburini (2007), Sen Gupta et al. (2009), Phipps et al. (2010), and Milker and Schmiedl (2012). The same transfer function, weighted by the mean, was applied on the vertical distribution range of planktic foraminiferal species to estimate minimum water-column depth. The depth ranges of planktic foraminifera were compiled from Rebotim et al. (2017). For comparison, we also computed depth estimates based on the plankton/benthic ratio (P/B) using the regression function by van der Zwaan et al. (1990).

Changes in organic matter content are often coupled with reduced oxygen concentrations at the seafloor (Jorissen et al. 1995; Koho et al. 2008). An increase in organic matter is evaluated by analyzing changes in the proportional abundance of foraminiferal species assigned to five ecological groups (EG1 to EG5). These groups were proposed based on different degrees of opportunistic species’ response to varying levels of organic matter enrichment (Alve et al. 2016; Jorissen et al. 2018). We used only species whose proportional abundance was >3%, because ecological information on rare species is often not available. For visualization in box plots, percentages were recalculated after discarding rare species (Dominici et al. 2008). For some quantitatively important species (e.g., Planulina ariminensis) not listed by Alve et al. (2016) and Jorissen et al. (2018), a tentative attribution to a category of EG was attempted based on literature evidence (e.g., Rasmussen 2005). In this study, for each EG, we report the proportion of taxa, including unlisted species that were tentatively attributed to a corresponding group. We denoted this as the “aggregated ecological group” (AEG). To visualize the results, we used the conceptual TROX model (Jorissen et al. 1995), representing the Terebratula samples in a scheme adapted from Koho et al. (2008) by incorporating the main benthic foraminiferal species found in our samples.

The ecological and environmental distribution of extant species of planktic foraminifera is associated with levels of primary productivity and temperature (roughly warm-oligotrophic vs. cold-eutrophic species) (Hemleben et al. 1989; Sierro et al. 2003). The taxonomy generally follows the concepts of Kennett and Srinivasan (1983) with consideration of other sources (e.g., Poore and Berggren 1975; Malgrem and Kennett 1977; Darling et al. 2006; Schiebel and Hemleben 2017). Productivity based on planktic foraminiferal assemblages was assessed by comparing samples containing Terebratula with those where it was absent. The attribution of the planktic species to a category (warm-oligotrophic vs. cold-eutrophic) was based on information provided by Spezzaferri et al. (2002), Sierro et al. (2003), and Incarbona et al. (2013). We performed a modern analogue technique (MAT) with the packages analogue and analogueExtra (Simpson 2007) to estimate sea-surface temperature (SST) from the Águilas Zanclean samples. As a training data set, we used the modern North Atlantic planktic foraminiferal data set from Kucera et al. (2005), available at the Pangaea repository. The fossil species were standardized with the closest extant relatives using the morphogroups proposed by Serrano et al. (2007). The exception was that we included Globorotalia hirsuta as a proxy for G. margaritae (e.g., Globoturborotalita rubescens was taken as a proxy for Globoturborotalita gr. apertura; Globorotalia inflata for G. puncticulata; Globigerinoides ruber [lumping pink and white types] for Globigerinoides obliquus and Globigerinoides extremus). We computed an NMDS ordination of the modern data set and the Águilas samples, and we superimposed the calculated SST isotherms on the ordination plot using the function ordisurf from vegan. The variance of the Kucera et al. (2005) data set explained by SST was calculated with canonical correspondence analysis (CCA).

To evaluate bottom-water temperatures, we selected 12 specimens of Terebratula calabra collected from different stratigraphic intervals and shell beds from the Águilas Basin for oxygen isotope analysis. The specimens were embedded in resin and cut along a longitudinal axis to produce thin sections, which were subsequently screened for diagenetic alteration with the cathodoluminescence microscope Technosyn 8200 MK II at the GeoZentrum Nordbayern, University Erlangen–Nuremberg. Selected samples of nonluminescent Terebratula shells were polished, etched in 5% HCl solution for 15 seconds (Crippa et al. 2016), sputtered with gold, and photographed under a scanning electron microscope (Jeol JSM 6400) at the University of Vienna. Carbonate powders extracted with a microdrill from the secondary layer somewhat posterior to the middle part of the shell were reacted with 100% phosphoric acid at 70°C using a Gasbench II connected to a Thermo Fisher Delta V Plus mass spectrometer at the GeoZentrum Nordbayern, University Erlangen–Nuremberg. All values are reported in per mil relative to VPDB. Reproducibility and accuracy were monitored by replicate analysis of laboratory standards calibrated by assigning δ¹³C values of + 1.95‰ to NBS19 and −47.3‰ to IAEA-CO9 and δ¹⁸O values of −2.20‰ to NBS19 and −23.2‰ to NBS18. Reproducibility for δ¹³C and δ¹⁸O was ± 0.09 and ± 0.08 (1 SD), respectively. To estimate temperatures from oxygen isotopes, we used the equation from O’Neil et al. (1969). We assumed a seawater δ¹⁸O_sw = +1.5‰ VSMOW for the western Mediterranean in the late Zanclean, based on Recent seawater composition in the eastern Mediterranean (Schmidt 1999; Rohling 2013). Sample S104 was identified as an outlier and was excluded from analysis.

The Q-mode NMDS ordination plot (Fig. 5A) shows that the samples with Terebratula occupy a transitional position between Terebratula-barren bottomset and foreset samples, as was already observed from lateral facies variations in the field (Fig. 2B). The exceptions are CA10 (a bottomset sample with rare Terebratula) and TB (Terebratula biostrome), which shows affinity with the samples at the base of the CA section. PERMANOVA finds a significant difference among bottomset, foreset, and Terebratula samples, but only ~19% of the total variance is explained by these groups. Permutation tests applied on PERMDISP find a significant although moderate compositional variation between the bottomset and foreset samples (permutest: F = 4.77, df = 2, p = 0.022), whereas a post hoc Tukey’s honest significant difference test for multiple comparisons among groups indicates a significant compositional difference only between bottomset and foreset samples (p = 0.013). The benthic foraminiferal composition of the samples containing Terebratula is characterized, in general, by the dominance of Bolivina gr. dilatata and Cassidulina carinata, followed by cibicidids and asterigerinids (Fig. 6). The TB sample is dominated by cibicidids. The equation of Báldi and Hohenegger (2008) yielded depth estimates of ~50 to 90 m for the Terebratula samples when including non-overlapping shallow-water species and the vertical distribution range of planktic foraminifera (Fig. 5B). These values are consistent with estimates based on sedimentological models. Removing non-overlapping shallow-water species yielded bathymetric estimates of ~90 to 125 m. The regression function by van der Zwaan et al. (1990) returned unrealistic depth estimates for the Terebratula samples (50–200 m), most probably because of onshore transportation of planktic foraminifera by currents (Murray 1976).

The box plots show that samples with Terebratula contain a slightly higher proportion of benthic foraminifera belonging to the EG of “sensitive species” (AEG1) compared with samples where Terebratula is absent (Fig. 7). The Terebratula-containing samples vary between ~35% to ~60% of AEG1 (in the TB sample), the latter being the maximum value in the data set. All these samples fall within the range shown by samples devoid of Terebratula with regard to the proportion of benthic species belonging to the EGs of “indifferent species” (AEG2), “tolerant” (AEG3), and “second-order opportunists” (EG4). Most Terebratula samples, however, display low proportions (~10%) of benthic species of EG4 (Fig. 7). The exception is sample CA10, a bottomset sample, which contains ~20% of EG4 species and rare Terebratula. In the whole data set, “first-order opportunists” (EG5) were absent.

The Terebratula samples fall within the ranges of cold-eutrophic and warm-oligotrophic species contained in the Terebratula-barren samples (Fig. 8A). Sample CA10 displays the maximum proportion of cold-eutrophic species (~60%), whereas TB yields a relatively low proportion (~30%) (Fig. 8A). Overall, the Terebratula samples are dominated by warm-oligotrophic species, except for sample CA10 (Fig. 8A). The NMDS ordination of the nannoplankton (Fig. 8B) shows that the sample with rare Terebratula is associated with a Calcidiscus leptoporus–Coccolithus pelagicus assemblage, whereas the sample with abundant Terebratula is included in the Reticulofenestra small assemblage. In six additional Terebratula outcrops, nannoplankton were rare in four samples and absent in the other two. The samples with rare nannoplankton, however, contain the Calcidiscus leptoporus–Coccolithus pelagicus assemblage.

The MAT analysis using data from Kucera et al. (2005) shows that SST ranges from ~16°C (for sample CA10) to ~22°C (samples TP1, TP3, and Terebratula biostrome) (Fig. 9). The Terebratula samples CA15 and TP2 fall within the range of 17°C to 19°C. CCA estimates that SST explains 39% of the variance (permutation test: F = 540.46, df = 1, p = 0.001) in the Kucera et al. (2005) data set.

Regarding the oxygen isotopes, the Terebratula samples displayed a good preservation of the secondary layer fibers (Fig. 10A–C). The exception is sample S86, which shows cracks in the fibers. All analyzed Terebratula samples were nonluminescent, except for the punctae in some specimens (Fig. 10D–F). The volume of nonluminescent shell material is far greater than that of the sediment infilling the punctae. We therefore assume that the oxygen isotopic signal represents primary, diagenetically unaltered values. The δ¹⁸O_calc of the 12 Terebratula samples varies between −0.32 and 1.16, with a mean value of 0.71 (Table 1). Assuming that the seawater δ¹⁸O in the Águilas Basin during the late Zanclean was + 1.5‰ by analogy with the warmer Recent eastern Mediterranean (Schmidt 1999; Rohling 2013), then the δ¹⁸O_calc values translate into temperatures between 17.1°C and 23.8°C (mean: 19.1°C). Three measurements, done along the posterior–anterior axis of an isolated Terebratula shell from the foreset facies (sample CA28) (measurements S81, S102, S103), yield δ¹⁸O_calc values of 0.71, 0.91, and 0.57 VPDB. The calculated temperatures are 19.1°C, 18.2°C, and 19.7°C, respectively. A fourth measurement from the anterior part of the shell (S104) is considered to be an outlier (Table 1).

Terebratula calabra in the studied sequence of the Águilas Basin is very rare or absent in the rhodalgal facies (lower shoreface; ~25–30 m depth) (Figs. 2A and 4D,E) and in the deepest fine-grained sandy facies (offshore) (Figs. 2A, 4C, and 5). The maximum population density occurs close to (but notably below) the shoreface–offshore transition in fine-grained sands at estimated paleodepths of about 60 to 90 m (Figs. 2A,B, 4B, and 5B). The benthic foraminiferal assemblage of the Terebratula samples (Fig. 6) characterizes an offshore environment (e.g., Sgarrella and Moncharmont Zei 1993; Rasmussen 2005; Milker et al. 2009; Frezza et al. 2010; Mojtahid et al. 2010).

Modern Terebratulida bathymetrically equivalent to Terebratula mostly occur in shallow rocky habitats in bays and fjords (Supplementary Table 1). Terebratula species, instead, were abundant in soft sediments: coarse to muddy sands (e.g., Gaetani 1986; Barrier et al. 1987; Dominici 2001; Gramigna et al. 2008). Terebratula, like certain other brachiopods (Rudwick 1961; Richardson 1981; Llompart and Calzada 1982), was potentially attached to ascidians, whose spicules are pervasive in section CA. We did not quantify the percentage of Terebratula shells affected by pedicle attachment traces (P. obliquus), but this trace is rare. This matches previous reports on several species of Terebratula, in which 1% or less were affected by P. obliquus (Taddei Ruggiero and Bitner 2008). Other members of Terebratulinae are known to facultatively form clusters, such as Pliothyrina (Bell and Bell 1872; Rudwick 1961) and Liothyrella (Foster 1974; Richardson 1981; Peck et al. 1997). A cluster of Terebratula calabra from the Águilas Basin was illustrated by García-Ramos and Zuschin (2019).

Among the dominant benthic foraminifera in the Terebratula samples, a first group of species (B. gr. dilatata, Bulimina aculeata, C. carinata) suggests that Terebratula in this basin was subject to seasonal inputs of labile organic matter, possibly associated with episodes of dysoxia at the seafloor (Barmawidjaja et al. 1992; Fontanier et al. 2003; Langezaal et al. 2006; Abu-Zied et al. 2008; Mendes et al. 2012). This interpretation is reinforced by the presence of the Calcidiscus leptoporus–Coccolithus pelagicus nannoplankton assemblage in several Terebratula samples, suggesting productivity pulses triggered by coastal upwelling (Silva et al. 2009; Auer et al. 2014). In contrast, the “Reticulofenestra small” nannoplankton assemblage in sample CA15 points to an opportunistic response to increased eutrophic levels, environmental disturbance, or water stratification related to continental runoff or riverine input (Wade and Bown 2006; Ćorić and Hohenegger 2008; Auer et al. 2014). This latter assemblage correlates with warmer-water SST indicated by planktic foraminiferal assemblages (Supplementary Fig. 2). A second group of benthic foraminifera (Heterolepa dutemplei, P. ariminensis, Cibicides gr. refulgens, Discorbinella bertheloti, Cibicidoides gr. pachyderma, Biasterigerina planorbis, and Cibicidoides lobatulus) suggests oligotrophic, well-oxygenated background conditions under the influence of strong bottom currents (Donnici and Barbero 2002; Schönfeld 2002; Szarek et al. 2006; Fontanier et al. 2008; Koho et al. 2008; Schweizer et al. 2009; Frezza et al. 2010; Buosi et al. 2012). Background oligotrophism and high oxygen levels are suggested by the proportional dominance of the “sensitive” and “indifferent” EGs of benthic foraminifera (Fig. 7) and warm-oligotrophic planktic foraminifera (Fig. 8A). Note here that the sample from the Terebratula biostrome (TB), however, is mostly characterized by oxiphylic species such as C. gr. pachyderma and C. gr. refulgens, which are suspension-feeding epizoans (Koho et al. 2008; Schweizer et al. 2009). Our data from the Águilas Basin suggest, overall, that Terebratula preferred oligotrophic and well-oxygenated habitats, under moderate to strong currents, but tolerated mesotrophic conditions and fluctuating concentrations of oxygen levels at the seafloor (Fig. 11).

Modern (1981–2010) water temperatures off Águilas (SE Spain) range from 13.5°C to 28°C (range = 14.5°C). The monthly averages are between 17°C and 22°C, and the annual mean is 19.4°C (Guijarro et al. 2015) (Fig. 12A). The temperatures calculated from the oxygen isotope ratios (min = 17.1°C; max = 23.8°C; mean = 19.1°C; Table 1), as well as MAT (min = 16.3°C; max = 21.9°C; mean = 19.7°C), are consistent with the modern average temperatures. The MAT results are similar for the Terebratula-barren and Terebratula-bearing samples (Fig. 12B). This suggests that temperature alone does not explain the presence/absence distribution patterns of Terebratula during the Zanclean. Given that Terebratula lived in offshore environments, the temperatures derived from oxygen isotope ratios would not be representative of sea-surface conditions. The assumption is that, for the Zanclean, Mediterranean SST was higher than today (Templado 2014; Tindall and Haywood 2015). In the Águilas Basin, during the late Zanclean, other paleoclimatic proxies are the species Clypeaster cf. aegyptiacus, Echinolampas spp., Hyotissa sp., Talochlamys ercolaniana, Hinnites crispus, Spondylus crassicosta, and Gigantopecten latissimus (García-Ramos and Zuschin 2019). These species are considered to be warm-temperate to subtropical taxa (Brébion et al. 1978; Ávila et al. 2015). For these species, Monegatti and Raffi (2001) proposed monthly average SST of 24°C–25°C for at least five to six months. The relatively narrow range of temperatures obtained from oxygen isotopes (6.7°C; 3.6 °C without the 23.8 outlier; Fig. 12B) suggests that Terebratula either preferentially lived below the thermocline (cf. Houpert et al. [2015] for the seasonal thermoclines in the modern Mediterranean) or grew during particular seasons of the year. “Terebratula” sp. from Styria (Langhian, Austria), T. cf. calabra from Guadix (Tortonian, Spain), and T. scillae from Gallina (Calabrian, Italy) (Bojar et al. 2004; Clark et al. 2016; Rollion-Bard et al. 2016) also display relatively narrow temperature ranges: 4.9°C, 3.9°C, and 2.4°C, respectively. An exception is Pliothyrina maxima, another member of Terebratulinae, from the Coralline Crag (Zanclean, UK), showing a large range of 11.1°C (Vignols et al. 2018). Importantly, however, the data in those studies are derived from a sclerochronological approach on one or two specimens, sometimes mixing signals from the primary and secondary brachiopod shell layers. Oxygen isotope ratios measured on the primary shell layer and posterior and anterior parts of the shell should be interpreted with caution, because they are likely affected by nonequilibrium isotope fractionation (e.g., Romanin et al. 2018). The similar brachiopod species Liothyrella uva from Antarctic waters is subject to temperature ranges of only 2°C–3°C, not surviving above 4.5 °C (Peck 2005). In contrast, the cool-temperate Liothyrella neozelanica experiences values between 8°C and 18°C (Lee 1991). The relatively narrow temperature ranges of different Terebratula species should be further investigated as a possible cause for their extinction during the dramatic climate changes of the late Pleistocene, which would have left them insufficient time to adapt (Peck 2007). Note also that Terebratula not surviving beyond the Jaramillo Subchron coincides with the onset of the strongest glacial–interglacial shifts during the Pleistocene (e.g., Rohling et al. 2014).

Terebratula had an optimum close to the OTZ, but why was the genus rare or absent in the shoreface and basinward beyond the OTZ? Our data suggest that oxygenation levels, food availability, or temperature—which have been invoked to explain the distribution of some brachiopod species (e.g., Tunnicliffe and Wilson 1988; Kowalewski et al. 2002; Tomašových et al. 2006; Peck 2007)—might not have acted as limiting factors in the Águilas Basin. This is because Terebratula appeared to be fairly tolerant to local variations in these parameters.

The consistent pattern in the Águilas Basin is the peak abundance of Terebratula in sediments devoid of coralline algae. Accordingly, light penetration might have exerted an important influence on the distribution of Terebratula in relation to grazing pressure (Noble et al. 1976; Witman and Cooper 1983; Asgaard and Stenthoft 1984; Asgaard and Bromley 1991; Tomašových 2008; Zuschin and Mayrhofer 2009; Radley 2010). The bioerosion assemblages in the Terebratula biostrome are interesting because they display affinity with the Gnathichnus ichnofacies (Bromley and Asgaard 1993; De Gibert et al. 2007). Bromley and Asgaard (1993) proposed the Entobia ichnofacies for deep tier–dwelling borings in littoral rocky substrates subject to long exposure. The Gnathichnus ichnofacies, instead, characterizes shallow-tier structures on briefly exposed substrates (i.e., shells) in deeper water (De Gibert et al. 2007). The outcrop from the Pliocene Roussillon Basin (De Gibert et al. 2007), for example, was a shoreface environment whose shell beds are mostly composed of ostreids and pectinids. There, the dominant traces are G. pentax (rasping traces produced by regular echinoids feeding on algae) and Radulichnus inopinatus (produced by the radular grazing activity of gastropods or polyplacophorans) (De Gibert et al. 2007). In contrast, the dominant macrobioerosion trace in the Terebratula biostrome is Entobia isp., whereas other traces such as G. pentax (Fig. 4L) are rare (Molinu et al. 2013). García-Ramos and Zuschin (2019) argued that the rare occurrence of the bioerosion trace Gnathichnus on Terebratula shells, coupled with the absence of Radulichnus, may indicate dim light or aphotic conditions, in line with interpretations elsewhere (Bromley 2005). This interpretation is consistent with the microendolith assemblages reported by Molinu et al. (2013) from the Terebratula biostrome. That assemblage is dominated by fungal traces (Saccomorpha clava, Orthogonum lineare, Flagrichnus isp.), whose producers are more common in relatively deep aphotic environments (e.g., Glaub 2004; Wisshak 2012). The microendolith trace R. clavigera is also present in the Terebratula biostrome. The producer of this trace characterizes the euphotic zone (Golubić and Radtke 2008), but can also be common in deep euphotic habitats (Wisshak 2012). Overall, the microendolith assemblages suggest irradiances around 0.01% or less (Wisshak 2012), which classifies the habitat from the Terebratula biostrome as a transition from upper to lower circalittoral (Cameron and Askew 2011). Such microendolith ichnoassemblages pointing to dim light are consistent with the absence of coralline algae in the Terebratula samples. The dysphotic zone (1–0.01% irradiance) can occur at depths between 40 and 120 m in some Mediterranean localities (Ballesteros 2006). Recent shells of the brachiopod Gryphus vitreus are affected by the endolithic green alga Ostreobium queketti (which is adapted to extremely low light conditions) to a depth of 130–135 m off Corsica (Emig 2018). Similar benthic foraminiferal assemblages as in the Terebratula biostrome (Fig. 6) occur at 70–80 m depth in the Strait of Bonifaccio (Buosi et al. 2012), which supports the above interpretation based on microendolith ichnoassemblages.

The rarity or absence of Terebratula in shoreface environments was discussed earlier, but in the Águilas Basin the distribution of this brachiopod in upper offshore (circalittoral) environments is not homogeneous (Fig. 2A). The peak density occurs in the foreset–bottomset transition (i.e., the lower rollover of the clinoforms). Why was Terebratula rare or absent in the upper foreset? The rhodalgal facies disappears basinward beyond the OTZ (coincident with the topset–foreset transition; i.e., upper rollover zone) (Fig. 2A). This suggests that the upper foreset was already under poorly lit conditions or other excluding environmental factors for coralline algae were present. Terebratula could have potentially colonized this habitat, which was safe from high grazing pressure, but it was absent. This pattern can be explained by peak sedimentation rates at the upper foreset (upper offshore, circalittoral) due to its proximity to the upper rollover, which is a threshold between high and low hydrodynamic conditions (from sediment advection and bypass at the topset to sedimentation at the foreset) (Driscoll and Karner 1999; Cattaneo et al. 2007; Mitchell 2012). Thus, the environmental conditions at the upper foreset likely surpassed the physiological capacity of Terebratula to cope with sedimentation (Williams et al. 2018). Sedimentation rates decrease gradually down the foreset (Mitchell 2012). Accordingly, the facies belt at the foreset–bottomset transition marks the threshold of sedimentation rates tolerated by Terebratula.

Terebratula density rapidly decreases beyond the lower rollover basinward (Figs. 2A and 3C). This is coincident with Gaetani (1986) and Barrier et al. (1987), who report T. calabra as being restricted to proximal circalittoral environments, whereas other Terebratulidina such as T. scillae, G. vitreus, and Stenosarina sphenoidea display a deeper optimum in the lower circalittoral and bathyal (Gaetani and Saccà 1985). A common pattern of subaqueous delta-scale clinoform systems is the association with offshore currents running parallel to the shoreline along the foreset–bottomset transition (e.g., Pomar et al. 2002; Cattaneo et al. 2007; Patruno and Helland-Hansen 2018). The occurrence of aggrading sandwave fields migrating basinward at the localities of Terreros and La Carolina (García-Ramos and Zuschin 2019) demonstrates the sustained development of cyclonic current systems during the late Zanclean in the Águilas Basin (Fig. 1C). This interpretation is supported by the dominance in the benthic foraminiferal assemblages of suspension-feeding species that favor vigorous currents and by the abundance in the Terebratula biostrome of clionaid sponges, which cannot cope with high levels of turbidity (Carballo et al. 1994; Taylor et al. 2003). This pattern does not seem unique to the Águilas Basin, because records of Terebratula at the lower rollover and adjacent bottomset of delta-scale clinoforms are described elsewhere (e.g., Llompart and Calzada 1982; Pomar and Tropeano 2001; Soria et al. 2003; Videt 2003; Gramigna et al. 2012; Massari and D’Alessandro 2012; Reolid et al. 2012). As brachiopods are facultatively active suspension feeders (La Barbera 1977; Wildish and Kristmanson 1997), the attenuation or disappearance of the along-slope currents basinward possibly affected the Terebratula paleocommunity; this reflects the physiological cost of shifting from passive to permanently active suspension feeding (La Barbera 1977; James et al. 1992). An analogous case of offshore and onshore decrease in population density as a function of bottom current velocity has been described for extant communities of the terebratulidine G. vitreus (Emig 1989; Emig and García-Carrascosa 1991).

Cenozoic to Quaternary Terebratulinae have often been found in relatively shallow-water environments, from boulders at the toe of beach cliffs (Aigner 1983; Dixon 2011; Betancort et al. 2014) to intertidal and very shallow subtidal gravelly bottoms (Diedrich 2012). They often occur adjacent to sandwave fields where tidal or other types of currents are present but attenuated (Barrier et al. 1987; Roetzel et al. 1999; Pomar and Tropeano 2001; Courville and Crônier 2003; Kroh et al. 2003; Bosselaers et al. 2004; Calvo et al. 2012; Reolid et al. 2012). Such conditions prevent burial of the brachiopod paleocommunity by migrating subaqueous sandwaves. These brachiopods have also been found close to submarine hard-bottom structures that they possibly colonized, including shallow-water submarine cliffs (Kroh et al. 2003; Pervesler et al. 2011). Smaller species, such as Maltaia maltensis and “Terebratula” styriaca, were able to inhabit crevices and sheltered microenvironments in coralline algal buildups (e.g., Bianucci et al. 2011; M. Harzhauser personal communication 2019), cavities in coral reefs (Barbera et al. 1995), or in parareefal environments (Conesa et al. 2007). In the Pleistocene, some species colonized dykes in seamounts and walls on paleocliffs at bathyal depths (Ietto and Bernasconi 2005; Titschack et al. 2005). Other species have also been reported from muddy bottoms at bathyal depths (Thomsen et al. 2005; Rögl et al. 2008), but it remains to be determined whether these assemblages were transported. Most occurrences, however, are associated with environments close to the lower shoreface–offshore transition and upper offshore settings, often co-occurring with bryozoans and/or acorn barnacles (De Porta et al. 1979; D’Alessandro and Iannone 1982; Gaetani 1986; Studencki 1988; Taddei Ruggiero 1996; Bitner and Pisera 2000; Montenat et al. 2000; Gramigna et al. 2008; Pavia and Zunino 2008; Puga-Bernabéu et al. 2008; Di Stefano and Longhitano 2009; Messina et al. 2009; Long and Zalasiewicz 2011; Giannetti et al. 2018, 2019; Crippa et al. 2019). Terebratulid colonization of the photic zone and their co-occurrence with coralline algae have sometimes been explained by the onset of eutrophic conditions triggered by upwelling, which can be deleterious for phototrophic and mixotrophic organisms (Brandano et al. 2016). Other authors interpreted that Terebratula, which is often typical of monospecific to paucispecific assemblages, was an opportunist able to colonize environments subject to disturbance associated with mesotrophic conditions (Massari and D’Alessandro 2012). Overall, the available evidence points to Terebratula preferring habitats where grazing disturbance was reduced because of poorly lit environments (e.g., Pedley and Grasso 2002; Brandano et al. 2015).

The late Zanclean deposits of the Águilas Basin record cycles of Terebratula paleocommunities that developed offshore (circalittoral) on fine-grained sediments deposited at the foreset–bottomset transition of subaqueous delta-scale clinoforms. These deposits therefore provide a good scenario to understand the paleoenvironmental distribution of this taxon in space and time. The analysis of benthic and planktic foraminiferal and nannoplankton assemblages suggests that, overall, Terebratula thrived in relatively warm, oligotrophic to mesotrophic, well-oxygenated environments influenced by strong bottom currents. The oxygen isotopes showed that Terebratula in this basin lived in a relatively narrow range of temperatures (6.7°C). Such narrow ranges have also been reported for other species, potentially helping explain their extinction during the abrupt climate changes of the late Pleistocene: these brachiopods may have been unable to adapt quickly enough to such high-amplitude seawater temperature fluctuations after the Jaramillo Subchron. The consistent occurrence of terebratulids in sediments devoid of coralline red algae, combined with the bioerosion assemblages in the Terebratula biostrome, suggest that the limiting factor affecting the onshore distribution of Terebratula was light penetration and the associated high grazing pressure. This would explain the virtual absence of Terebratula in shoreface environments. Higher sedimentation rates at the shoreface–offshore transition also excluded the Terebratula paleocommunity in the upper foreset. In contrast, further offshore beyond the foreset–bottomset transition, we conclude that the attenuation or disappearance of along-slope currents was responsible for the lack of Terebratula populations.

We thank G. Romero Sánchez from the Paleontological Heritage Service of the Community of Murcia (Spain) for granting permits to conduct paleontological fieldwork in the study area. We also thank P. Bukenberger and J. Souto-Derungs (University of Vienna) for their help with SEM images and M. Bošnjak (Croatian Natural History Museum, Zagreb) for suggesting improvements to an earlier version of the article. E. Borghi, G. Díaz-Medina, M. Harzhauser (Natural History Museum of Vienna), and A. Giannetti (University of Alicante) are acknowledged for providing paleoenvironmental information on Terebratulinae from certain localities; E. Borghi and G. Crippa (University of Milan) also for sending literature. Special thanks to A. Tomašových (Slovak Academy of Sciences) for his help with technical aspects. D.A.G.-R. wishes to thank Y. Sun and E. Jarochowska (University of Erlangen–Nuremberg) for support in Erlangen. The authors thank M. Stachowitsch (University of Vienna) for proofreading and improving the language. We thank the editor J. Crampton for constructive comments during the review process. The insightful reviews of two anonymous reviewers substantially improved the final version of this article.