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

Studies of paleocommunities and trophic webs assume that multispecies assemblages consist of species that coexisted in the same habitat over the duration of time averaging. However, even species with similar durability can differ in age within a single fossil assemblage. Here, we tested whether skeletal remains of different phyla and trophic guilds, the most abundant infaunal bivalve shells and nektobenthic fish otoliths, differed in radiocarbon age in surficial sediments along a depth gradient from 10 to 40 m on the warm-temperate Israeli shelf, and we modeled their dynamics of taphonomic loss. We found that, in spite of the higher potential of fishes for out-of-habitat transport after death, differences in age structure within depths were smaller by almost an order of magnitude than differences between depths. Shell and otolith assemblages underwent depth-specific burial pathways independent of taxon identity, generating death assemblages with comparable time averaging, and supporting the assumption of temporal and spatial co-occurrence of mollusks and fishes. INTRODUCTION 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. MATERIAL AND METHODS Study Area and Target Species 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 Repository1; Fig. 1). This is an open shelf under the sedimentary input of 1GSA Data Repository item 2020173, details on methods and radiocarbon ages, is available online at http://www.geosociety.org/datarepository/2020/, or on request from editing@geosociety.org. Downloaded from https://pubs.geoscienceworld.org/gsa/geology/article-pdf/doi/10.1130/G47210.1/4974096/g47210.pdf by guest on 28 March 2020 2 www.gsapubs.org | Volume XX | Number XX | GEOLOGY | Geological Society of America 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 waveinduced 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 waveand 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]; GoodmanTchernov 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). Shell and Otolith Dating 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 graphitetarget 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. Age-Frequency Distributions 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 KolmogorovSmirnov (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 λ1 (initial phase of high disintegration in the TAZ not affected by burial) to a slow loss λ2 (a function of slower disintegration in the socalled 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 λ1 to λ2 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 oneand 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. RESULTS 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 Figure 1. Map of the Levantine Basin, easternmost Mediterranean Sea, with location of three collection sites off Ashqelon, southern Israel. Depth contours in m. 34.25° E SG10 SG30 SG40


INTRODUCTION
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(Kosnik et al., , 2013Krause et al., 2010;Tomašových et al., 2014Tomašových et al., , 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 , 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.

MATERIAL AND METHODS Study Area and Target Species
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 1 ; 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 waveinduced 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 dif-fer from bivalves with shell structure similar to Donax . 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).

Shell and Otolith Dating
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 graphitetarget 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.

Age-Frequency Distributions
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 λ 1 (initial phase of high disintegration in the TAZ not affected by burial) to a slow loss λ 2 (a function of slower disintegration in the socalled 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 λ 1 to λ 2 occurs at rate τ and can reflect the rate of diagenetic stabilization in the SZ (Tomašových et al., 2014. 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 oneand 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.

RESULTS
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 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 λ 1 were variable, they indicated that both shells and otoliths can   disintegrate rapidly after death, with yearly or decadal half-lives. The estimates of λ 2 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).

Effect of Burial on AFDs
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 λ 1 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 λ 2 (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 λ 2 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, λ 1 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 λ 1 , 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).

Consequences for Paleocommunity Studies
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., 2016Albano et al., , 2018Ritter 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.

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