Temporal scales, sampling designs and age distributions in marine conservation palaeobiology Free

Supplementary material: R language scripts that replicate analyses are available at https://doi.org/10.6084/m9.figshare.c.6487005

Conservation palaeobiology uses the data on geobiological dynamics of past ecosystems to guide conservation decisions, to predict biotic responses to future environmental changes, and to contribute to the basic eco- and evolutionary foundation of conservation biology (Jackson et al. 2001; Dietl et al. 2015; Kidwell 2015; Rick et al. 2016; Barnosky et al. 2017; Tyler and Schneider 2018; Fordham et al. 2020). Fossil assemblages in outcrops and cores, together with surface accumulations of dead skeletal remains (i.e. death assemblages), provide a direct insight into rates and magnitudes of ecosystem responses to natural and anthropogenic environmental changes at time-scales that greatly exceed the duration of ecological surveys. High temporal extent of the stratigraphic records with fossil and death assemblages is especially important given the long history of human transformation of marine and terrestrial systems lasting for centuries to millennia in many parts of the world (Jackson et al. 2001; Lotze et al. 2006; Kidwell 2015; Stephens et al. 2019). In addition, deep-time records, particularly those from the Cenozoic and thus dominated by extant taxa, provide an opportunity to assess the resilience of marine ecosystems during past episodes of global warming, deoxygenation and acidification, representing past analogues of future environmental states.

The approaches used in conservation palaeobiology rely on integration of geohistorical records with ecological data on modern ecosystems. This integration occurs either (1) directly as in studies where assemblages observed alive are compared with death assemblages preserved in seabeds and in sediment cores; or (2) indirectly, when inferences from the deep-time fossil record are used to constrain models of ecosystem responses to future climatic or oceanographic states. Thus, one of the fundamental challenges in applying fossil data to conservation problems is the disparity in the temporal scale between palaeoecological and ecological records. Palaeolimnological and palynological records deposited in freshwater systems over the past c. 10 000 years tend to have a high stratigraphic resolution, and therefore serve as excellent ecosystem and climate archives (Smol 1992; Jackson et al. 2000; Last and Smol 2001; Brewer et al. 2012; Wolfe et al. 2013), as is the case with their deep-time analogues (Wang et al. 2019). Depositional processes in marine systems, especially in anoxic deep-sea and semi-enclosed basins, in river-dominated marginal marine environments, and in coral reefs, which are frequently characterized by persistently-high sedimentation or accretion rates, can also generate successions with high stratigraphic resolution, both over the past centuries (Osterman et al. 2005; Thibodeau et al. 2006; Rabalais et al. 2007; Yasuhara et al. 2007; Tsujimoto et al. 2008; Gooday et al. 2009; Roff et al. 2015; Wingard et al. 2017; Irizuki et al. 2018) and millennia (Baumgartner et al. 1992; Finney et al. 2000; Wilkinson et al. 2014; Bhattacharya et al. 2019; Hage et al. 2022; Palmer et al. 2022), as well as in older stratigraphic intervals (Miller and Eriksson 1997; Campbell and Nesbitt 2000). However, the marine stratigraphic record on open continental shelves is typically characterized by temporally-variable sedimentation rates, with phases of omission or very slow sedimentation, and tends to be more bioturbated than the record of lakes, deltas or silled marine basins (Barrell 1912, 1917; Schindel 1980, 1982; Anders et al. 1987; Sadler and Strauss 1990; Bianchi and Allison 2009). This incompleteness and low resolution of sedimentary successions deposited in marine systems can be partly compensated by their temporal persistence over hundreds of thousands to millions of years, providing windows into long-term ecosystem responses to perturbations analogous to the present-day or future environmental changes (Yasuhara et al. 2008, 2017, 2020a; Zuschin et al. 2014; Caswell and Frid 2017; Burke et al. 2018; Dominici et al. 2018; Frieling and Sluijs 2018; Keller et al. 2018; Lewandowska et al. 2020; Spalding and Hull 2021).

In marine systems, the temporal scale of death assemblages (here, accumulations of dead skeletal remains located in the mixed layer in the uppermost part of the sedimentary column affected by bioturbation and physical mixing) and fossil assemblages (here, subsurface accumulations corresponding to Holocene or older historical layers in cores or outcrops) in the stratigraphic record typically varies by several orders of magnitude (Flessa and Kowalewski 1994; Kidwell 1997, 2013; Olszewski 1999). The distinction between death and fossil assemblages is not based on their age but rather on whether they did or did not transit into the historical layers, below the depth of biotic or physical mixing and beyond the active input of newly-dead shells from living communities. This distinction also underlies conceptual differences between various sampling designs used in conservation palaeobiology. The temporal scale of death and fossil assemblages can be reduced to two main components: (1) stratigraphic resolution, i.e. temporal separation among increments (sampled sedimentary layers); and (2) time averaging of fossil assemblages (their temporal grain), i.e. the age range of skeletal remains within increments, generated by mixing of remains of organisms that died at different times (also called microstratigraphic acuity, Schindel 1982; Behrensmeyer et al. 2000; or depositional resolution, Kowalewski 1996; Kowalewski and Bambach 2008). Both components vary significantly through time and among marine depositional environments as documented (1) by the negative dependence of sedimentation rates on the timespan over which a given sediment thickness was deposited (Sadler 1981) and (2) by time averaging of centimetre-scale stratigraphic increments varying from decadal to multi-millennial or even longer scales (Scarponi et al. 2013; Tomašových et al. 2018, 2022a; Zimmt et al. 2022). These empirical disparities in scale are also predicted by sequence-stratigraphic models (Kidwell 1986, 1991; Holland 2000) and by taphonomic and bioturbation models (Olszewski 2004; Tomašových et al. 2014). The variability in temporal scale is further exacerbated in analyses of assemblages derived from various types of records, e.g. when comparing palaeoecological with ecological (Kidwell 2001, 2013; Hohenegger et al. 2018; Jonkers et al. 2022) or archaeological data (Newsome et al. 2007; Rick and Lockwood 2013).

This variability in temporal scale directly determines whether the transformation of the stratigraphic record of ecosystem changes to their original chronological sequence is possible or not and what sampling design is the most effective in capturing the magnitude and rate of ecological turnover before and after the anthropogenic impacts. For example, as high time averaging coarsens temporal resolution, the proportion of dead specimens sourced by a new post-impact community state is expected to be small relative to the proportion of pre-impact specimens in the same time-averaged death assemblage accumulating in the surface mixed-layer (Kidwell 2008). Therefore, time averaging triggers high inertia of death assemblages to recentmost regime shifts as observed in the composition of the living community. Uppermost core increments formed by death assemblages located in the mixed layer thus do not necessarily fully capture post-impact community states. Low stratigraphic resolution and high time averaging can further obfuscate the timing and magnitude of temporal turnover in ecosystem composition (Scarponi and Kowalewski 2007; Tomašových and Kidwell 2010; Hunt 2012; Nawrot et al. 2018; Zimmt et al. 2021), systematically modify the apparent onset or offset of events (Guinasso and Schink 1975), and generate ecologically misleading co-occurrence patterns (Peng and Broecker 1984; Paull et al. 1991; Sepulcre et al. 2017; Ausín et al. 2019) or food-web reconstructions (Roopnarine and Dineen 2018; Shaw et al. 2021). Owing to the changes in temporal scale within cores or in comparisons between non-averaged living and time-averaged death or fossil assemblages, the nature of the marine stratigraphic record can thus generate apparent compositional or diversity shifts in fossil assemblages in the absence of true changes in the original past living communities (Tomašových and Kidwell 2009). If the variation in temporal scale is not accounted for in palaeobiological analyses, such cross-scale analyses can lead to ill-informed palaeoecological inferences (Miller 1986; Kowalewski et al. 1998; Holland and Patzkowsky 2002) and consequently incorrect conservation or management decisions.

Here, we outline some of the commonly used sampling designs in marine conservation palaeobiology, as exemplified by the studies included in this Special Publication volume (Nawrot et al. 2023), and place them in a common conceptual framework. We discuss how changes in temporal scale imposed by different study designs affect the interpretation of palaeoecological data and suggest research strategies that can mitigate this problem.

Palaeobiologists use various records to measure the magnitude of turnover between the baseline (typically Holocene or late Pleistocene) and novel (present-day or post-impact) ecosystem states (Williams and Jackson 2007; Willis et al. 2010; Yasuhara et al. 2012a; Kosnik and Kowalewski 2016). These records can be used to estimate the historical range of variability in ecosystem composition, assess the magnitude and rate of change induced by anthropogenic impacts, or set targets for restoration and evaluate the degree of ecosystem recovery. The history of past ecosystems as observed in the pre-Pleistocene stratigraphic record – often referred to as deep-time record – further represents a natural eco-evolutionary laboratory that allows better understanding of long-term responses to past environmental changes, which in turn enable predictions for the future (Saupe et al. 2014; Sibert et al. 2016; Pimiento et al. 2017; Salvatteci et al. 2022; Yasuhara and Deutsch 2022; Yasuhara et al. 2022). These two complementary research goals – setting the context for present-day conditions and using past environmental perturbations as potential analogues of future environmental change – are often listed as specific to the near-time and deep-time approach, respectively (Dietl and Flessa 2011; Dietl et al. 2015).

