A new approach for investigating spatial relationships of ichnofossils: a case study of Ediacaran–Cambrian animal traces

Abstract

Trace fossils record foraging behaviors, the search for resources in patchy environments, of animals in the rock record. Quantification of the strength, density, and nature of foraging behaviors enables the investigation of how these may have changed through time. Here, we present a novel approach to explore such patterns using spatial point process analyses to quantify the scale and strength of ichnofossil spatial distributions on horizontal bedding planes. To demonstrate the utility of this approach, we use two samples from the terminal Ediacaran Shibantan Member in South China (between 551 and 543 Ma) and the early Cambrian Nagaur Sandstone in northwestern India (between 539 and 509 Ma). We find that ichnotaxa on both surfaces exhibited significant nonhomogeneous lateral patterns, with distinct levels of heterogeneity exhibited by different types of trace fossils. In the Shibantan, two ichnotaxa show evidence for mutual positive aggregation over a shared resource, suggesting the ability to focus on optimal resource areas. Trace fossils from the Nagaur Sandstone exhibit more sophisticated foraging behavior, with greater niche differentiation. Critically, mark correlation functions highlight significant spatial autocorrelation of trace fossil orientations, demonstrating the greater ability of these Cambrian tracemakers to focus on optimal patches. Despite potential limitations, these analyses hint at changes in the development and optimization of foraging at the Ediacaran/Cambrian transition and highlight the potential of spatial point process analysis to tease apart subtle differences in behavior in the trace fossil record.

Introduction

Trace fossils capture the activity of extinct organisms, providing important insights into the evolution of behaviors such as burrowing and foraging (e.g., Plotnick 2012). Foraging animals develop strategies that minimize their energy expenditure when searching for resources, which are not typically distributed uniformly (Bond 1980; Brown 2000). Optimization depends on the extent of resource heterogeneity—in a resource-rich area, suboptimal searching strategies incur little or no cost (Codling et al. 2008). Identifying the extent of heterogeneity in trace fossils can provide important context regarding the behavioral ecology of the tracemakers through geological time, as well as the distribution of resources in a particular environment. Current metrics for describing horizontal trace fossils are primarily limited to the bedding-plane bioturbation index (BPBI), which describes the extent of bioturbation through assignment of a score from 1 to 5, where 1 indicates no bioturbation and 5 indicates maximal bioturbation (Miller and Smail 1997). Further work to estimate the percentage of bioturbation on a surface uses grid maps and the presence/absence of trail intersections (Marenco and Hagadorn 2019). However, these metrics do not quantify the relative patchiness of ichnofossils or provide the methodological power to compare to random distributions to a suite of other models.

One alternative approach to investigate such patchiness is to use spatial point process analyses (SPPA) to quantify the spatial scale and strength of spatial patterns with respect to the positions of horizontal trails on bedding planes (Illian et al. 2008; Wiegand and Moloney 2013). SPPA in the fossil record have focused primarily on the use of nearest neighbor analyses (e.g., Leighton and Schneider 2004; Shroat Lewis et al. 2014), where the distances between each pair of fossil specimens are calculated and compared with expected random patterns to determine whether nonrandom spatial patterns occur. However, nearest neighbor analyses only capture small spatial scales; if all specimens are within 10 cm of each other, for example, then the maximum resultant spatial pattern described will be 10 cm. This spatial limitation means that it is not possible to capture patterns at different spatial scales, so complex patterns are often not detected. Recently, SPPA have been successfully applied to body fossils to elucidate Ediacaran reproduction (Mitchell et al. 2015), competition dynamics (Mitchell and Kenchington 2018), taphomorph identification (Mitchell and Butterfield 2018), morphological characters (Mitchell et al. 2018), and the drivers of community dynamics (Mitchell et al. 2019, 2020), as well as facilitation within Silurian coral communities (Dhungana and Mitchell 2020) and Jurassic crinoid colonization (Hunter et al. 2020). Such studies were possible due to an extensive set of ecological models that, in modern settings, correspond to different underlying processes for sessile organisms, such as mutual habitat associations (Getzin et al. 2008), dispersal limitation (Lin et al. 2011), facilitation (Lingua et al. 2008), and density-dependent mortality (Fey et al. 2019). These models enable comparisons between the distributions of in situ fossil communities with those controlled by known ecological and biological processes today to determine the most likely control(s) of their spatial patterning.

Few studies have used SPPA to investigate the distribution of trace fossils. SPPA of the spatial positions of the Ediacaran body fossil Kimberella and its trace scratch marks Kimberichnus are suggestive of avoidance of pre-grazed patches (Mitchell et al. 2020). SPPA have also been used to visualize and quantify the distribution of predatory traces on shelled invertebrates, enabling the resolution of different modes of drilling predation (Rojas et al. 2020). In principle, this approach can be extended to trails, trackways, and burrows, whereby the trace fossils are described by line segments (Fig. 1). The SPPA of their midpoints can be used to determine the nature (attractive/clustering or repelling/segregation), spatial scale, and strength of certain behaviors. Unfortunately, work using SPPA on movement and/or foraging ecology has not been used to investigate the trajectories of extant or fossil organisms, so there are few models for comparison with trace fossil spatial patterns. This limited set of movement models means that SPPA modeling fitting comparisons cannot be performed, and instead the SPPA of trails and trackways focus on describing relative lateral heterogeneity and optimization of the foraging strategies employed. For example, the trace midpoints of a single randomly moving organism (Fig. 1A) will produce a significantly different spatial point pattern than an organism responding to spatial heterogeneity (Fig. 1B,C). In the latter case, there are points that are both closer together and further apart than the random situation, that is, a higher variance of interpoint distances. Nonrandom patterns may indicate the ability to focus on good-quality resource patches. Comparisons of the spatial distributions of different trace midpoint patterns on a single surface can further test whether tracemakers responded to similar underlying heterogeneities, because they would have focused on the same area/patches, thus creating traces in this same area (Fig. 1). Therefore, two tracemakers that respond to the same heterogeneities will result in midpoint spatial patterns of similar scale and nonrandom bivariate spatial distributions. Such analyses also provide insights into possible time averaging: substrate heterogeneities are unlikely to persist for longer than ecological timescales, so nonrandom spatial patterns between ichnotaxa suggest that they co-existed within similar timescales (cf. Mitchell et al. 2015).

