Radiocarbon dating supports bivalve-fish age coupling along a bathymetric gradient in high-resolution paleoenvironmental studies

Paleoecological inferences about co-occurrence patterns and niche overlap assume that species present in the same sedimentary layer are of similar age and therefore potentially interacted with each other (Lyons et al., 2016). However, fossils preserved together within a single stratum can represent organisms that lived at vastly different times due to condensation, bioturbational mixing, and physical reworking (Kowalewski, 1996; Kidwell, 2013). Multiple examples of co-occurring shells of mollusks and brachiopods have been shown to differ significantly in median ages and time averaging (Kosnik et al., 2009, 2013; Krause et al., 2010; Tomašových et al., 2014, 2019). Such differences can be generated by intrinsic factors like between-species variation in skeletal durability (Kosnik et al., 2007; Kowalewski et al., 2018), in timing and duration of shell production (Tomašových et al., 2016), or in propensity to out-of-habitat transport. These intrinsic factors can be modulated or overwhelmed by gradients in extrinsic factors that influence burial and disintegration, such as sedimentation rates (Krause et al., 2010) and pore-water chemistry (Best et al., 2007).

Although paleoecological analyses are increasingly focused on whole ecosystems (Villéger et al., 2011; Roopnarine and Angielczyk, 2015), no studies have assessed time averaging of co-occurring species belonging to phyla with different ecosystem functions. Here, we quantified time averaging and modeled disintegration and burial of suspension-feeding bivalve shells and predatory nektobenthic fish otoliths along a 10–40 m depth gradient on the Mediterranean Israeli shelf to test the hypothesis that species co-occurring in the same death assemblage but subject to different intrinsic factors did not temporally co-occur in the original biological community. The shells and otoliths of our target species have comparable size and durability, but undergo different pathways after death. Infaunal bivalves are more likely to die and be buried in situ. In contrast, otoliths can be deposited far from the life location because they either originate from predated fish through feces (Nolf, 1995), implying that their final location depends on the predator range, or carcasses are made buoyant by bacterial decay gases and transported to the surface where they drift away (Elder and Smith, 1988), especially at temperatures greater than 16 °C (year-round in most temperate to tropical seas). Suspension-feeding bivalves and predatory fishes can further respond differently in terms of their population fluctuations to variation in environmental factors such as nutrient regimes due, e.g., to top-down controls of the trophic web. These differences can generate major variation in the structure of time averaging (determined by median ages and indicators of age range). In contrast to our expectation, we found that both taxa possessed very similar median ages and interquartile age ranges and that differences in age structure were smaller within depths than between depths. These results suggest that mollusks and fishes co-occurred temporally and spatially, and they point to the prevalence of depth-specific taphonomic and burial pathways independent of taxon identity.

We collected death assemblages with a Van Veen grab sampler at 10, 30, and 40 m depth off Ashqelon (southern Israel), Eastern Mediterranean, in autumn 2016 (Table DR1 in the GSA Data Repository¹; Fig. 1). This is an open shelf under the sedimentary input of the Nile River (Inman and Jenkins, 1984) and with seawater temperature between 17 °C and 30 °C. The fair weather wave base at 15–25 m (Hyams-Kaphzan et al., 2008) and strong wave-induced counterclockwise longshore currents that reach ∼35 m depth transport siliciclastic sands from the Nile Delta northwards (Golik, 1993; Avnaim-Katav et al., 2015) and limit the deposition of fine-grained sediments. This bathymetric gradient thus coincides with a decline in the exposure to wave- and current-driven erosion and reworking, leading to a decline in grain size from sands to muddy sands up to muds and to a higher net sedimentation rate (2.4 mm/yr at 40 m, exceeding the rates at 10 m [0.4 mm/yr] and 30 m [0.2 mm/yr]; Goodman-Tchernov et al. [2009] and our own unpublished data based on sediment cores at the two deepest sites). We dated the bivalve Donax semistriatus and the benthic fish Ariosoma balearicum (conger eel) at 10 m, and the bivalve Corbula gibba and a multispecies gobiid assemblage at 30 and 40 m (Table DR2). The targeted shells and otoliths are aragonitic (Degens et al., 1969). Although Corbula has conchiolin layers, which retard shell dissolution in waters undersaturated in calcium carbonate and increase shell strength, its taphonomic pathway does not differ from bivalves with shell structure similar to Donax (Gallmetzer et al., 2019). Corbula and Donax ranged in length from 2.5 to 5.6 mm (median 3.7 mm) and from 3.0 to 16.0 mm (median 5.0 mm), respectively. Otoliths ranged in length from 1.5 to 7.5 mm (median 3.5 mm). Their taphonomic pathways are poorly known, but otoliths are regarded as durable remains (Nolf, 1985).