Three types of conceptually different, but not mutually exclusive sampling designs were used in analyses of anthropogenic impacts and in the detection of baseline and novel ecosystem states on the basis of palaeoecological and ecological data: (1) live–dead, (2) live–fossil and (3) dead–fossil analyses (Fig. 1). They differ in (1) whether the surface death assemblages, actively sourced by input of shells from living assemblages and collected by grabs, boxcores or multicorers, represent the youngest or the oldest ecosystem state, and (2) whether temporal scales of assemblages compared in these pairwise analyses are similar or not. Death assemblages that were formed in a marine environment but are now exposed on emergent reef flats or on deltaic plains (due to the infill of accommodation space, uplift or eustatic sea-level fall) and thus lack active input of recently dead shells, conceptually belong to the fossil assemblages that are otherwise typically located in historical layers of sediment cores (Fig. 1b). Additional or more complex sampling designs can arise, for example, when living assemblage surveys are also available from the late nineteenth or early twentieth century (e.g. Zu Ermgassen et al. 2012; Thurstan et al. 2013), potentially predating the timing of major anthropogenic impacts and the deposition of the late twentieth century sediments sampled in cores or in surface death-assemblages.

This type of sampling design includes comparisons between surface death assemblages and living assemblages. The pre-impact states are expected to be at least partly captured by grab-collected or top-core death assemblages located in the mixed layer, time-averaged to several centuries or millennia during the Late Holocene (e.g. Pandolfi and Minchin 1996; Greenstein and Pandolfi 1997; Zuschin et al. 2000; Kidwell 2007; Zuschin and Ebner 2015; García-Ramos et al. 2016; Martinelli et al. 2016; Clark et al. 2017; Dietl and Smith 2017; Rillo et al. 2019). These types of death assemblages are typically contrasted with living assemblages sampled during ecological surveys conducted in the twentieth or twenty-first centuries at the same or nearby locations. Beach death assemblages or the records of cetacean strandings, sourced by landward transport from intertidal, subtidal or pelagic habitats, can be contrasted with contemporary regional-scale living assemblages compiled on the basis of surveys in intertidal and subtidal habitats (Warwick and Light 2002; Pyenson 2010; Archuby et al. 2015; Rojas and Martínez 2020; Bhattacherjee et al. 2021).

On the one hand, strong shifts in the composition of assemblages of organisms with durable skeletal remains, such as foraminifers, corals, or molluscs, that occurred since the nineteenth and twentieth centuries, when coupled with slow sedimentation and high sediment mixing will typically generate death assemblages strongly dominated by pre-impact cohorts (so-called taphonomic inertia; Kidwell 2008). Top-core increments are thus frequently inert to recentmost changes owing to their time averaging. In other words, older pre-impact cohorts are so abundant that they swamp any contribution of new post-impact cohorts in terms of their proportional abundance. Impact-driven changes in the composition of living assemblages, increasing the relative contribution of post-impact cohorts in the surface death assemblages, thus propagate into the death assemblage composition slowly (Cameron 1995; Kidwell 2008) or with delay, like when a death assemblage accumulates by transportation of remains from adjacent habitats (Hammerman et al. 2022). Impacted benthic ecosystems tend to be characterized by higher live–dead mismatch than more pristine ecosystems because death assemblages are not immediately diluted by inputs from new, post-impact community states (Kidwell 2007, 2008; Leshno et al. 2015; Gilad et al. 2018; Michelson et al. 2018). Similarly, assemblages of planktic foraminifers collected by moored sediment traps differ from top-core death assemblages as a function of the amount of temperature change at the given location (Jonkers et al. 2019), indicating that surface death assemblages are sufficiently inert to recentmost input and represent pre-industrial baseline states. On the other hand, strong shifts in the composition of communities will naturally propagate more quickly to death assemblages when coupled with high sedimentation, low mixing and/or fast disintegration of skeletal remains. However, the input of the post-impact cohorts to surface death assemblages inevitably counteracts the inertia at least to some degree even when sedimentation rates are very slow because younger cohorts are exposed to disintegration for shorter time than pre-impact cohorts, generating right-skewed postmortem age distributions (Olszewski 1999; Kidwell 2002). Some temporal autocorrelation between death and living assemblages can be expected when individuals in living assemblages represent surviving members of recentmost recruit cohorts that were partly incorporated into a death assemblage or are descendants of individuals that form a death assemblage (Tomašových and Kidwell 2011; Feser and Miller 2014). Observations of fine-scale spatial structure or patchiness preserved in time-averaged surface death assemblages can be partly expected from such contribution of recentmost cohorts although they can also indicate temporal persistence in the dynamic of benthic communities (e.g. Miller 1988; Arkle and Miller 2018; Casebolt and Kowalewski 2018; Hyman et al. 2019). Surface death assemblages thus generally have a capacity to simultaneously capture signatures of both older assemblage states, shaped by pre-impact cohorts if the timing of ecosystem shift occurred recently, and the youngest assemblage states. Among the studies included in this volume, Kokesh and Stemann (2022)  and Meadows et al. (2023)  provide examples of the live–dead approach and are discussed in more detail below.