Figure 1. Schematic drawing showing different movement paths (blue lines), midpoints of the movement paths (stars), and background patchy environment (gray areas). Each midpoint is identified by the length and orientation of the corresponding line segments (“marks”). A, Random MCF: The moving organism does not respond to the underlying environment. The lengths of each segments, their orientations, and the distances between their midpoints are random. B, Nonrandom length movement: The tracemaker detects a preferable underlying substrate and so makes shorter path segments to maximize the time spent in the preferred resource. Marked lengths of paths and distances between midpoints are nonrandom, whereas marked orientations of paths are random. Short segments occur close together. C, Nonrandom length and nonrandom orientation MCF: The tracemaker detects a preferable underlying substrate and makes both shorter path segments and sharper turns to maximize the time spent in the preferred resource. Crucially, it turns around when it leaves its preferred area. Mark lengths and orientations of paths are nonrandom. Short segments and those with tight turns are more likely to occur closer together.

Mark correlation functions (MCFs) are a suite of functions that quantify how the marks (here segment lengths and absolute orientations) of a point (here trace fossil midpoints) vary with spatial scale (Velázquez et al. 2016). Under an optimal foraging strategy, a tracemaker will adjust both the length and orientation of each segment to focus on a preferred resource (Fig. 1C vs. Fig. 1A,B). If it strays out of the resource area, it quickly turns back, resulting in high directionality of orientations (Buhler and Grey 2017). Therefore, the directionality of trace fossil orientations is expected to be spatially heterogeneous, mirroring the underlying resource, and accompanied by shorter segment lengths to maximize time spent in favorable resource patches (Kitchell 1979; Koy and Plotnick 2007, 2010). The optimization of foraging requires two different aspects: the ability to find quality resources, assessed using midpoint SPPA; and the ability to focus on these, for which MCFs are used. Combining these methods provides a novel approach for identifying the extent of foraging optimization of tracemakers.

To explore the power of SPPA methodologies to identify differences in behavior, we have applied them to two ichnofossil-rich horizons, from a critical time in the evolution of animal mobility—the Ediacaran/Cambrian transition. Trace fossils from this interval reveal crucial insights into the development of mobility, novel feeding behaviors, and the putative appearance of crown-group metazoans (Seilacher et al. 2005; Chen et al. 2013, 2019; Buatois and Mángano 2016; Evans et al. 2019; Xiao et al. 2019). Indeed, adaptions to optimize foraging strategy—defined here as the exploitation of any resource distributed on the seafloor—have been suggested as potential drivers of early bilaterian evolution and diversification (Koy and Plotnick 2007, 2010; Budd and Jensen 2017; Mitchell et al. 2020). As such, methods that enable the determination of the extent of lateral resource heterogeneities and analyses of how animals exploit such resources are crucial for understanding the adaptation and progression of foraging behaviors during the Ediacaran and into the Cambrian.

Here, we demonstrate how spatial analyses of trace fossils can be used to elucidate foraging behavior by analyzing morphologically simple horizontal trails on a bedding plane from the late Ediacaran Shibantan Member (Xiao et al. 2021) and another from the early Cambrian Nagaur Sandstone of the Nagaur Group in Rajasthan, NW India (Pandey et al. 2014). These horizons yield trace fossils with comparable gross morphologies (and therefore likely underlying behavior) making them ideal candidates for this comparative study and comparison with other terminal Ediacaran traces (Jensen and Runnegar 2005).

Geological Setting

The Shibantan Member of the upper Ediacaran Dengying Formation, Yangtze Gorges area, South China, contains a high abundance of ichnofossils. These are preserved within 100–150 m of dark gray, thin-bedded, laminated peloidal limestone deposited in subtidal environments above the storm wave base and regularly impacted by current activity (Meyer et al. 2014; Xiao et al. 2020, 2021). Depositional age is constrained between 550 and 543 Ma (Dornbos et al. 2004; Condon et al. 2005; Xiao et al. 2021; Yang et al. 2021; Fig. 2). Shibantan ichnofossils document extensive trails, resting or dwelling traces, undermat mining for oxygen and/or food (Xiao et al. 2019), and putative trackways of bilaterian animals with paired appendages (Xiao et al. 2021). Bedding-plane bioturbation intensity reaches 20%–40% (Meyer et al. 2014), comparable to some pre-trilobite Cambrian carbonates (Dornbos et al. 2004), although lower than the maximum observed for post-trilobite Cambrian carbonates (Fan et al. 2021). These ichnofossils are interpreted as recording the activity of organisms burrowing in and out of microbial mats, likely mining for oxygen and/or food (Xiao et al. 2019), which we refer to broadly as “resources,” and may have been spatially heterogeneous (Darroch et al. 2021) and temporally dynamic in Ediacaran oceans (Wood et al. 2015). The bedding plane analyzed in this study was excavated at a horizon ~70 m above the base of the Shibantan Formation, covering an area of 0.51 m². The excavated slabs were flipped and reassembled to view both part and counterpart (Droser et al. 2019).

Figure 2. Geological map of the Yangtze Gorges area (A) and stratigraphic column of the Dengying Formation (B). Inset map shows major tectonic units, with red rectangle marking area shown in A. Star marks field site. Fossil ranges are updated from Xiao et al. (2021) and Liang et al. (2020). The two possible correlations of the 551.1 ± 0.7 Ma tuffaceous bed are based on An et al. (2015) and Zhou et al. (2017) and are denoted by two blue arrows. Mbr, Member; Fm, Formation; HMJ, Hamajing Member; SSF, small shelly fossil; Cry., Cryogenian; Cam., Cambrian. See Xiao et al. (2021) for sources of radiometric dates.