Shells and otoliths were dated by accelerator mass spectrometry (AMS), using powdered carbonate targets (Bush et al., 2013). Four shells were also analyzed using the standard graphite-target method to assess the accuracy of carbonate targets, but we did not conduct such an assessment for otoliths due to their small size. Shell carbonate-target ages had a small offset for the youngest samples (Table DR6). However, we used the measured carbonate-target ages for both types of samples because a correction for the offset did not change our results (see the Data Repository and Table DR7). Radiocarbon ages were converted to calendar ages (see the Data Repository). We report all ages in calendar years before 2016 CE, the year of sample collection.

We computed measures of central tendency (median age), dispersion (interquartile range [IQR]), and skewness of the age-frequency distributions (AFDs) for each assemblage. The median ages were compared with the Wilcoxon test, and the shapes were compared with the Kolmogorov-Smirnov (K-S) test, obtaining p values with Monte Carlo simulations due to the small sample size (R package dgof; Arnold and Emerson, 2011). To determine the dynamics of shell loss from the sampled surface sediment layer (here assumed to represent the well-mixed taphonomically active zone [TAZ]) by disintegration and/or burial, we used three models: (1) a one-phase exponential model, defined by a single instantaneous per-individual loss rate λ, (2) a Weibull model with a gradual temporal decline in the loss rate, and (3) a two-phase exponential model that accounts for an abrupt temporal decline from a fast loss λ₁ (initial phase of high disintegration in the TAZ not affected by burial) to a slow loss λ₂ (a function of slower disintegration in the so-called sequestration zone [SZ] and the net rate of shell burial). The SZ can represent patches of sediment with less corrosive conditions within the TAZ or layers immediately below the TAZ; shells from these layers can be still exhumed back into the TAZ by burrowers or by storms, and thus be incorporated into AFDs as measured here. The decline in loss rate from λ₁ to λ₂ occurs at rate τ and can reflect the rate of diagenetic stabilization in the SZ (Tomašových et al., 2014, 2016). The otolith AFD at 10 m showed a drop in skeletal production in the past few centuries and a minimum age of 322 yr (Fig. 2), and the one- and two-phase exponential models were thus adjusted at this site with a parameter setting for the termination of production at 400 yr. We used the Akaike information criterion with correction for small sample sizes (AICc) to identify the best model. We computed the half-lives in the TAZ for the three models according to Tomašových et al. (2016) and the confidence intervals of AFD summary statistics and model parameters with a bootstrapping procedure with 10,000 iterations. We used Spearman rank correlations to assess relations between parameters and median ages and IQRs. All analyses were conducted in the R statistical environment.

The median ages of shell and otolith AFDs were similar within each depth, but they changed strongly and in parallel between depths, with median age equal to ∼600–700 yr at 10 m, increasing to ∼1400–1500 yr at 30 m, and declining to 13–200 yr at 40 m (Figs. 2 and 3; Table 1). At 10 m, both AFDs were right-skewed, and shape (K-S D = 0.5, p = 0.08) and median age (Wilcoxon W = 125, p = 0.62) did not differ between shells and otoliths. The mode of the otolith AFD was at ∼500 yr, followed by a sudden decline in the most recent centuries with a minimum age of 322 yr. This decline in production in the latest few centuries led to a smaller IQR in otoliths (156 yr) than in bivalves (962 yr). At 30 m, the shape (K-S D = 0.2, p = 0.98), the median ages (Wilcoxon W = 115.5, p = 0.66), and IQRs (1344 and 1499 yr) did not differ between shells and otoliths. At 40 m, the shape (K-S D = 0.5, p = 0.103) and the median age (Wilcoxon W = 33.5, p = 0.093) also did not differ between the two taxa, although differences in IQRs (292 and 18 yr for shells and otoliths, respectively) were slightly larger than at 30 m. To summarize, the differences in median age and IQR between depths exceeded several centuries up to 1000 yr, whereas within-depth differences were less than a few centuries (Fig. 3).

A one-phase exponential model best explained the dynamics of loss of four assemblages, and a two-phase exponential model best explained the other two (Table DR5). However, differences between one- and two-phase models in these four assemblages in AICc were less than six to eight units; i.e., model likelihoods were effectively equal, allowing us to compare assemblages on the basis of the two-phase model (Fig. DR1). Although estimates of λ₁ were variable, they indicated that both shells and otoliths can disintegrate rapidly after death, with yearly or decadal half-lives. The estimates of λ₂ showed that through time, both otoliths and shells shifted to half-lives equal to 500–1200 yr at 10 m, to ∼1200 yr at 30 m, and to 100–200 yr at 40 m (Fig. DR1).