These analyses typically include comparisons between subsurface fossil assemblages preserved in historical layers and living assemblages. The pre-impact ecosystem states can be captured in subsurface core increments below the mixed layer (Fig. 1a) but also in uplifted Pleistocene and Holocene outcrops (Fig. 1b) that are not affected by any input of cohorts from post-impact states (Pandolfi and Jackson 2006; O'Dea et al. 2014, 2020; Dillon et al. 2021; Rivadeneira and Nielsen 2022). Analyses that are based on sediment cores collected on land or on sediment cores collected in marine settings, focusing on the Middle Holocene core segments (Toth et al. 2019) or avoiding core tops affected by bioturbation (e.g. Barbieri et al. 2020), also fall into this category. The Holocene and older fossil assemblages can be contrasted with living assemblages sampled in the twentieth or twenty-first centuries from similar habitats in the same region. In contrast to the natural standardization of spatial resolution of assemblages obtained through grab-based sampling for live–dead analyses, this approach requires application of consistent methods that allow standardized estimation of diversity and composition in both living marine communities and fossil assemblages. In this volume, Ivkić et al. (2023)  assess the efficiency of two standard plotless transect methods (line intercept transects and point intercept transects with different interval length) and the most widely used plot method (photoquadrats) on diversity and community composition of late Pleistocene coral reefs of the Red Sea. They argue that plotless methods perform better than photoquadrats because the effect of time averaging is minimized. Point intercept transects with 20 cm intervals are most efficient when quantitatively sampling fossil coral reefs and the results can be directly compared to those gained with plotless methods in modern reefs without introducing a bias caused by sampling design.

The third type of sampling design includes comparisons between surface death assemblages and fossil assemblages from the subsurface layers. These analyses are typically not limited to pairwise comparisons but extend to analyses of all increments in the mixed layer and below it. In contrast to the live–dead analyses, dead–fossil analyses may not reach the youngest ecosystem state and thus may not fully capture the post-impact state. Many palaeobiological studies focused on within-core changes in composition of fossil and death assemblages representing the late-Holocene highstand sea-level phases (Aronson and Precht 1997; Edgar and Samson 2004; Armenteros et al. 2016; Cramer et al. 2017, 2020; Handley et al. 2020; Hong et al. 2021). One advantage of this approach is that the significant decline in time averaging between death and living assemblages is avoided in these comparisons (Kowalewski et al. 2015), although systematic increase in time averaging is still expected to occur between mixed-layer and historical-layer assemblages (Tomašových et al. 2023). We note that differences in the composition detected in live–dead and live–fossil comparison can be unrelated to ecological changes as they can be generated by intraspecific differences in durability (Cummins et al. 1986; Pandolfi and Greenstein 1997; Ford and Kench 2012) or lifespan (Edinger et al. 2001; Kidwell and Rothfus 2010; Cronin et al. 2018). The live–dead and live–fossil comparisons thus need to take these taphonomic effects into consideration; dead–fossil analyses can be expected to be less sensitive to these differences.

Bauder et al. (2022)  used this type of design and documented temporal constancy in the composition of foraminiferal assemblages between fossil and surface death assemblages in a sediment core and grab samples from the One Tree Reef (OTR) lagoon on the Great Barrier Reef over the past c. 400 years. However, they found that the composition of surface death assemblages from other lagoons (more impacted by human activities) differs from the OTR assemblages, indicating that the OTR lagoon remains relatively pristine and still reflects pre-industrial baseline states. Scarponi et al. (2023)  also used a conceptually similar approach to compare Mid-Late Holocene fossil assemblages sampled from the subsurface marine stratigraphic record of the Po Plain with surface death assemblages collected in shallow-subtidal environments of the NW Adriatic Sea. They found a significant compositional shift towards novel states recorded by the nearshore death assemblages, characterized by lower species richness, higher beta diversity, and increased abundance of the opportunistic bivalve Corbula gibba, scavengers and deposit feeders, most likely reflecting human transformation of coastal habitats and introduction of non-native species.

Rather than simply exhibiting two states before and after the impact, the real-world ecosystem dynamics tend to be more complicated, with multi-state shifts, successional stages, collapses followed by recoveries, or by transient states (Finney et al. 2010; Yasuhara et al. 2012a, b; Hughes et al. 2013; Seddon et al. 2014). For example, ecosystems may recover and return to baseline conditions after transient novel states due to active remediation or restoration efforts. Molluscan communities dominated by chemosymbiotic lucinids colonized organic-rich muds on the Palos Verdes Shelf (southern California) in the 1970s–90s after the onset of wastewater treatment enforced by new regulations; these lucinid communities were replaced eventually by more functionally diverse molluscs with Nuculana, partly returning to the state prior to the wastewater contamination (Leonard-Pingel et al. 2019; Tomašových et al. 2019a). Some subset of the surface death assemblages in this system is thus dominated by post-impact states, and consequently the contrast between living and death assemblages measures the strength of recovery rather than the earlier pollution-driven loss in functional diversity, especially at sites with faster sedimentation. Although the extraction of such complex history on the basis of surface death assemblages alone is difficult, their composition can still be informative when analyses are integrated with other sources of data such as geochronological age dating or long-term monitoring of living assemblages (Bizjack et al. 2017; Albano et al. 2018). For example, Kokesh and Stemann (2022)  compared surface death assemblages with living assemblages to assess the importance of short-term invasion of the Kingston Harbour (Jamaica) by the Asian green mussel Perna viridis. This mussel invaded the Caribbean Sea in the 1990s, expanded its geographical range, but declined in abundance since the 2010s. Their analyses indicate that surface death assemblages faithfully record the spatial variation in past invasion intensity and that relative abundances of other species in the living molluscan community have not returned to the pre-invasion state.

Ecological assessments of environmental pollution or sediment toxicity reveal fine-scale sensitivity of foraminifers, ostracods or molluscs to sewage, oil spills, organic enrichment or oxygen depletion (Bandy et al. 1964; Martin 2000; Schönfeld et al. 2012; Alves Martins et al. 2019). When live-collected data are scarce due to low standing population densities or lack of long-term monitoring data, death assemblages that integrate ecological history over longer timespans can be used to constrain sensitivity of species to environmental stressors (Kidwell 2013). For example, in this volume, Mamo et al. (2023)  document sensitivity of combined live and death assemblages of benthic foraminifers to eutrophication and oxygen depletion, exceeding the tolerance of corals that cannot tolerate highly eutrophic conditions. Benthic foraminifers thus can track the whole spectrum of trophic conditions. This approach, in which pristine states are based on the least contaminated segments of gradients, naturally extends to detection of anthropogenic impacts on the basis of microfossil assemblages in sediment cores (Alve 1995; Murray 2006; Scott et al. 2007; Dolven et al. 2013).

The use of benthic indices based on various community attributes (Benthic response index, Bentix, AMBI, M-AMBI, foram-AMBI) in the assessment of ecological quality of marine ecosystems on continental shelves reflects this exceptional sensitivity of benthic groups to anthropogenic impacts (Borja et al. 2000; Smith et al. 2001; Simboura and Zenetos 2002; Alve et al. 2016; Jorissen et al. 2018), and these indices are increasingly used in live–dead studies (Dietl et al. 2016; Smith et al. 2020). This approach assumes that a spatial gradient in community composition observed today accurately captures the baseline state. However, this implicit assumption is put into question by the findings that most present-day benthic ecosystems in which these indices were defined are not pristine and instead frequently reflect shifting baselines (Jackson 1997; Lotze et al. 2006). Therefore, this space-for-time substitution approach, in which living assemblages collected across space substitute for the lack of their temporal coverage, can miss the compositional state typical of non-impacted conditions. However, such pre-impact states are still captured in the stratigraphic record. For example, Hesterberg et al. (2020) documented that prehistoric oysters from archaeological middens were larger than those observed in oyster reef at locations in the Gulf of Mexico that were considered to be pristine. As a consequence, even in the absence of challenges associated with the detection and separation of baseline and novel states, longer temporal archives represented by the geohistorical data are needed to validate ecological assessment studies and to identify the most sensitive species that can be rare in present-day polluted ecosystems, but were important members of past communities (e.g. O'Brien et al. 2021). For example, Leshno et al. (2016) found that the spatial variation between polluted and control stations in live–dead mismatch positively correlates with the variation in AMBI or Bentix indices, meaning that high mismatch is associated with a shift from better status in death assemblages to worse ecological status in living assemblages. However, palaeoecological archives also show that low-diversity assemblages inhabiting naturally stressed environments can generate indices that suggest low ecological quality status even in the absence of anthropogenic impacts (Barbieri et al. 2020). In this volume, Smith et al. (2022)  found that time-averaged death assemblages can track baseline benthic community structure in terms of its ecological quality status as measured by the multivariate-AMBI index (Muxika et al. 2007). Living assemblages were exposed to environmental change and transformed to death assemblages in simulations of local communities and metacommunities. Reference conditions based on time-averaged death assemblages produced the same remediation decision that would have been made using a perfect knowledge of living assemblages for the entire duration of the simulation with 500 generations. This efficiency occurs because time averaging induces inertia and death assemblages thus retain the memory of pre-impact states. In contrast, remediation decisions were often incorrect when the reference conditions used in the multivariate-AMBI index were based on non-averaged living assemblages from the last ten generations of the simulation, illustrating the susceptibility of these assessments to shifting baseline syndrome.

The assessment of ecosystem state and range of variability in ecosystem composition naturally depends on the temporal duration of a study. The turnover exhibited by pre-Holocene fossil assemblages exposed to natural disturbances thus can be compared to turnover observed in present-day communities affected by anthropogenic impacts at multiple time-scales. Among the different drivers that have been found responsible for significant ecosystem change in the geological past, three covarying processes are of particular interest for conservation purposes: deoxygenation, acidification and abrupt warming. In this regard, pre-Holocene intervals allow assessments of not only extirpation or extinction phases, but also of recovery over longer time-scales. The warming events that were coupled with oxygen depletion and altered seawater chemistry during the Early Jurassic oceanic anoxic event (Caswell and Coe 2013; Caswell and Frid 2013; Danise et al. 2021; Atkinson et al. 2023), during the Paleocene–Eocene thermal maximum (Tian et al. 2021; Hupp et al. 2022) or during the Mid-Miocene Climatic Optimum (Steinthorsdottir et al. 2021), provide deep-time analogues of ecosystem responses to carbon cycle perturbations (Foster et al. 2018). For example, in this volume, Caswell and Herringshaw (2023)  found that the recovery of benthic communities and bioturbation is significantly prolonged in the wake of the Toarcian oceanic anoxic event in the Yorkshire Basin, in contrast to faster recovery on open shelves and in oceanic environments. The stratigraphic record of benthic communities based on trace fossils v. skeletal remains is decoupled, with faster recovery in the depth of mixed layer and diversity of trace fossils as opposed to abundance or diversity of body fossils.

Deep-ocean sediment cores encompassing hundreds of thousands of years of sedimentation have facilitated the palaeoceanographic and micropalaeontological studies of the Cenozoic record. They document that changes in temperature and oxygen controlled marine ecosystems on million-year, millennial and centennial time-scales as observed during the mid-Pliocene Warm Period and during the late Pliocene and Pleistocene high-frequency climate cycles (Bassetti et al. 2010; Huang et al. 2018; Yasuhara et al. 2020b). This interval includes the last major glacial–interglacial cycle, with ice sheets reaching their maximum positions between 26.5 and 19 ka, followed by a deglaciation and an abrupt rise in sea level (Clark et al. 2009). This phase represents the last among eight major glacial and interglacial cycles of about 100 kyr duration over the past 730 kyr (middle and late Pleistocene), triggering as many cycles of absolute sea-level variation of 130 m (Spratt and Lisiecki 2016). Late Pliocene and early Pleistocene climatic cycles were less severe, but nevertheless important, with periodicities of 41 and 23 kyr (Lisiecki and Raymo 2005) and cycles of eustatic variation of about 50 m (Miller et al. 2005, 2012; Rohling et al. 2014; Dumitru et al. 2021).

Considering that most present-day marine invertebrate species have stratigraphic ranges spanning millions of years, and thus persisted in the face of this eustatic and climatic variability, stratigraphic palaeobiology of Pliocene and Pleistocene high-frequency depositional sequences is a useful means to understand long-term ecosystem dynamics and to test consequences of climatic changes on marine ecosystems. On the one hand, the comparison of fossil assemblages periodically re-occuring in the same part of depositional sequences (Pandolfi 1996; Kowalewski et al. 2015), and those sampled across glacial–interglacial or stadial–interstadial cycles (Cannariato et al. 1999; Kitamura et al. 2000; Tager et al. 2010; Scarponi et al. 2022), generally indicate resilience of benthic communities to natural environmental changes. On the other hand, high-frequency and intense climate changes can coincide with significant range shifts that can lead to assembly of communities that do not have any analogues today or to extinctions of species with narrow thermal niches (Roy et al. 1995, 1996; García-Ramos et al. 2020; Orzechowski and Finnegan 2021). Species expand poleward during interglacials (Valentine and Jablonski 1993; Rasmussen et al. 2003; Aguirre et al. 2005; Monegatti and Raffi 2007; Garilli 2011; Yasuhara and Danovaro 2014; Albano et al. 2022) and retreat equatorward during glacials (Raffi 1986; Dominici 2001; Garilli 2011; Crippa et al. 2019; Borghi and Garilli 2022). Threshold climatic effects entail changes in the dynamics of ocean circulation and bottom-water oxygen fluctuations, affecting the deep sea biota of enclosed epicontinental basins, such as the Sea of Japan during the Middle Pleistocene (Mid-Brunhes Event, MIS 10–11, Gallagher et al. 2015; Huang et al. 2018) or the eastern Mediterranean during the last interglacial (Sapropel S5, Rohling et al. 2006; Marino et al. 2007; Capozzi and Negri 2009; Rodríguez-Sanz et al. 2017). In a scenario that predicts fast temperature increase for most oceanic regions in the coming decades (Mitchell et al. 2017), deep-time studies foresee its long-term, structural effect on the tropical ecosystems (Kiessling et al. 2012) and on the native fauna in temperate regions, further modified by recent biotic invasions from tropical seas (Albano et al. 2021).

Climatic models have focused on the study of the mid-Piacenzian Warm Period (mPWP) as the most recent past interval of sustained global warming (Haywood et al. 2013), with mean global temperatures higher by c. 2–3°C (Vega et al. 2020) and a global sea level higher by 16.2–17.4 m than today (Dumitru et al. 2019). The deep-time study of zooplankton diversity patterns in the central and North Atlantic Ocean has shown that diversity–temperature relationships have been remarkably constant through the mPWP and during the subsequent climatic and eustatic cycles, suggesting that species diversity is rapidly reorganized as the ranges of species respond to temperature changes on ecological time-scales (Yasuhara et al. 2012b). Similarly, Dominici and Danise (2022)  found that the effects of the mid-Pliocene warming did not change the Mediterranean onshore–offshore gradient in the composition of molluscan communities in the long term. However, the effect of warming disproportionately affected deeper outer-shelf and bathyal environments, consistent with the results of other studies on the response of the deep-sea benthos to major climatically driven oceanographic changes (Yasuhara et al. 2008, 2014). Danise and Dominici (2022)  found that the temporal change in diversity of bivalves was rather gradual as some species survived in temporary refugia in the warmer, eastern Mediterranean. More importantly, they found that species extinction in the last 5.3 Ma depended on occurrence frequency, geographical range and habitat specialization. Epifaunal, mobile pectinids, including species with high abundance and broad geographical ranges, exhibited higher extinction than infaunal, siphonate venerids and lucinids owing to their higher specialization to carbonate habitats undergoing loss and fragmentation, and to their high metabolic rates in relation to larger body size. The selectivity of climate-driven extinctions and extirpations in the geological past further inform estimates of modern extinction risk for poorly studied marine species and provide insights into macroevolutionary processes (e.g. Harnik et al. 2012; Finnegan et al. 2015; Saupe et al. 2015; Collins et al. 2018; Reddin et al. 2022).

One of the major challenges in integrating ecological data from living communities with death and fossil assemblages, and in transforming patterns and processes documented in the deep-time fossil record to present-day conservation problems is the disparity in the temporal scale. The importance of the scaling problem is magnified by three factors. First, stratigraphic resolution and time averaging can be decoupled and their empirical estimation is challenging (Kowalewski and Bambach 2008). Second, both aspects of the temporal scale can change systematically downcore (Tomašových et al. 2023). Finally, stratigraphic resolution and time averaging can vary systematically within depositional sequences in response to tectonic subsidence and global eustatic changes (Holland 2000; Scarponi et al. 2013). Below we focus on the first two factors and outline their implications for near-time conservation palaeobiology. More detailed discussion of the importance of the sequence stratigraphic architecture and facies shifts in determining palaeoecological patterns in more ancient successions can be found elsewhere (e.g. Kidwell 1991; Holland 2000; Hannisdal 2007; Patzkowsky and Holland 2012; Holland and Patzkowsky 2015).

Both stratigraphic resolution and time averaging are determined by geological and ecological processes in the mixed layer, i.e. by sedimentation, bioturbation and disintegration (Kidwell and Bosence 1991; Olszewski 2004; Tomašových et al. 2014, 2023). Analyses of ecological data are limited by the total temporal durations of time-series assembled over the past years or decades (Willis and Birks 2006; Wolkovich et al. 2014). The duration, or temporal grain of sampling events and temporal separation between them in monitoring time-series of living assemblages is typically fully controlled in ecological surveys. In contrast, these components of temporal scale cannot be manipulated when assessing fossil or death assemblages. The estimation of temporal scale is challenging even for sedimentary successions for which age models are available, and the mixing-induced age overlaps between adjacent stratigraphic increments do not have any equivalent in ecological time-series. First, interpolating the net sedimentation rate derived from age models to undated time intervals is biased by scale dependency of sedimentation rates (Fig. 2, Sadler 1981). Second, estimation of stratigraphic resolution and time averaging on the basis of age models requires not only that the sedimentation rates are accurately measured but also that the mixing depth is estimated (Soetaert et al. 1996; Walbran 1996; Boudreau 1998; Teal et al. 2010; Solan et al. 2019; Tomašových et al. 2019a; Díaz-Asencio et al. 2020; Edelman-Furstenberg et al. 2020; Song et al. 2022). Below, we show that simple age models can provide some approximation of stratigraphic resolution and time averaging. However, an accurate estimation of both components of temporal scale in the Holocene–Anthropocene record ultimately requires that age models are coupled with post-mortem age-frequency distributions of skeletal particles dated by amino-acid racemization, ¹⁴C or U–Th.

The relationship between stratigraphic resolution and time averaging is not random (Fig. 2) because time averaging will tend to be low if stratigraphic resolution (and sedimentation rate) is very high, regardless of mixing depths and the rate of disintegration in the taphonomically active zone (TAZ; Davies et al. 1989a; Powell 1992; Olszewski 2004; Berkeley et al. 2014; Tomašových et al. 2014, 2019a; Ritter et al. 2019). If sedimentation rate is extremely low and skeletal disintegration in the TAZ is very slow, low stratigraphic resolution will be coupled with high degree of time averaging. However, if sedimentation rate is extremely low these two components can also be fully decoupled, for example, when two rapidly deposited layers entombing weakly averaged assemblages are separated by a longer hiatus (at the scale of millions of years), because no skeletal remains deposited during the hiatus phase are preserved (Tomašových et al. 2020a).

Stratigraphic gaps can be neglected when they are shorter than the time-scale of the ecological process of interest. When the mixing depth can be estimated, time averaging can be approximated from the estimate of net sedimentation rate as a transit time needed to advect skeletal remains below a given increment or below a mixed layer (Wheatcroft 1990; Wheatcroft and Drake 2003). For example, when sedimentation is continuous over millennial scales (e.g. 1 m deposited over 10 000 years, leading to 0.01 cm a⁻¹, with four 5 cm-thick increments sampled every 20 cm) and the depth of mixing is limited (shallower than the thickness of 5 cm-thick increments), time averaging of fossil assemblages can be estimated as the transit time of skeletal particles through a 5 cm-thick increment, i.e. it will be equal to 500 years (Fig. 2a). In the same scenario, the stratigraphic resolution of two increments separated by 20 cm is equal to 2000 years. The same sediment column deposited over the total duration of 10 000 years can incorporate hiatuses exceeding 1000 years (Fig. 2b), leading to temporal separation between increments equal either to c. 3100 years (intervals with hiatuses) or c. 1100 years (intervals without hiatuses). Time averaging of increments between hiatuses will be lower (equal to c. 312.5 years) than in the first scenario because the short-term sedimentation is faster between the omission phases (0.016 cm a⁻¹). In this hiatus-rich scenario, the estimate of time averaging based on the long-term sedimentation rate overestimates the true time averaging (Fig. 2b). In contrast, if bioturbation is 10 cm deep (Fig. 2c), the time needed for shells to be advected below 5 cm is longer, and time averaging will be 1000 years (while stratigraphic resolution is reduced to 1500 years). In this case, the estimate of time averaging based on the long-term transit time underestimates its true magnitude (e.g. Tomašových et al. 2019a).

To summarize, a fossil record with a single long-term net sedimentation rate can be formed by assemblages differing in stratigraphic resolution and time averaging owing to variability in mixing depth and scale-dependency of sedimentation rates. In spite of these complexities, the estimates of time averaging based on the distribution of shell ages v. the estimates approximated from the long-term transit time can be an order-of-magnitude similar (Scarponi et al. 2013; Tomašových et al. 2022a). Simple age models for cores or outcrops thus still can be informative about the temporal scale of fossil assemblages even in the absence of postmortem age-frequency data. As the estimation of age distributions is increasingly more difficult in pre-Holocene successions, approximation of the temporal resolution of the stratigraphic record based on other sedimentological, taphonomic or ichnological criteria is a key step in palaeobiological analyses (Kidwell 1986; D.J. Davies et al. 1989b; Dominici 2004; Kemp and Sexton 2014; Crossley and Clark 2015; Trabucho-Alexandre 2015; N.S. Davies and Shillito 2021; Agiadi et al. 2022; Tomašových et al. 2022a).

The magnitude of time averaging in Holocene death and fossil assemblages can be directly assessed by dating multiple individual skeletal remains using U–Th, ¹⁴C and amino-acid methods (Kaufman and Manley 1998; Allen et al. 2013; Clark et al. 2014; Bright et al. 2021). However, most numerical estimates of time averaging are based on postmortem age distributions of a single species, and thus capture only within-species time averaging. When palaeobiological questions are directed towards the temporal fate of a single species in a region or over its larger geographical extent (Powell et al. 2017, 2020; Bergström et al. 2022), within-species time averaging is an appropriate target. However, multi-species or whole-assemblage time-averaging is of main interest in the studies of ecosystem dynamics that target changes in community composition or diversity. Only a few studies have assessed time averaging of two or more co-occurring species of molluscs and brachiopods (Kosnik et al. 2009, 2013; Krause et al. 2010; Tomašových et al. 2019a, b) or corals (Lybolt et al. 2011; Hammerman et al. 2021), and even less have targeted species markedly differing in durability or belonging to different phyla (Kowalewski et al. 2018; Albano et al. 2020; Nawrot et al. 2022). Although the outcomes of these studies are contingent on local conditions – in particular sedimentation rates and mixing depth – they typically detect significant age offsets between species and demonstrate that time averaging in multi-species assemblages can exceed time averaging of individual species. Death and fossil assemblages composed of a mixture of benthic taxa, typically represented by individuals inhabiting the seafloor at the scale of several tens of metres, and nektonic or pelagic taxa, with individuals distributed in the water column over larger regional scales, can further generate inherent difference in their spatial resolution. In this volume, Leonhard and Agiadi (2023)  discuss the formation and preservation of otolith death assemblages, with overview of applications of stable isotopes, elemental ratios or sclerochronological methods that are informative about the ecological dynamics of past fish populations.

Three types of changes in temporal scale can be expected to occur in analyses that use palaeoecological records or compare palaeoecological and ecological data: (1) changes in stratigraphic resolution and in time averaging within cores or outcrops; (2) changes in temporal scale at the transition between the surface mixed layer and the historical layers; (3) disparity between temporal scale of time-averaged palaeoecological data and snapshots of living assemblages represented by ecological samples.

1. Systematic changes in stratigraphic resolution and in time averaging within cores or outcrops result from a first-order control of sea-level change on sedimentation rate, leading to the formation of hiatuses or condensed sediments with high levels of skeletal alteration at maximum flooding surfaces (Kidwell 1986; Banerjee and Kidwell 1991; Benvenuti and Dominici 1992; Scarponi et al. 2017; Zecchin et al. 2021; Zimmt et al. 2022). Stratigraphic successions on siliciclastic continental shelves deposited during the Holocene thus typically exhibit significant increase in time averaging of fossil assemblages near the maximum flooding surfaces, followed by declining time averaging in the highstand systems tract under regimes characterized by sediment progradation (Scarponi et al. 2013, 2017). Trends in time averaging in Holocene carbonate sequences are expected to be more complex than in siliciclastic systems, as carbonate accretion rates depend not only on the accommodation space, but also on the type of carbonate factory and the associated trophic and oceanographic regimes that affect ecology of carbonate producers (Pomar and Haq 2016; Reijmer 2021).

2. Changes in temporal scale within the mixed layer and eventually between the mixed layer and the historical layers can be generated by three distinct processes. First, skeletal remains in the lower parts of the mixed layer are exposed to stochastic mixing induced by burrowers for longer, leading to higher time averaging near the base of the mixed layer relative to its uppermost part near the sediment–water interface (Tomašových et al. 2023). Second, the top-core increments in the mixed layer are frequently geologically transient, especially on shelves subjected to storms or other sources of erosion and winnowing, contributing to scale dependency of sedimentation rates (Sadler 1981; Frignani et al. 2005; Katz et al. 2022). This effect will also generate the smallest time averaging in the surface, geologically transient increments. Third, the overall downcore decline in porosity and water content owing to increasing sediment compaction is also expected to generate downcore increase in time averaging (Taranu et al. 2018; Stegner et al. 2019). In this volume, the downcore increase in time averaging was detected by Berensmeier et al. (2023)  in a study of a Holocene core in the Northern Adriatic Sea. They found that the degree of time averaging of assemblages formed by shallow-infaunal bivalves shifted from centuries in the uppermost increment in prodelta silts to millennia in the condensed shell bed corresponding to the transgressive sand sheet. As age-homogeneity is not developed in the surface 10 cm, as opposed to age-homogeneity in the subsurface shell bed, the authors hypothesize that this trend in time averaging was also enhanced by a recent decline in the depth of mixing related to the twentieth century eutrophication.

3. Time-averaged palaeoecological and non-averaged ecological samples used in live–dead analyses represent fundamentally different types of data in terms of their temporal resolutions (Kidwell and Bosence 1991; Martin 2004; Olszewski and Kidwell 2007; Tomašových and Kidwell 2009; Kidwell 2013). The contrast in temporal scale between death assemblages time-averaged to centuries or millennia (Flessa et al. 1993; Carroll et al. 2003; Kidwell et al. 2005; Kosnik et al. 2007) and instantaneous or analytically averaged living assemblages thus typically exceeds at least one order of magnitude (Kidwell 2013).

The ecological attributes that are compared between the assemblages in three sampling designs (Fig. 1) or in any time-series can be defined by species- or genus-level composition, functional traits, phylogenetic relatedness, body-size structure, or other aspects of pre- and post-impact ecosystems (Cintra-Buenrostro et al. 2005; Simões et al. 2009; Dietl et al. 2016; Harnik et al. 2017; Haselmair et al. 2021; Lawing 2021; Pruden et al. 2021; Agiadi et al. 2023). An assemblage that differs in composition from preceding assemblages in a temporal series (Tomašových and Kidwell 2011; Pandolfi et al. 2020; Mottl et al. 2021) by an amount that exceeds baseline variability is indicative of a major ecological shift and a novel state. However, the magnitude of the compositional shift is expected to increase with decreasing stratigraphic resolution and to decline with increasing time averaging. On the one hand, the magnitude of ecological change will be overestimated if gaps between fossil assemblages are not detected or when the composition of non-averaged living assemblages reflects a short-term variability rather than a genuine emergence of a novel state. Such changes in temporal scale of assemblages can lead to false positives (i.e. to detection of a novel state when none is present). On the other hand, the magnitude of ecological change will be underestimated when pre- and post-impact states are mixed in death assemblages, and can lead to false negatives when detecting novel states.

On the one hand, the differences in diversity between time-averaged death assemblages and snapshot-like living assemblages are expected to occur even in the absence of any actual temporal change in ecology of skeletal producers (Peterson 1977; Staff and Powell 1988; Olszewski and Kidwell 2007; Tomašových and Kidwell 2009). For example, live–dead analyses consistently found that local and regional inventory diversity (i.e. alpha and gamma diversity) increase, while compositional turnover between samples (beta diversity) decreases in death assemblages relative to the live communities from which they were sourced. These differences are not driven by differences in sample size between living and death assemblages (Tomašových and Kidwell 2009). The smaller beta diversity observed in fossil or death assemblages relative to living assemblages thus does not need to reflect any temporal change in compositional turnover (Kidwell and Tomašových 2013). These empirical results are qualitatively and quantitatively consistent with metacommunity models that predict changes in alpha and beta diversity only as a consequence of differences in time averaging between living and death assemblages. On the other hand, true changes in spatial or temporal turnover can be muted by high time averaging of fossil or death assemblages. The apparently low temporal variability in the composition of fossil assemblages frequently observed in deep-time and near-time stratigraphic sections, contrasting with higher variability observed among non-averaged living assemblages, can be partly explained by the effect of time averaging (Tomašových and Kidwell 2010). Anthropogenic impacts can homogenize spatial variability in species composition, reducing spatial beta diversity of living post-impact communities relative to their pre-impact counterparts. However, as beta diversity of death assemblages sourced by pre-impact communities will be also reduced by their time averaging, the estimates of beta diversity in live and death assemblages can be similar, masking the signature of recent biotic homogenization in live–dead analyses.

To control for these effects, one approach is to standardize temporal resolution by pooling assemblages into bins of equal duration and separated by equal time intervals. However, this binning approach can lead to a loss of stratigraphic resolution, which eventually becomes lower than allowed by age distributions of skeletal remains (Burke et al. 2019; Pandolfi et al. 2020; Mottl et al. 2021; Staples et al. 2022), and does not account for mixing of pre- and post-impact cohorts within bins. When scales cannot be standardized, the cross-scale comparisons (either between death and fossil assemblages within a core or between death/fossil and living assemblages) can still be informative about the magnitude of turnover when time averaging and stratigraphic resolution of the assemblages are estimated. Although many important ecosystem attributes such as diversity are highly sensitive to changes in temporal scale, as predicted by species–time relationships (Preston 1960; McKinney and Frederick 1999; Rosenzweig 2001; Adler and Lauenroth 2003), other characteristics such as mean species composition are less sensitive to scaling (Tomašových and Kidwell 2011). The challenging effects of variable temporal scale can be further resolved when age distributions of surface death assemblages or fossil assemblages in cores are available, as they (1) provide direct information on time averaging and stratigraphic resolution and (2) allow reconstructing the original time-series from a time-averaged stratigraphic record using stratigraphic unmixing.

The live–dead and dead–fossil approach can capture and identify the subsets of compositional shifts because the surface death assemblages are compared against distinct reference states, i.e. to older (fossil) assemblages in dead–fossil analyses or to younger (living) assemblages in live–dead analyses. To assess the effect of time averaging on the accurate detection of the magnitude of the past ecosystem shift and to confront the mixed nature of surface or top-core death assemblages, it is useful to map their time averaging against the timing of the major ecosystem shift, set to occur in AD 1950 in Figure 3. This mapping shows pros and cons of analyses that use the surface death assemblages either as tracers of initial states or as end-points of ecosystem change. The shift from pre-impact to post-impact states is buried deeply below the sediment–water interface when sedimentation rate is high (Fig. 3a). In this scenario, the within-core trend in composition (dead–fossil comparison) alone captures the transition from pre- to post-impact community state. The surface death assemblage does not incorporate any cohorts from pre-impact community states, and a live–dead comparison is thus not informative about the ecological shift (Fig. 3a). Surface death assemblages, which formed in deltaic environments with high sedimentation rate, can lack taphonomic inertia and conform to this first scenario. Such assemblages consist largely of cohorts a few years old and almost fully correspond to post-impact states. For example, the stratigraphic depth at which the mid-twentieth century regime shift is preserved is located 90 cm below the sediment–water interface in the Po prodelta (Tomašových et al. 2018). Under intermediate sedimentation rate (Fig. 3b), the within-core trend also captures the transition from a pre- to a post-impact community state but the top-core assemblage consists of a mixture of pre- and post-impact cohorts (moderate taphonomic inertia). Both the dead–fossil and the live–dead comparison underestimate the magnitude of the compositional shift in this second scenario (Fig. 3b). In contrast, the pre-impact community states remain exposed at the sediment–water interface under low sedimentation rate (Fig. 3c). In this third scenario, the dead–fossil comparison strongly underestimates the overall compositional shift because the taphonomic inertia is high and the surface death assemblage is thus strongly dominated by cohorts from the pre-impact state. As the top-core death assemblage effectively captures the pre-impact community state in this third scenario, live–dead analysis is more informative about the magnitude of the compositional shift induced by anthropogenic impacts than the dead–fossil analysis (Fig. 3c).

Time averaging of the fossil and surface death assemblages and the timing of the compositional shift relative to the timing of sampling thus represent the two key variables that determine whether the surface death assemblages are dominated by pre-impact and post-impact cohorts. The live–dead approach will be most efficient when surface death assemblages are highly time-averaged whereas the dead–fossil approach will be most effective when their time averaging is limited and sediment cores possess high stratigraphic resolution. This difference is taken into account in live–dead analyses where surface death assemblages are taken to be conservative tracers of anthropogenic impacts, i.e. the lack of live–dead mismatch can reflect the lack of inertia and does not indicate the absence of compositional shift, whereas a strong mismatch is a robust signature of recentmost anthropogenic impacts (Kidwell 2007). Taphonomic inertia in marine environments contrasts with palaeolimnological studies where surface lake sediments are sampled every few years to assess changes in the standing community, i.e. death assemblages are expected to track changes in source communities, and inertia is assumed to be negligible (Smol 1992).

In this volume, Albano et al. (2023)  documented very low time averaging of epifaunal bivalves on rocky substrates from the Mediterranean Israeli shelf and from seagrass beds (Posidonia oceanica) from the Crete shelf, with seven out of eight samples possessing decadal-scale interquartile age range. This limited time averaging contrasts with higher time averaging observed in soft sediments. The seagrass substrate is represented by the so-called matte formed by the dense system of roots and rhizomes of Posidonia oceanica, which entangles shells and limits bioturbation. Rocky substrates are probably affected by fast physical removal of shells and their transport to adjacent habitats. Surface death assemblages on rocky substrates and Posidonia matte substrates thus can be expected to be characterized by low inertia and low ability in capturing older community states.

However, taphonomic inertia is also a function of the number of individuals entering the death assemblage per time, and if production of shells in a new ecosystem state is small, like when total standing density of shell-bearing organisms declines, even death assemblages with limited time averaging can be inert as the new, post-impact cohorts are numerically rare. For example, Meadows et al. (2023)  assessed the efficiency of surface death assemblages to trace rapid climate change in Alaska on the basis of a 15-year time-series of living assemblages. They found that surface death assemblages are dominated by deposit-feeding nuculanid bivalves, whereas living assemblages are dominated by mixed-feeding tellinid bivalves. Although age data and skeletal preservation indicate very low, multi-decadal time averaging, and thus low potential for inertia (Meadows et al. 2019), live–dead mismatch is high at locations where the bivalve biomass and abundance has declined over time. This study thus underscores the role of shell input in enhancing the taphonomic inertia as death assemblages are not diluted by the recentmost shells sourced from post-impact living assemblages. Surface death assemblages from arctic environments thus represent a short-lived but accurate ecological memory of conditions prior to the recentmost warming.

Here, we discuss two approaches that can be used to account for the changes in temporal scale when inferring ecosystem dynamics at local or regional spatial scales. First, when scales cannot be standardized, we suggest sampling designs that combine dead–fossil (within-core) with live–dead analyses alone. These designs thus should integrate (i) fossil assemblages from cores or outcrops, (ii) surface death assemblages located within the mixed layer and (iii) living assemblages. Second, frequency distributions of postmortem ages collected from surface death assemblages or subsurface core increments allow scale standardization. Age distributions are shaped by disintegration in the TAZ, by mixing, by burial of skeletal remains below the TAZ by new sedimentation (or reef accretion), and by temporal changes in original input of skeletal remains from living assemblages to the mixed layer (Tomašových et al. 2014, 2023). The disintegration–burial dynamic predicts a declining abundance of increasingly older cohorts in surface death assemblages, which is thus a null expectation and not necessarily an ecological signal of recentmost increase in input of skeletal remains into a death assemblage. However, if this effect of disintegration and burial is accounted for, we suggest that age distributions from surface death assemblages or sediment cores can be used to infer temporal changes in abundance and mortality of skeletal producers.

Most surface death assemblages in marine systems, especially when deposited or time averaged over the duration that encompasses the ecological regime shift, will typically represent some mixture of pre- and post-impact states. Therefore, the combination of subsurface fossil assemblages (not affected by young cohorts from post-impact states), surface death assemblages (mixtures of pre- and post-impact states), and post-impact living assemblages can be most informative in analyses of baselines and novel states. Several studies effectively used this approach (Edinger et al. 2001; Alin and Cohen 2004; Schönfeld and Altenbach 2005; Mendes et al. 2013; Casey et al. 2014; Kusnerik et al. 2020). On the one hand, the increments that are sufficiently deep below the present-day mixed layer may cover the baseline states better than the surface death assemblages, and thus sediment cores penetrating beyond the mixing depth can be most efficient in capturing pre-impact states. Studies that use the Holocene outcrops are thus also not affected by inertia. On the other hand, living assemblages are not affected by any mixing, and when sampled through time or over a broader area, can better capture post-impact states than uppermost increments of sediment cores alone (Tomašových et al. 2019b).

To document the power and efficiency of this integrated approach, we visualize the shift in molluscan community composition that took place in the late twentieth century in the Gulf of Trieste (northern Adriatic Sea): the shift is visible but muted when comparing subsurface fossil assemblages with top-core death assemblages (Gallmetzer et al. 2017; Tomašových et al. 2017), whereas the addition of living assemblages illuminates the actual magnitude of the compositional change. The multivariate analysis of two 1.5 m-long cores from the Bay of Panzano in the Gulf of Trieste collected at 12 m water depth in 2013 compared with surveys of living assemblages collected with Van Veen grabs between 1985 and 2013 at water depths between 10 and 30 m visualizes the integration of all three types of records (Fig. 4). This survey is comparable to the intermediate scenario in Figure 3, with decadal time averaging in the surface death assemblages, centennial time averaging in the subsurface fossil assemblages and long-term sedimentation rate of c. 0.2–0.4 cm a⁻¹ (Tomašových et al. 2017, 2018). A compositional shift occurring in the twentieth century preserved c. 15 cm below the sediment–water interface, is muted relative to the shift observed between the surface death assemblages and living assemblages (Fig. 4). The total compositional shift towards high dominance of the hypoxia-tolerant bivalve Corbula gibba (from c. 20 to 70%) is instead best documented by the difference between subsurface (fossil) and living assemblages. The full integration that incorporates data from cores and living assemblages thus shows that the magnitude of the compositional shift is higher than expected on the basis of within-core analysis or dead–fossil analyses alone.

Similarly, the sampling design in Poirier et al. (2022)  incorporates different types of records, including living assemblages from tidal flats, surface death assemblages from active beach ridges and fossil assemblages from older chenier ridges in the southern English Channel. Owing to depositional conditions that led to migration of chenier ridges, these death and fossil assemblages provide a regional-scale, time-averaged signal of the subtidal and intertidal ecosystems, whereas living assemblages are informative about present-day intertidal communities. They found that the back-barrier fossil ridges that were largely formed during the Early Medieval Period document pre-industrial baseline conditions. Death assemblages from active beach ridges are characterized by millennial-scale time averaging, and thus are dominated by older cohorts that swamp the signature of the twentieth century colonization by invasive species and by species associated with shellfish farming. As active beach ridges are not strongly diluted by recentmost input of subtidal species, differences in composition between fossil and active ridges record a decline in grazing and carnivorous species related to the loss of seagrass beds and macroalgal belts occurring over the past centuries. Poirier et al. (2022)  suggest that this change reflects warming that followed after the Little Ice Age cold phase and the recent anthropogenic warming.

The time-scales exhibited by the stratigraphic record are ultimately determined by sedimentation, mixing and disintegration rate of skeletal remains in the mixed layer and by temporal variability in these processes. A transformation of stratigraphic depth to time typically involves squeezing or stretching distances among increments as informed by an age model (Ripepe and Fischer 1991; Meyers 2012; Li et al. 2018; Hohmann 2021; Lougheed 2022). However, as shown in Figure 1, first, interpolation of ages to undated increments is sensitive to the mismatch between long-term and short-term sedimentation rates, i.e. short-term variability in sedimentation rate is missed by interpolation between levels that are directly age-dated (Sadler 1981; Schumer et al. 2011; Tipper 2015; Trampush and Hajek 2017). Second, a single age model characterized by a single long-term net sedimentation rate can be associated with distinct stratigraphic resolution and increment-specific time averaging owing to differences in the depth of biogenic mixing. In high-resolution studies, deconvolution approaches can transform stratigraphic patterns to temporal patterns by incorporating the structure and depth of the mixed layer (Berger and Heath 1968; Guinasso and Schink 1975; Schiffelbein 1985; Dolman et al. 2021; Liu et al. 2021). These approaches either assume that the mixed layer is vertically homogeneous or require parameterization of the downcore profiles in mixing rates. They do not account for non-local mixing and do not resolve any ecological changes that are averaged out within the increments themselves. On the other hand, ecological attributes of individual skeletal particles from the Holocene sediments can be ordered along a temporal axis based on their absolute ages estimated with amino-acid racemization, ¹⁴C or U–Th, regardless of their position in the stratigraphic column (Tomašových et al. 2017, 2019a, b). Increment-specific age distributions obtained in this way thus not only provide estimates of time averaging and insights into disintegration-burial dynamics, but also provide key insights into the depth-to-time transformation and the ecosystem history of carbonate producers (Toth et al. 2012; Clark et al. 2017; Leonard et al. 2020a, b).

Assuming a steady-state input of dead remains into a surface death assemblage, frequency distributions of postmortem ages are shaped by disintegration of remains in the TAZ, by burial and exhumation induced by burrowers in the mixed layer, and by burial of skeletal remains by new sedimentation from terrigenous input or by reef accretion. Under such conditions, age distributions of skeletal remains in surface death assemblages sourced by active constant input from living populations, will be right-skewed, i.e. dominated by the youngest cohorts. Their shapes will follow an exponential distribution (i.e. the slope of a distribution follows a straight line when frequencies of age cohorts are log-transformed) when disintegration rate in the TAZ does not change with shell age. The shapes of age distributions will be more complex, heavy-tailed or L-shaped when disintegration rate declines with shell age and/or older skeletal remains are preferentially exhumed to the TAZ from the underlying, shell-rich zones (Tomašových et al. 2014, 2023). Under constant input of dead individuals into surface death assemblages, the mode of such distributions will be formed by the recentmost cohorts because disintegration rate in the TAZ, burial below it, and transport to other habitats, act longer on older cohorts (Flessa et al. 1993; Olszewski 1999, 2004). The right-skewed age distributions thus do not simply signify increasing abundance or mortality towards the time when the given death assemblage is sampled. Although the rate of input of shells from living to death assemblages is rarely constant and any stochastic oscillation in standing abundance or in timing of mortality will generate deviations from right-skewed age distributions, the overall tendency for higher abundance of younger cohorts can be expected when disintegration rates in the TAZ and/or burial rates below it are not negligible. Several analyses of coral production or mortality over the past decades or centuries suggested recent changes in these parameters on the basis of age distributions of surface death assemblages as these are increasingly dominated by young cohorts (Yan et al. 2019; Chen et al. 2021; Hammerman et al. 2021). Although these studies also used other lines of evidence in their overall inference about the recentmost ecological dynamic of coral communities (e.g. direct documentation of bleaching events), they did not account for the disintegration-burial effect in the assessments of age distributions in surface death assemblages. The importance of this effect depends on the magnitude of the rates of burial and disintegration. When these rates are negligible, the increase in abundance of youngest cohorts may reflect a true increase in abundance and/or mortality. In real-world scenarios, however, some combination of ecological variability and non-negligible disintegration and burial rate can be expected. Some independent information about the overall loss of older remains via disintegration, accretion or transport is thus of key importance in interpreting the overall shape of surface age distributions in palaeoecological studies. For example, when net sedimentation rate is very low as estimated by other means, and thus the residence time of skeletal remains in the mixed layer is primary determined by disintegration, the burial term can be neglected.

Estimates of disintegration rates are not only useful to understand the shape of age distributions but also because increasing rates of bioerosion and dissolution driven by eutrophication and ocean acidification compromise coral-reef functioning and can lead to net reef erosion (Andersson and Gledhill 2013). For example, the present-day accretion rates of shallow fore-reef coral habitats in the Caribbean Sea or in the Tropical Eastern Pacific can be an order of magnitude lower than Holocene net accretion rates (Perry et al. 2013; Wizemann et al. 2018; McNicholl et al. 2020; Enochs et al. 2021; Morris et al. 2022). Increased coral rubble cover can negatively impact coral recruitment, and rubble mobilization abrades and smothers corals (Kenyon et al. 2022). These processes leading to disintegration and erosion leave signatures not only in the shape of age distributions, but also in preservation states of skeletal remains themselves (Lescinsky et al. 2012; Dawson et al. 2014; Wizemann et al. 2015; Kuffner and Toth 2016; Ramírez-Viaña et al. 2021). Edelman-Furstenberg (2023)  found a geographical gradient in molluscan preservation that ranges from oxygen-depleted, upwelling environments that limit metazoan predators and borers, minimizing surface shell loss, to well-oxygenated environments with higher abundance of predators and microboring animals that all foster disintegration. Molluscan taphofacies with unusually low levels of skeletal damage from borers and encrusters can represent sensitive tracers of conditions affected by Anthropocene deoxygenation (Tomašových et al. 2021).

To reconstruct the initial number of individuals that entered the death assemblages, the effect of overall loss of dead remains on the final shape of age distributions can be accounted for by dividing the empirical number of individuals in age bins by the survival function of a disintegration–burial model. For example, when loss rate produced by disintegration (in the TAZ) and burial (below the TAZ) is constant in time and does not change with postmortem age, instantaneous rate of loss (λ) of skeletal remains is equal to the inverse of the mean of their age distribution (simple exponential model in Fig. 5a). The survival function of this simple exponential model is equal to e^−λ*age, and λ is a sum of disintegration and burial. When loss rate declines with postmortem age of skeletal remains, the survival function depends on the mechanisms of slowdown in disintegration. For example, the sequestration model (Fig. 5b, Tomašových et al. 2014) includes three parameters: (1) initial disintegration rate λ₁ in the TAZ, (2) sequestration rate τ at which disintegration declines and skeletal remains move from the TAZ into the so-called sequestration zone (SZ), and (3) loss rate λ₂ related to disintegration rate in the SZ and to burial into underlying increments (Fig. 5b). This model successfully explains the shape of age distributions formed by shells of molluscs in surface death assemblages (Tomašových et al. 2014). The sequestration parameter τ can mechanistically correspond to net rate of burial of skeletal remains below the TAZ by new sedimentation (or by burrowers, left column in Fig. 5b) and/or to rate of diagenetic stabilization (right column in Fig. 5b). The survival function of this sequestration model is equal to αe^−λ₂*age + (1−α)e^{−(λ₁+τ) ×age}, where α = τ/(τ + λ₁ + λ₂). As this sequestration model allows for zero disintegration rate of skeletal remains in the increments below the TAZ, it is also useful in correcting for disintegration in unmixing of age data from both surface and subsurface increments of sediment cores. This model can be extended from two zones (TAZ and SZ) to a stochastic transition-rate matrix, allowing explicit addition of mixing rates and prediction of downcore changes in the shape of age distributions (Tomašových et al. 2023).

In Figure 6, we show three age distributions based on surface coral death assemblages with right-skewed shapes and a tendency towards lower abundance of older cohorts that exceed several decades or centuries in age (grey histograms in Fig. 5). The highest abundance of the youngest cohorts in surface death assemblages is a null expectation of disintegration and burial dynamic under a temporally constant input of dead remains to the mixed layer, regardless of whether disintegration does or does not change with postmortem age of skeletal remains. The proportion of preserved coral remains decreases with the age of a cohort because older remains are exposed to disintegration for longer. The input of old cohorts thus would be equal to the input of young cohorts when disintegration of coral remains in the TAZ occurs at decadal scales in the Saudi Arabian Red Sea or at the Weizhou Island, or when it would be even slower at the Luhitou Peninsula (grey lines representing initial abundance of corals in Fig. 6). Although coral disintegration rates are geographically variable (Brown et al. 2021), they tend to be fast in the TAZ, with disintegration occurring at monthly and yearly scales (Molina-Hernández et al. 2022; Morais et al. 2022). Therefore, the contribution of disintegration to the abundance of cohorts as observed in age–frequency distributions should not be neglected.

The rarity of older cohorts, as observed in the Red Sea or South China Sea, thus does not necessarily imply any major recent increase in mortality of corals relative to former time intervals. The argument invoking the ecological decline in abundance of skeletal producers, based on the declining abundance of older cohorts in age distributions observed in surface death assemblages, thus needs to be validated on the basis of other evidence. For example, El Niño bleaching events that negatively affected the growth of massive corals and induced the shift towards rubble-dominated environments were documented in the Red Sea over the past three decades (Cantin et al. 2010; Hammerman et al. 2022), and monitoring studies also documented ecological decline in abundance of Acropora since 1980s and growth rate of Porites significantly declined in the late twentieth century in the South China Sea (Yu et al. 2019; Kang et al. 2021).

With some caveats, internal modes in age distributions of surface death assemblages can be diagnostic of past maxima in abundance or mortality of skeletal producers over the duration of time averaging in the mixed layer (Lybolt et al. 2011; Tomašových et al. 2016). If abundance of carbonate producers declines through time towards present-day sampling at a slow rate, the timing of maximum abundance preserved in the resulting age distribution will still be pulled towards the sampling time by the effect of the disintegration rate (Tomašových et al. 2016). If carbonate producers become locally extirpated, or mortality of long-living taxa such as massive corals stops, no new dead skeletal remains enter into a death assemblage and the pulling effect of the disintegration rate will be inactive. Therefore, major fluctuations in abundance occurring rapidly can have unusually large potential to be preserved in surface death assemblages. Multiple analyses of molluscs, brachiopods or corals show that surface postmortem age distributions are frequently multimodal and reflect variability in production at decadal, centennial and millennial scales. For example, surface age distributions document high volatility in bivalve abundance (Albano et al. 2016) and mortality rate of shallow-water corals at decadal scales over the past century (Clark et al. 2017; Chen et al. 2021; Hammerman et al. 2021), millennial-scale switches and shutdowns in production (Yan et al. 2019), or pervasive regime shifts in functional composition on benthic communities lasting for more than c. 100 years on the California shelf (Tomašových and Kidwell 2017) and in the Gulf of Trieste (Tomašových et al. 2019b).

Age distributions formed by cold-water corals in bathyal environments show internal modes in abundance separated by gaps that can exceed c. 10 000 years in duration (Fig. 7). This extensive time averaging, exceeding 10 000 years, indicates that both long-term net sedimentation rates and disintegration rates of aragonitic corals are very low in these environments. Disintegration models that invoke sequestration of skeletal remains via diagenetic stabilization can be more realistic than simple disintegration models in these environments. Incorporating this sequestration dynamic into the reconstruction of the original skeletal input can subdue or minimize the importance of the recentmost peak in abundance (Fig. 7). In this volume, Tomašových et al. (2022b)  documented a bimodal age distribution formed by the epifaunal brachiopod Gryphus vitreus at mid-bathyal depths, indicating two main, temporally narrow peaks in abundance of this species over the past two millennia. The high abundance of this species in the past contrasts with its present-day rarity in the southern Adriatic Sea. Excellent preservation of still articulated shells that are several centuries old indicates that their disarticulation and overall disintegration rates are very low. Thus, the difference in shape between the raw age distribution and the reconstructed trajectory in abundance of this brachiopod is negligible.

Several studies also found remarkable hiatuses in production of shallow-water corals over the past millennia as observed in age distributions of fossil assemblages extracted from sediment cores. Such hiatuses in coral production may still be produced by temporal changes in preservation. Rodriguez-Ruano et al. (2022)  estimated that it would take 135–800 years for bioerosion to remove 2300 years-worth of reef framework growth. However, reef framework in Panamá is stabilized by sediment that prevents the activity of bioeroders; only the upper c. 1 m of the open framework produced during several decades was vulnerable to erosion, supporting the hypothesis that increased ENSO-driven climatic variability was responsible for the hiatus in reef framework off Pacific Panamá.

The conditions that determine the transition of surface death assemblages through the mixed layer and into the historical layers (i.e. addition of new sediment, mixing, and disintegration) influence the shapes of age distributions in successively deeper increments. The age distributions that are right-skewed in a surface increment become more symmetric in subsurface increments owing to mixing effects, with the mode of distributions defined by the rate of long-term net sedimentation or accretion (Tomašových et al. 2023). The controls of the shape of fossil age distributions are thus more complex than those of the shape of age distributions from surface death assemblages, although parameters such as disintegration, burial and mixing can still be estimated with stochastic transition-rate matrices (Tomašových et al. 2023). Interpreting age distributions over the duration of core deposition is not based on the shape of the age distribution in any specific increment, but on the shape of the whole-core age distribution. When sedimentary increments consist of time-averaged assemblages formed by overlapping age cohorts, stratigraphic unmixing informed by a succession of increment-specific age distributions can transform the stratigraphic trend in preserved abundance to the original chronological signal in species abundance (Tomašových et al. 2017, 2019a, b).

The unmixing procedure consists of four steps. We show these steps using the stratigraphic record of the infaunal bivalve Timoclea ovata collected in a 1.55 m-long core at 44 m water depth in the NE Adriatic Sea off Brijuni Islands (Schnedl et al. 2018; Gallmetzer et al. 2019). The mixed layer is about 20 cm thick and the downcore increase in median age of the bivalve indicates that the long-term (millennial-scale) net sedimentation rate was very slow at the coring site (<0.01 cm a⁻¹, Tomašových et al. 2022a). Age distributions in 5 cm increments show that both death assemblages of this bivalve in the mixed layer and subsurface assemblages below it are time-averaged to a few millennia. This species is numerically abundant in the middle and lowermost part of the core (Fig. 8a). The unmixing procedure starts with: (1) estimation of age distributions in 13 increments on the basis of age data, using amino acid racemization calibrated with ¹⁴C (n = 300, with 20–30 individuals dated per increment); (2) interpolation of the age distributions (Fig. 8b) to undated increments on the basis of adjacent increments with dated shells (this interpolation step can be justified when subsurface age distributions are symmetrical and normal-shaped); and (3) inferring the final shape of the distribution in each increment by sampling N shell ages from each increment-specific, directly observed or interpolated age distribution (N = total abundance of T. ovata in each increment) (Fig. 8c). This procedure thus leads to the whole-core age distribution represented by all shells of T. ovata preserved in the core (grey histogram in Fig. 8d). In the fourth step, similarly as in examples of age distributions from surface death assemblages, the unmixing proceeds with the estimation of disintegration rate in the mixed layer (Fig. 8e). The effect of burial below the core can be neglected because its base is marked by a ravinement surface and no shells of T. ovata are located below it (Fig. 8f).

The age distribution in the uppermost 12 cm is L-shaped and fits well with a sequestration model where disintegration rate declines with shell postmortem age. Although a simple exponential model (Fig. 5a) can be useful when accounting for a temporally-constant disintegration of skeletal remains in surface death assemblages, or their burial below them, the assumption of temporal constancy in disintegration is violated in sediment cores that encompass increments below the TAZ. Therefore, the sequestration model (Fig. 5b), with three parameters that specify rates of disintegration within and below the TAZ and rate of their sequestration below the TAZ, provides a natural model for unmixing of abundance variability of skeletal producers from sediment cores. The sequestration-model estimate of rate at which disintegration declines (τ = 0.0013 a⁻¹) is comparable to the estimate of the overall net sedimentation rate. The division of the whole-core distribution with the survival function of the sequestration model (disintegration rate [λ₁] = 0.036 a⁻¹, sequestration rate [τ] = 0.0013 a⁻¹, and final disintegration rate or burial below the core [λ₂] is set to zero) removes the recentmost peak in abundance and leads to the reconstructed variability in abundance of T. ovata in 200-year cohorts (Fig. 8g). The estimates of abundance in these age bins are a function of population density and lifespan of the species. If lifespan can be estimated, these values can be scaled to yearly standing densities (Tomašových et al. 2017). The unmixing procedure in the R language script is shown in the Supplementary material.

In the absence of unmixing, the raw age distribution, which is constrained by the number of dated specimens and thus does not account for stratigraphic variability in abundance (white histogram in Fig. 8d), would miss the major peak in abundance of T. ovata in the past (Fig. 8e). The peak may represent habitat tracking, i.e. the spatial shift of a population to preferred habitats in response to environmental changes (Brett et al. 2007), as the increase in abundance from 8000 to 6000 years ago corresponds to the phase when sea level rose rapidly and reached the present-day position in the northern Adriatic Sea. The ecological cause of the abundance decline of T. ovata starting at 4000 years ago may be related to the development of the bryozoan meadows, which are typical of the eastern Adriatic late Holocene (McKinney 2003).

Other unmixing analyses documented significant fluctuations in production of benthic molluscs. For example, shallow-infaunal bivalves Gouldia minima and Nuculana taphria significantly declined in abundance over the past one or two centuries in the northern Adriatic Sea (Tomašových et al. 2019b) and on the southern California shelf (Tomašových et al. 2019a), respectively. Whole-core age distributions also constrain the timing of their abundance increase that coincides with palaeowater depths at the coring sites corresponding to modern bathymetric preferences of these subtidal species. Analyses of age distributions of shallow-water corals preserved in sediment cores indicate local- to regional-scale turnoffs in production that lasted for few millennia in the tropical Eastern Pacific (Toth et al. 2015, 2019), at the Great Barrier Reef (Perry and Smithers 2011; Leonard et al. 2020a, b), and in the Florida Keys (Toth et al. 2018). In the absence of age distributions, the detection of these dynamics on the basis of the raw stratigraphic record in molluscan or coral abundance would be limited owing to condensation and mixing affecting these cores.

The inferences of past ecological dynamics on the basis of the fossil record are inherently constrained by (1) the differences in the temporal scale of fossil assemblages in sediment cores or outcrops, death assemblages in the surface mixed layer and living assemblages in ecological time-series, and (2) the challenges in transformation of depth to time owing to stratigraphically variable temporal resolution and scale dependence of sedimentation rates. Stratigraphic patterns rarely accurately mirror chronological trajectories, and when such perfect fidelity is assumed, palaeoecological inferences will lead to ill-informed interpretations and thus incorrect recommendations for conservation and restoration of marine ecosystems. Whether surface death assemblages capture predominantly pre-impact or post-impact community states depends on sedimentation, mixing and disintegration in the mixed layer and on the timing of the ecological shift relative to the time of sampling. On the one hand, surface death assemblages do not necessarily capture baseline states accurately because they are still affected by the input of recentmost skeletal remains from the novel community states. On the other hand, within-core analyses comparing subsurface fossil assemblages with surface death assemblages can also underestimate the magnitude of ecological shift because the time-averaged death assemblages can be dominated by older cohorts, thus not accurately capturing post-impact states. However, integrating surface death assemblages with both fossil assemblages and living assemblages can resolve these challenges and untangle the dynamics of past ecosystem shifts. The assessment of postmortem age distributions allows direct estimation of stratigraphic resolution and time averaging. Even more importantly, age distributions can be used to transform time-averaged stratigraphic records to chronological palaeoecological patterns, leading to findings of high volatility in species abundance, long-term oscillations in production or detection of abrupt regime shifts.

We thank Moriaki Yasuhara and Mark Whiteman for scientific and editorial comments. This research contributes to the objectives of Q-MARE (a PAGES working group).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

AT: conceptualization (lead), writing – original draft (lead), writing – review & editing (equal); SD: conceptualization (equal), writing – original draft (supporting), writing – review & editing (equal); RN: conceptualization (equal), writing – original draft (supporting), writing – review & editing (equal); MZ: conceptualization (equal), writing – original draft (supporting), writing – review & editing (equal).

This research was supported by the Slovak Research and Development Agency (grant no. APVV 17-0555), and the Slovak Scientific Grant Agency (grant no. VEGA 2/0106/23).

The data and scripts are available in the Supplementary material.