The early Cambrian Nagaur Sandstone of the Nagaur Group near Dulmera, Rajasthan, NW India (Fig. 3), contains abundant trace fossils, including Cruziana, Diplichnites, Monomorphichnus, Rusophycus, and Treptichnus, among others (Pandey et al. 2014; Sharma et al. 2018). Ichnofossils are mostly preserved on the soles of sandstone beds, which are characterized by the abundance of trough cross-beds, planar cross-beds, and wave ripples. The fossiliferous sandstones were inferred to have deposited in shoreface environments (Pandey et al. 2014). The maximum depositional age of the Nagaur Sandstone is constrained by the youngest detrital zircons of ~540 Ma (McKenzie et al. 2011), and the unit is variously considered as Cambrian Stage 2 or Stage 4 based on trace fossil assemblages (Pandey et al. 2014; Hughes 2016). A conservative age estimate is thus between 539 and 509 Ma. We selected a 0.039 m² area of Nagaur Sandstone for this study.

Figure 3. Geological map of the Bikaner-Nagaur basin (A) and stratigraphic column of the Marwar Supergroup (B). Inset map in A shows general location of the basin in India. Black star marks fossil locality near Dulmera in NW India. Bb, Boulder bed; Dol, dolostone; Lst, limestone; Gr, Group; Pr, Period; SS, sandstone. Note: stratigraphic column is not to scale, but the Nagaur SS is reported to be approximately 500 m in total thickness, with the underlying Bilara Group approximately 300 m thick (Pandey et al. 2014). After figs. 1 and 2 of Pandey et al. (2014) and courtesy of Dhirenda Pandey.

Materials and Methods

Trace Fossil Markup

The Shibantan bed sole surface was chosen because there was a reasonable area of bedding plane with abundant ichnofossils and limited evidence of preservational or erosional heterogeneity (Fig. 4A; see also Xiao et al. 2021: fig. 6a). Taxonomy of Ediacaran trace fossils can be highly controversial (Darroch et al. 2021), and so here we recognize two distinct morphogroups consisting of relatively large (~1 cm) and small (~1 mm) horizontal traces, with no overlap in width. The smaller ichnofossils are Helminthoidichnites-like, and the larger ones are characterized by poorly preserved spiral structures reminiscent of Streptichnus (Jensen and Runnegar 2005; Xiao et al. 2021; Fig 4E). Importantly, even if this preferred taxonomy is incorrect, the statistical analyses we outline will not be impacted, as these large trace fossils represent similar behaviors, likely made by the same progenitor.

Figure 4. Ediacaran and Cambrian bedding surfaces. A, Shibantan surface analyzed in this study. B, Close-up of black box in A, showing small Helminthoidichnites-type trails (white oval) and Streptichnus with spiral structures (black arrows) indicative of this ichnotaxon. C, Nagaur surface analyzed in this study, with Treptichnus- and Rusophycus-type traces visible. D, A similar surface in the Nagaur Sandstone, illustrating small (white oval), circular (white square), and Treptichnus-type trace fossils. Lines in B and D illustrate the difference between connected (black) and unconnected (white) mark-up strategies of larger horizontal trails. E, Shibantan Streptichnus specimen with the radiating branches. F, Rusophycus from the Nagaur Group.

The Nagaur bed-sole surface similarly contains limited preservational or erosional heterogeneity (Fig. 4C). We mapped four distinct types of trace fossils on this bedding plane. Of primary interest are larger horizontal traces, similar to Treptichnus pedum previously identified from this unit (Pandey et al. 2014), and smaller horizontal trails ~1 mm wide (Fig. 4D, white oval). The larger trails may not contain all the diagnostic features of T. pedum, so we classify them more broadly as Treptichnus isp. These exhibit a probing behavior broadly similar to that of Streptichnus from the Shibantan. Smaller trails were regarded as Planolites isp. (Pandey et al. 2014). Circular traces are of similar diameter to the width of the larger horizontal traces, and the two are often associated. A fourth group of bilobed traces represent Rusophycus and/or Cruziana (Pandey et al. 2014; Fig. 4F).

The outline of each bed and individual trace fossils were marked as vector lines within Inkscape 0.92.3 (Fig. 5). Areas with low fossil density, interpreted as artifacts of differing erosion and/or poor photographic contrast, were excluded (Fig. 5B). Critically, such areas are much larger than the centimeter-scale heterogeneities investigated here, and thus would not affect this study. Most of the traces can be clearly marked; however, two different modes were employed for larger horizontal trails (Streptichnus and Treptichnus). Initially, each individual trace was marked irrespective of adjacent ones (unconnected trails; Fig. 4B,D). Then, in instances where traces appear to represent a series of segments/steps that can be reliably ascribed to the same progenitor, they were connected or simplified to a single segment to focus on the overall paths of the tracemakers (connected trails; Fig. 4B,D). While traces are continuous objects, rather than points, these spatial methods have already been established for modern organisms, such as root systems, so are appropriate to apply to trace segments (Kitchell 1979; Koy and Plotnick 2010). Each single line represents a single action of the tracemaker, so that as long as consistent labeling approaches are applied across each surface, this marking-up approach results in a representation of the behaviors of the tracemakers across the bedding planes.

Figure 5. The mark-up process. A, The Shibantan surface. B, The surface with the area marked up in red, the large ichnotaxon in black, and the small ichnotaxon in green. C, Close-up of area showing the details of how the large and small ichnotaxa were marked up. D, Close-up area showing the midpoints of each trace segment.