The major role in shaping the differences in time averaging along the depth gradient was played by the gradients in shell burial rather than by gradients in initial disintegration for two reasons. First, the rate of initial disintegration λ₁ did not correlate with median ages or IQRs (Spearman r for median = −0.37, p = 0.49; r for IQR = 0.03, p = 1), although it was variable due to constraints from small sample size and temporally variable production. Additionally, the rate of loss λ₂ (which reflects burial and/or later-stage disintegration) correlated negatively with median ages (r = −0.94, p = 0.006) and IQRs (r = −0.88, p = 0.02). Second, faster burial at the deepest station was suggested by a sedimentation rate (2.4 mm/yr) that was one order of magnitude greater than at 10 m (0.4 mm/yr) and 30 m (0.2 mm/yr). Medians and IQRs of both shells and otoliths positively correlated with the inverse of sedimentation rates (Spearman r for median = 0.96, p < 0.002; r for IQR = 0.84, p = 0.037), and λ₂ correlated with sedimentation rate (r = 0.89, p = 0.02), indicating that this parameter indeed captured the time scale of shell burial below the TAZ, being greater at 10 and 30 m (millennial scale) than at 40 m (centennial scale). In contrast, λ₁ did not vary with sedimentation rate (r = 0.24, p = 0.64). At 40 m, initial loss was also facilitated by more aggressive pore waters; shells were more brittle than at shallower stations, and none was in pristine condition despite their young age. Accordingly, the two-phase exponential model showed the largest λ₁, also indicating that shells disintegrated at a higher rate at this site. The between-depth differences in sedimentation rates correlated with sediment distribution: Sands at 10 m pass to muddy sands and muds at 30 and 40 m (Table DR1) and thus can reflect stronger bypassing and winnowing of fine-grained sediments driven by wave-induced longshore currents at shallower depths (Stanley, 1989). Sea-level rise did not constrain the maximum age of our studied species along the transect because even the 10 m site was flooded ∼8000 yr ago, and the Holocene sea-level rise was limited to 2 m in the past 5000 yr (Sivan et al., 2004). This prevalence of extrinsic factors can be counteracted only by major differences in durability, leading to differences in median ages and IQR by multiple orders of magnitude (Kowalewski et al., 2018).

Quantification of the time resolution of marine fossil assemblages has been limited so far to primary consumers, such as foraminifera, mollusks, and brachiopods (e.g., Kidwell et al., 2005; Krause et al., 2010; Albano et al., 2016, 2018; Ritter et al., 2017). Fishes are important components of marine trophic webs, and their otoliths preserve species-specific morphology and great abundance in marine and lake sediments (Nolf, 1985) and offer paleobathymetric, paleoclimatic, and other paleoenvironmental information (Agiadi et al., 2018). Comparisons between shells and otoliths are meaningful because of their comparable durability due to similar size, mineralogy, and microstructure. However, otolith postmortem pathways can lead to significant out-of-habitat transport, while infaunal bivalves are more likely to remain buried in situ after death. Additionally, the changes in water and nutrient discharge of the Nile River over the entire Holocene (Hassan et al., 2012; Sun et al., 2019) can generate differences in temporal production of groups at different trophic levels. However, our findings indicate that these two intrinsic factors were ultimately negligible and that depth-specific burial pathways independent of taxon identity dominantly contributed to the formation of assemblages with comparable time averaging. Age distributions were not homogenized by cross-shelf transport processes on the high-energy shelf, and these organisms thus spatially and temporally co-occurred in the original living assemblage.

Sampling in Israel was conducted in the framework of the project “Historical ecology of Lessepsian migration” funded by the Austrian Science Fund (FWF) P28983-B29 (principal investigator: Albano). We thank Bella S. Galil for her support throughout the planning and running of the Lessepsian project, and Erin Dillon for discussions. Otolith dating was supported by research grant PA-RG201803 from the Palaeontological Association (UK) to Agiadi. Shell dating was supported by a grant by the University of Vienna (Austria) to Zuschin. Tomašových was supported by the Slovak Research and Development Agency (APVV17–0555). Jan Steger, Danae Thivaiou, Katherine Whitacre, and Jordon Bright prepared samples, which were analyzed for radiocarbon at the University of California−Irvine Keck Accelerator Mass Spectrometer Laboratory. Jan Päßler analyzed the sediment granulometry. Itay Katzman and the crew of the Mediterranean Explorer vessel helped throughout the field work. Three anonymous reviewers provided useful comments on the initial manuscript.