Line data were extracted from Inkscape using a custom script written in Haskell (https://github.com/egmitchell/traces). Analyses were performed in R (R Core Team 2017) using spatstat (Baddeley and Turner 2005).

Quantification of Lateral Heterogeneity

Previous paleontological work on foraging behavior has focused on quantifying the lengths and orientations of trace fossils with comparisons to L-systems and Lévy foraging distributions (Plotnick 2003; Sims et al. 2014). Such analyses can be used to investigate foraging behavior by comparing trace segment lengths and orientations to known models. Of these models, one of the best known is Lévy flights, which are thought to optimize foraging strategy (Viswanathan et al. 1996), although the extent to which these accurately describe observed foraging is still very much debated (Edwards et al. 2007; Sims et al. 2008; Humphries et al. 2013; Raichlen et al. 2014). A Lévy distribution predicts that relatively short trace segment lengths represent foraging in areas with plentiful resources, with rare long segments representing the search for nondepleted resource patches (Edwards et al. 2007). This distribution is thought to maximize search efficiency, because the animal(s) will spend more time in high-resource areas, then move quickly far away once the resource is depleted. Lévy distributions are parameterized by a factor μ that quantifies how right-skewed the distribution is, with optimum foraging thought to occur at μ = 2 (Humphries and Sims 2014). Random or minimally optimized strategies can be modeled as Brownian motion and are adequately efficient when resources are abundant and uniformly distributed (Sims et al. 2012). Correlated random walks involve a correlation between successive orientations of each trace segment (Patlak 1953), so produce a direction bias locally, but with diminished directionality over broader spatial scales, thus producing uniform segment orientations (Benhamou 2006). These approaches are used in movement ecology for extant species to analyze long trails covering substantial spatial scales (Humphries and Sims 2014). While such traces are found in the fossil record on surfaces of comparable spatial scales, far more common are relatively small surfaces with multiple partial segments of a longer trail, making such analyses difficult to directly compare (Plotnick 2012). As such, alternative approaches are needed to provide insights into the evolution and development of foraging strategies in the fossil record.

Outside the paleontological record, the level of heterogeneity (LH*) has been applied across a wide range of systems such as plant communities, crime locations, or earthquakes to quantify the degree of patchiness within different systems (Shu et al. 2019). Other methods exist for examining similar spatial patterns (see Shu et al. 2019); however, we prefer LH*, because it provides a single measure for relative heterogeneity over all relevant spatial scales and is easily implemented using spatstat. Spatial patterns are described using point distributions, whereby a homogeneous Poisson process describes randomly distributed points within a mapped area (Illian et al. 2008). The LH* is calculated by summing the difference in nearest neighbor distances between the observed points (here fossil specimens) and the expected random nearest neighbor distribution as described by a homogeneous Poisson model over all spatial scales up to the maximum nearest neighbor distance (Shu et al. 2019). The higher the LH*, the greater the differences between the homogeneous pattern and the observed pattern, so the greater the level of heterogeneity. Therefore, calculation of the LH* for ichnotaxa can be used to evaluate the relative lateral heterogeneity of different tracemakers.

Density plots of the trace surfaces were made to aid visualization of the relative homogeneities, a smoothed density function was used to calculate the relative density across the surfaces rather than plotting individual points/traces. To quantify the level of homogeneity of traces on the Shibantan and Nagaur surfaces, two different analyses were performed. First, the homogeneity test of Hosking and Wallis was used to test whether the spatial distributions were significantly different from what would be expected from a homogeneous distribution (Hosking and Wallis 1993). For the Hosking and Wallis test, the heterogeneity measure approximates a normal distribution with zero mean and standard deviation of 1 and is calculated using the sample coefficients of variation, skewness, and kurtosis. Homogeneity is defined when the heterogeneity measure is less than 1. Second, the LH* was calculated between different ichnotaxa and surfaces. Segment lengths and orientations were extracted for each ichnofossil to determine how they were spatially distributed and enable comparisons to known foraging distributions. Third, the skewness and kurtosis of the distributions of the segment lengths were quantified and tested for normality and lognormality using the Shapiro-Wilk tests and compared with Lévy distributions. The parameter for the best-fit Lévy distribution was found, and the fit of Lévy distributions to normal and lognormal distributions was compared using Akaike information criterion (AIC) values. The proportion of ichnofossils in each orientation cohort and their means were calculated using Mclust (Fraley and Raftery 2017).

Spatial Analyses

The spatial distributions of each ichnotaxon were described by taking the midpoints of each segment, then calculating pair correlation functions (PCFs), which describe spatial patterns. The PCF of a point pattern is calculated as the density of points within a given radius of each point, for all points. Calculating the mean density for a range of different radii (e.g., from 0 to 2 m) allows the changes in point density over the relevant given spatial scale (generally half the minimum width of the mapped area) to be described. The spatial relationship between any pair of ichnotaxa was similarly examined using bivariate PCFs. These determine how the relative density of points (here the ichnofossil midpoints) vary over the mapped area. When PCF = 1 points are distributed randomly within the mapped area, PCF > 1 indicates aggregation or clustering, while PCF < 1 suggests greater spacing (segregation) than expected (Illian et al. 2008).

Similar to PCFs, to establish whether observed patterns are significantly nonrandom, MCF data were compared with a homogeneous Poisson model and 999 Monte Carlo simulations of homogeneous Poisson models, with simulation envelopes defined as the highest and lowest 5% of generated data (Wiegand and Moloney 2013). To further aid comparisons with the random model (homogeneous Poisson model), the goodness-of-fit test p_d was calculated, where p_d = 0 corresponded to no model fit (nonrandom distribution) and p_d = 1 corresponded to a perfect model fit (random). Significant excursions outside the simulation envelope and poor fit under a homogeneous Poisson model (p_d < 0.1) signify nonrandom distributions.

To investigate how trace length and orientation varied over the mapped area, MCFs were used to test for spatial autocorrelation. MCFs are a type of SPPA in which the spatial patterns of a quantitative variable (the mark, here the segment length or orientation angle) for a set of spatial points are tested to see whether the variable exhibits nonrandom spatial behavior, that is, spatial autocorrelation (Grabarnik et al. 2011). When the lengths are the marks for the spatial data, the MCF is calculated as the function (Eq. 1) k _LL(r):

(1)$$k_{LL}( r ) = \displaystyle{{E_{ou}[ {M( O ) M( u ) } ] } \over {E[ {M, \;{\rm \;}{M}^{\prime}} ] }}$$

where E_ou is the conditional expectation, given that there are points at location O with mark M(O) and u with mark M(u) separated by a distance r; E is the expectation; and M and M′ are marks drawn randomly from the distribution of marks. Where MCF = 1, the values associated with each point exhibit no spatial autocorrelation (random distribution). The interpretation of MCF ≠ 1 depends on the function itself, but for all cases indicates spatial autocorrelation (nonrandom behavior). For our analyses, MCF > 1 indicates trace fossils more likely to be aggregated with other specimens with similar marks (positive spatial autocorrelation), and MCF < 1 more likely to be segregated (negative spatial autocorrelation).

To establish whether or not the observed pattern is significantly nonrandom, the observed pattern was compared with a homogeneous Poisson model and 999 Monte Carlo simulations of homogeneous Poisson models were created, with the simulation envelopes defined as the highest and lowest 5% of simulated data (Wiegand and Moloney 2013). To determine the extent to which the observed MCF patterns are described by the random model (homogeneous Poisson model), the goodness-of-fit test p_d was calculated, with p_d = 0 corresponding to no model fit (marks are not randomly distributed) and p_d = 1 corresponding to a perfect random fit, so that the mark correlation function (Eq. 1) exactly describes the relationships between marks and specimen position. Significant excursions outside the simulation envelope where the MCF was not a good fit for the data (p_d >> 0.1) are interpreted as nonrandom.

Results

The Shibantan surface contains a total of 1301 ichnofossils within the sampled area of 0.51 m². Of these, 825 are classified as larger horizontal trails (Streptichnus), reduced to 424 trails when connected. The other 476 trails belong to the small ichnotaxon (Helminthoidichnites). For the Nagaur surface, a total of 1762 ichnofossils were mapped within the sampled area of 0.039 m². These include 1116 large trace fossils (Treptichnus, reduced to 1051 when connected), 476 smaller horizontal trails (Planolites), 189 circular traces, and 25 patches of arthropod scratchings (i.e., Rusophycus and/or Cruziana). All ichnotaxa, apart from Nagaur arthropod traces, exhibited significant lateral heterogeneity, with none exhibiting Lévy distributions (Table 1).

Table 1. Summary metrics for the spatial heterogeneity and segment lengths of analyzed trace fossils. The H-W (Hosking and Wallis) test assesses whether the traces are significantly nonhomogeneous; LH* is a relative measure of heterogeneity; Lévy μ is the parameter for Lévy distributions, where μ = 2.71 is considered an optimal foraging strategy; normality test is the Shapiro-Wilk test to determine whether the segment length distributions are normally distributed.

Lateral Heterogeneity

Density plots (Fig. 6) of the Shibantan surface show notably different patterns between the small and large ichnotaxa, with significant nonhomogeneous patterns (both p < 0.001) and with specimen densities varying up to 20% for the small morphogroup and 7% for the large (Fig. 6A,E, Table 1). The small ichnotaxon exhibited 3.48× higher levels of heteogeneity than the large connected ichnotaxon and 1.58× the levels of heterogeneity for the large unconnected traces (Table 1).

Figure 6. Density plots and spatial distributions of ichnotaxa given by pair correlation functions (PCFs) of trace fossil midpoints. Density maps and corresponding PCF plots are shown in pairs. A, B, Shibantan large ichnotaxon (PCF: solid line, connected; dashed line, unconnected). C, D, Nagaur large ichnotaxon (PCF: solid line, connected; dashed line, unconnected). E, F, Shibantan small ichnotaxon. G, H, Nagaur small ichnotaxon. I, J, Nagaur circular ichnotaxon (PCF of large connected ichnotaxon is plotted for comparison). K, L, Nagaur arthropod traces. M, N, Shibantan bivariate distribution between large connected and small ichnotaxa. O, P, Nagaur bivariate distribution between large connected and small ichnotaxa. Scale bar next to density map is percent coverage, with range depending on the ichontaxon. PCF plots: blue lines, Shibantan; purple, Nagaur; gray area, 999 Monte Carlo simulations with a 5% threshold. PCF > 1 is aggregation; PCF < 1 is segregation; and PCF = 1 is random.

Density plots of the simple horizontal trails on the Nagaur surface (Fig. 6C,G) show notably different patterns between the small and large ichnotaxa, with significant nonhomogeneous patterns (both p < 0.001) and with varying specimen densities up to 14% for the small morphogroup and 20% for the large (Fig. 6C,G, Table 1). The small ichnotaxon exhbited 3.97× higher levels of heteogeneity than the large unconnected ichnotaxon and 2.39× the levels of heterogeneity for the connected large traces (Table 1). The spatial heterogeneity of the circular traces was significantly nonhomogeneous (p < 0.001), with varying specimen densities up to 0.14% (Fig. 6I, Table 1). The arthopod traces were not significantly heterogeneous (p = 0.219), which could be in part due to the relatively small number of these trace fossils (Fig. 6K).

Distribution of Trace Lengths

For the Shibantan, the small ichnotaxon was more strongly skewed in terms of segment length, but the large connected ichnotaxon had a greater kurtosis (Fig. 7A–C, Table 1). Neither ichnotaxon was normally distributed (p << 0.1; Table 1). The smallest ichnotaxon exhibited a lognormal distribution (p << 0.001; Table 1) but the large connected ichnotaxon did not (p >> 0.1; Table 1). For the Nagaur, the large ichnotaxon was more strongly skewed with greater kurtosis (Fig. 7G–I, Table 1). None of the ichnotaxa were normally distributed (p << 0.1; Table 1), and the large ichnotaxon was not lognormally distributed (p << 0.1; Table 1).

Figure 7. Shibantan data are presented in the top two rows (A–F), Nagaur data in the bottom two rows (G–L). First (A–C) and third (G–I) rows are histograms of segment lengths, and second (D–F) and fourth (J–L) rows are rose diagrams of segment orientations. The three columns are for small ichnotaxon (A, D, G, J), large (unconnected) ichnotaxon (B, E, H, K), and large connected ichnotaxon (C, F, I, L), respectively.

The best-fit Lévy distributions for large and small horizontal trails on both surfaces were not close to the optimum of 2 (Sims et al. 2019; Table 1). AIC showed a worse fit for Lévy distributions than normal distributions, and modeled Lévy distributions had much higher kurtosis and skewness. Thus, these distributions were not well described by Lévy distributions, so such comparisons have limited value. Both Brownian walks and correlated random walks produce uniform distributions with skewness and kurtosis much lower than the observed distribution, so these models are unlikely to describe any of the observed data.

Spatial Analyses

For the Shibantan, the midpoints of both ichnotaxa exhibit significantly nonrandom spatial distributions (Fig. 6B,F). The distributions of large connected and unconnected trails are the same at spatial scales greater than 1.5 cm, with unconnected trails showing higher aggregation (PCF_max = 2.0) than connected trails (PCF_max = 1.5; Fig. 6B). The large ichnotaxon has a single scale of aggregation at 1–8 cm of ~1.5× the random density (p_d = 0.002; Fig. 6B). Small traces are significantly aggregated (p_d = 0.002) with ~3× the random density under 1 cm and 1.5× the random density under 8 cm (Fig. 6E,F). Bivariate midpoint PCF shows a single scale of aggregation at 1–8 cm of ~1.5× the random density (p_d = 0.025; Fig. 6M,N).

For the Nagaur, the midpoints of small, large, and circular ichnotaxa exhibited significantly nonrandom spatial distributions (Fig. 6D,H,J). Spatial distributions for the large connected and unconnected trails are the same, >0.5 cm, showing similar levels of aggregation (PCF_max ~ 2.0), but with the unconnected trails significantly segregated below 2 mm (p_seg = 0.001; Fig. 6D). Segregation is most likely due to how connected versus unconnected trails were identified—when trails are less than 2 mm apart, it is difficult to distinguish between them, thus, they were likely counted as a single trail. The large ichnotaxon is aggregated at two spatial scales, under 0.5 cm and between 0.5 and 2.3 cm of ~1.5× the random density (p_d = 0.001; Fig. 6C). Small traces are significantly aggregated (p_d = 0.001) with ~8× the random density under 0.5 cm and 2× the random density under 1 cm (Fig. 6G,H). The circular traces showed a significantly nonrandom distribution (p_d = 0.001), closely matching the aggregation of the large ichnotaxon, 0.5 cm (Fig. 6I,J), and consistent with the hypothesis that these were produced by the same organism but record slightly different behaviors. Arthropod traces were randomly distributed (p_d = 0.267; Fig. 6K,L). The bivariate midpoint PCF exhibits significant segregation between small and large trails across most spatial scales, with a single scale of aggregation between 2.75 and 3 cm of ~1.5× the random density (p_d = 0.001, p_seg = 0.001; Fig. 6O,P).

Small and large ichnotaxa on both surfaces exhibit significant autocorrelation of segment lengths. The Shibantan ichnotaxa are aggregated at similar segment lengths under 2.0 cm (small p_d = 0.011, large p_d = 0.001; Fig. 8A,C), with large unconnected traces showing a stronger correlation at smaller distances of <1 cm (MCF_min = 0.6 vs. 0.8; p_d = 0.001; Fig. 8A). The Nagaur ichnotaxa are aggregated at similar segment lengths under 0.5 cm (small p_d = 0.001, large p_d = 0.016; Fig. 8B,D), with large unconnected traces showing a slightly stronger correlation (MCF_min = 0.58 vs. 0.62; p_d = 0.016; Fig. 8B).

Figure 8. Mark correlation functions (MCFs) of spatial patterns of the lengths (A–D) and orientations (E–H) of traces. MCF < 1 is aggregation; MCF > 1 is segregation; and MCF = 1 is random. Blue lines are for the Shibantan and purple for Nagaur. For the large ichnotaxon, solid lines are the connected traces and dashed lines are unconnected.

In the Shibantan, the orientations of the large ichnotaxon did not show significant spatial autocorrelation (p_d = 0.997), although the large unconnected traces are borderline significant (p_d = 0.105; Fig. 8E). The small ichnotaxon did not exhibit any significant autocorrelation (p_d = 0.770; Fig. 8G). The lack of evidence for autocorrelation contrasts Nagaur ichnotaxa, which all show significant orientation autocorrelations (Fig. 8F,H). The Nagaur large ichnotaxon is significantly aggregated <1 cm (p_d = 0.012) with a weaker signal from the unconnected traces (p_d = 0.086). The Nagaur small ichnotaxon also exhibited significant autocorrelation (p_d = 0.082) with spatial aggregation of orientations <1 cm and significant segregation between 1.3 and 2.2 cm (Fig. 8H).

Discussion

Potential Limitations

The main focus of this study is to provide a proof of concept that SPPA can be used to reliably identify patterns in trace fossils exposed on horizontal bedding planes. We chose to analyze two terminal Ediacaran–early Cambrian slabs that had been previously documented by Z.C. and S.X. A clear limitation is that the two beds may not be representative of each geological period, with a variety of potential taphonomic factors contributing to the density and distribution of trace fossils on each surface, rather than the hypothesized evolutionary shifts in animal behaviors highlighted below. The two surfaces, although both from deposits with evidence for regular current activity in a shallow-marine environment, are interpreted to represent different depositional settings, given their different lithologies (carbonate vs. siliciclastic rocks). Therefore, paleoenvironmental factors may contribute to the different patterns identified. For example, animals can only leave horizontal burrows while a bed is exposed or very shallowly buried, so time of exposure will exert a significant control on the number of trace fossils on any given surface. If sufficient time is allowed and resources are plentiful, burrows will begin to overprint one another, leading to time averaging and altering our ability to resolve the full extent of past behavior. Importantly, these taphonomic uncertainties do not negate the utility of this method. Application to a broader suite of horizontal bedding planes with abundant trace fossils across the Ediacaran/Cambrian boundary will provide further tests of the hypotheses presented here.

Trace fossils provide essential insight into the behaviors of ancient organisms; however, there are also many limitations associated with the ichnofossil record. Of particular importance here is a lack of understanding of exactly what behaviors are being recorded by horizontal trails. Even with confident assignment to Streptichnus, it is unclear whether individual segments represent systematic probing (the preferred interpretation), similar to Treptichnus pedum, or are instead related to the movement of the tracemaker (Hughes 2016; Darroch et al. 2021). Although likely related to the presence of an organic mat, smaller, indistinct trails may represent movement and/or feeding (Jensen et al. 2006). Such ambiguities also leave open the question of exactly what resources (e.g., oxygen vs. food) were being exploited.

Analyses of trace segment length and orientation have been a key focus in extant movement ecology (e.g., Houghton et al. 2006; Humphries and Sims 2014; Kölzsch et al. 2015). In this study, we did not have long single traces or time dimensions, but instead had multiple shorter traces that limit the applicability of methods developed in extant movement ecology (for an overview, see Edelhoff et al. 2016) and further motivate the need to explore alternative approaches, such as the ones presented here. Comparisons with extant movement ecology could be made if the partial trails included in our dataset represent a random subsample of the trails of a single tracemaker population. Theoretical model comparison would be needed to test this hypothesis and to assess how robust such subsamples are, but preservational effects, such as possible biases toward the preservation of shorter trails, complicate theoretical modeling. An inflated abundance of shorter trails would increase the tendency of segment length distributions toward Lévy distributions, which are highly skewed. If the validity of such approaches could be established or much longer traces could be used, the comparison of segment length distributions to the best-fit Lévy model would enable the comparison of the relative extent of optimization for ichnotaxa.

Animal Behavior across the Ediacaran/Cambrian Boundary

Despite the abovementioned limitations, all ichnotaxa investigated represent energetically demanding behaviors, likely with some benefit gained by their progenitor. The nonrandom distribution of most trace fossils identified here suggests that the tracemakers responded to heterogeneously distributed resources, regardless of what resources were targeted. We note that both redox conditions and microbial mat distribution in Ediacaran–Cambrian shallow-marine environments were spatially heterogenous and temporally dynamic, as evidenced by redox proxy data (e.g., Wood et al. 2015) and microbially induced sedimentary structures as well as textured organic surfaces (e.g., Buatois et al. 2014; Droser et al. 2019; Darroch et al. 2021). Thus, resource heterogeneity is expected. Perhaps most importantly, the increased sophistication and specialization observed here would result in greater fitness of the trace-making animals in environments with heterogeneous resource distribution.

In ichnology, previous metrics for describing the extent of lateral bioturbation (Miller and Smail 1997; Marenco and Hagadorn 2019) have limitations in detecting patchiness at smaller spatial scales (Marenco and Hagadorn 2019), cannot distinguish between significantly nonhomogeneous patchiness, and thus are unable to quantify the spatial scale of any such patchiness. Spatial patterns are often hard to identify with the naked eye, especially when assessing the significance of a pattern (Law et al. 2009). Our use of the Hosking and Wallis homogeneity test enabled us to statistically assess whether the heterogeneity of each ichnotaxon was significantly different from a homogeneous pattern. The LH* metric indicates that the trails of the small ichnotaxon had greater heterogeneity than those of the large ichnotaxon on both surfaces. These two metrics provide a method to compare the extent of patchiness across beds with any level of coverage and to identify significant heterogeneities at much lower ichnofossil abundance. The ability to detect heterogeneities across a wide range of percent coverage indicates that these methods are applicable throughout the fossil record. Analysis of the Cambrian surface indicates distinct spatial differences between the arthropod scratches and other traces. This result is expected, given that these arthropod traces are interpreted to represent distinct behaviors (namely not foraging), especially because circular traces were likely related to the larger treptichniids and thus have a very similar spatial distribution to the treptichniids (Fig. 6J). Understanding how these metrics change beyond the early Cambrian with the increased depth and complexity of bioturbation observed through the Phanerozoic will also provide meaningful context for the patterns described here.

Aggregated spatial distributions of trace fossil midpoints from the Shibantan suggest preferential focus on optimal habitats. The smaller morphogroup showed more complex spatial patterns, with a strong (2.8×) aggregation under 2 cm and a weaker (1.5×) aggregation at larger spatial scales (Fig. 6B). These differences suggest two distinct external factors responsible for such aggregation. Similar scale and strength of aggregation between midpoint spatial distribution of the large and small ichnotaxa as well as their bivariate distribution suggest common resource exploitation. Some optimization was detected through the MCF analyses of segment lengths, which indicates decreased segment length in quality patches to optimize time spent there, but a lack of ability to turn back once they left the preferred area. Taken together, these analyses indicate the Shibantan large ichnotaxon represents relatively simple foraging behavior, likely responding to heterogeneities of a single resource with limited evidence for optimization via avoidance or return to optimal resource patches (cf. Fig. 1B).

The Nagaur trace fossils exhibited greater spatial aggregation than those from the Shibantan, with the small ichnotaxon showing significantly stronger (8×) aggregation than the large ichnotaxon (2×). Both Nagaur ichnotaxa have two scales of spatial aggregation, a stronger one under 0.5 cm and a weaker one at 0.5–1.5 cm for the small ichnotaxon and 0.5–2.5 cm for the large ichnotaxon. Differences in the PCF distributions of the small and large ichnotaxa suggest that they targeted different resources, further exemplified by strong segregation of the Nagaur bivariate PCF distribution (Fig. 6P). While some small-scale segregation may be preservational (e.g., large tracemakers erasing small traces), this process cannot explain the segregation over the majority of spatial scales. Thus, spatial distributions suggest that the Nagaur tracemakers had a greater degree of niche specialization than their Shibantan counterparts.

Strong directionality of small and large trace fossil segments on the Shibantan and Nagaur beds (Fig. 7, Table 2) runs against many models of movement ecology. Lévy flights, Brownian motion, and correlated random walks all describe more uniform distributions (Auger-Méthé et al. 2015). Strong directionality in extant taxa can arise when organisms reverse direction after straying out of food sources (Kölzsch et al. 2015). These changes in directionality could be driven by a variety of signals, including those from nonlocal sources (Fagan et al. 2017). Here, we would expect the orientations to have a spatial distribution that corresponds to the spatial distribution of food sources, because nonlocal sources would likely be operating over larger spatial scales, not providing the heterogeneity seen at the spatial scales we have described. By testing for the presence of spatial autocorrelation of orientations, and through comparisons of the spatial scales of any significant autocorrelation, we can further establish whether tracemakers do reverse direction after straying out of a high-quality resource patch.

Table 2. Bimodal orientation distributions of the trace orientations.

In the Shibantan, the orientations of the trace fossils do not exhibit significant spatial autocorrelation (Fig. 8E,G). Nonrandom spatial distributions of midpoints suggest preferential directionality was not focused on the sources of midpoint aggregations. As such, it is unlikely that the strong directionality exhibited in Ediacaran trace fossils analyzed here is due to resource distribution. More likely, the directionality is related to currents, as observed, for example, in fossil horseshoe crab trails (Buhler and Grey 2017). These patterns in the Shibantan show a second cohort of orientations approximately 160° to the majority of trails. Future studies targeting similar trace fossils associated with paleocurrent indicators (e.g., ripple marks) could test this hypothesized relationship.

Sophistication in the Nagaur ichnotaxa is captured in the MCF segment length analyses (Fig. 8B–D). Like the Shibantan trace fossils, MCFs indicate that small segment lengths likely represent focused behavior in beneficial resource patches (Fig. 8A–D). These Cambrian trace progenitors appear to have been able to detect optimal areas, only making minimal changes in direction to maximize time within the patches. In sharp contrast to the Shibantan, significant autocorrelation of trace orientations in the Nagaur (Fig. 8F,H) is consistent with tracemakers changing direction dependent on the spatial distribution of resources. This spatial autocorrelation of the Nagaur small ichnotaxon suggests that it was able to detect when it was leaving the resource patches and to make a significant change in orientation to enable a return to those patches (cf. Grabarnik et al. 2011). Small spatial–scale aggregations suggest the maintenance of similar paths in preferred resource patches. Significant segregation of the small ichnotaxon at medium spatial scales further suggests that the small tracemaker changed direction consistently.

The evolution of optimization and the capability to search for patchy resources may be crucial to clarifying early bilaterian ecology, as the motivation to find these resources may be a driver of early animal adaptations (Koy and Plotnick 2007, 2010; Plotnick et al. 2010; Budd and Jensen 2017). The ubiquity of Precambrian organic mats likely declined in the Phanerozoic with the advent and diversification of vertical bioturbating animals, providing new pressures for optimized foraging in the Cambrian that are absent or entirely different in the Ediacaran. Resulting adaptations could potentially start a feedback of increasing diversification due to increased heterogeneity through increased foraging (Mitchell et al. 2020). Despite evidence for heterogeneity in the Shibantan, the two ichnotaxa studied demonstrate a relatively limited capacity for resource-focused foraging in the Shibantan relative to the Nagaur. Further, the Nagaur ichnotaxa showed more behavioral complexity, with different niches occupied by horizontal grazers and intertaxon avoidance. Due to the limited sample sizes, it is not clear whether these samples are representative or to what extent taphonomy and/or paleoenvironment (among other factors) influenced our results, necessitating further work on surfaces across the Ediacaran/Cambrian transition to test potential evolutionary patterns. However, the fact that these patterns match previous observations of an increase in behavioral complexity across this boundary suggests that the methods employed here revealed a real evolutionary signal.

Further ground truthing based on analyses of younger deposits will continue to establish the utility of the methods described in this paper. Extending this work to a broader range of ichnotaxa and more bedding-plane surfaces will continue to test the hypotheses suggested here and identify the range of foraging efficiencies exhibited by early animals. However, increased ichnodisparity suggests that there was pressure to optimize foraging by the terminal Ediacaran, potentially contributing to changes in diversity at the Ediacaran/Cambrian transition. Only through analyses of trace fossil assemblages across this interval will we be able to detect and quantify how the development of complex trace morphologies relates to optimization of foraging. The ability to recognize differences between these Ediacaran and Cambrian behaviors suggests that spatial analysis is a powerful method for interrogating the trace fossil record.

Conclusions

Analysis of two trace fossil–rich surfaces from the Ediacaran Shibantan Member and the Cambrian Nagaur Sandstone shows that ichnotaxa on both surfaces exhibit lateral heterogeneity, likely in response to unevenly distributed resources. In the Ediacaran, both ichnotaxa investigated appear to have exploited the same resource. The small ichnotaxon also showed added aggregation at smaller spatial scales, indicating the preference for a second, distinct resource. Bivariate spatial patterns suggest that these Shibantan tracemakers likely did not interact directly with each other but shared similar resource requirements. Despite these patterns, foraging strategies were comparably limited, with only the smaller Shibantan ichnotaxon seemingly able to maximize time in presumed high-quality resource patches and neither expressing the ability to change orientation to concentrate on particular resources. In contrast, early Nagaur trace fossils exhibit notable increases in optimization of foraging behavior, specifically the ability of Nagaur tracemakers to adjust segment length and orientation to focus on preferred resources.

This study provides preliminary but quantitative data that Shibantan tracemakers exhibited limited evidence of resource optimization and that their Nagaur counterparts developed more sophisticated foraging strategies. Further work is needed to establish the extent to which these patterns reflect environmental difference between the two analyzed ichnofossil assemblages or behavior-evolutionary innovation across the Ediacaran/Cambrian transition. More broadly, the analytical techniques employed in this study provide a novel methodology for quantifying the lateral heterogeneity of bioturbation, the methods to test for spatially induced turning, and approaches to assess how tracemakers interact with lateral resource heterogeneity in their local environments.