Introduction

Drilling predation on shelled marine invertebrates represents a unique opportunity to quantify biological interactions in ancient and modern marine environments. Predatory drillholes represent direct records of predatory attack and can potentially document diverse aspects of predator–prey interactions (e.g., Kitchell et al. 1981; Kitchell 1986; Dietl and Alexander 2000; Kowalewski 2002; Hoffmeister et al. 2004; Rojas et al. 2017), although several assumptions and caveats need to be considered (e.g., Kowalewski 2002; Leighton 2002; Harper 2006; Klompmaker et al. 2019; Smith et al. 2019). Drillholes can provide valuable behavioral information regarding selectivity of predatory attacks in terms of prey species, prey size class, or drilling location on the prey skeleton (e.g., Kitchell 1986; Kowalewski 2002). Here, we focus on site selectivity and illustrate a reproducible approach that can be applied to quantify spatial patterns in the distribution of predation traces on modern and fossil prey remains.

Modern studies (Dietl and Alexander 2000; Chiba and Sato 2012) as well as investigations based on the fossil record (Kelley 1988; Dietl et al. 2001) have shown that drilling predators can display high levels of spatial stereotypy, likely reflecting prey handling during the attack and/or prey morphology (Kabat 1990; Kingsley-Smith et al. 2003; Rojas et al. 2015). Whereas many studies have documented the spatial distribution of drillholes for various prey groups, approaches used to assess those spatial patterns have varied widely, including qualitative descriptions (Negus 1975) and a multitude of quantitative approaches (Kelley 1988; Kowalewski 1990; Anderson et al. 1991; Dietl and Alexander 2000; Hoffmeister and Kowalewski 2001). These methodological inconsistencies hamper our ability to compare outcomes across studies or conduct meta-analyses based on previously collected data. In addition, most of the methods for analyzing site selectivity in drilling predation were primarily developed to test the null hypothesis of a random distribution of traces across arbitrarily defined sectors of the prey skeleton using goodness-of-fit, chi-square or Kolmogorov-Smirnov statistics (Kelley 1988; Kowalewski 1990; Anderson et al. 1991), and diversity metrics such as the Shannon-Weaver index (Dietl et al. 2001). These approaches fail to exploit the high-resolution information locked up in the spatial relationship between drillhole locations and depend critically on an arbitrary grid system used to divide the prey skeleton into sectors (Johnson 1984; Kowalewski 2004). Recently, the nearest neighbor index, which requires no artificial partitioning scheme, was introduced by Alexander et al. (2007) and implemented by Casey et al. (2015) to test the null hypothesis that predation traces on the prey skeleton occur in a completely random fashion (complete spatial randomness [CSR] hypothesis). However, this approach was not designed to identify patterns that may occur at different spatial scales or distances measured on the prey skeleton, for example, aggregation at short distances and segregation at large distances.

The recognition that predation traces on prey skeletons are spatially explicit and can be mapped enables the development of a proxy for site selectivity in drilling predation based on spatial point process modeling (Baddeley et al. 2016). Here, we introduce the spatial point pattern analysis of traces (SPPAT), an approach for visualizing and quantifying the distribution of predation traces on shelled invertebrate prey, which includes improved collection of spatial information inherent to drillhole location, improved visualization of spatial trends, and distance-based statistics for hypothesis testing (see also Clapham et al. [2003], Mitchell and Butterfield [2018], and Mitchell et al. [2019] for examples that use a similar spatially explicit approach to describe the distribution of well-preserved Ediacaran fossils). We illustrate the SPPAT approach through case studies on museum samples of the Plio-Pleistocene venerid bivalve Lirophora latilirata (Conrad, 1841) from the Atlantic coastal plain of the United States, which has been previously used in studies on drilling predation (Hattori et al. 2014; Klompmaker and Kelley 2015); drilling data from laboratory-based feeding trials of the tropical eastern Pacific naticid Notocochlis unifasciata (Lamarck, 1822) preying upon the venerid bivalve Iliochione subrugosa (W. Wood, 1828); and modern beach-collected samples of I. subrugosa from Central America. The SPPAT approach provides tools for visualizing distribution patterns of predation traces on shelled invertebrate prey that could not have been retrieved from current methods and uses distance-based functions to quantify those patterns.

Material and Methods

Case Studies.—We used fossil samples, modern beach-collected samples, and laboratory feeding trials to illustrate the utility of the SPPAT approach in assessing site selectivity of drilling predators. First, fossil specimens of the abundant, commonly drilled venerid bivalve L. latilirata from the Plio-Pleistocene of the Atlantic coastal plain of the United States were obtained for site-selectivity analysis by conducting a survey of the Invertebrate Paleontology Collections in the Florida Museum of Natural History (FLMNH), University of Florida (specimen acronym UF). The museum survey was focused on bulk-collected specimen lots (41 lots containing 1376 specimens) that were associated with large collections (Table 1). Drilled specimens (n = 166; Table 1) of L. latilirata were identified by examining shells under a binocular microscope for the presence of complete naticid-like drillholes, which are typically beveled in shape, and categorized as the ichnospecies Sedilichnus paraboloides (Bromley 1981) (formerly Oichnus paraboloides). We used the L. latilirata case study to illustrate how the SPPAT approach can be used to quantify large-scale, temporal shifts in site-selective behavior of drilling predators. The stratigraphic details of the fossil samples used in this study can be found in Supplementary Table S1.

Table 1.

Drilling data on museum samples of Lirophora latilirata and Iliochione subrugosa compiled in this study. *Stratigraphic context follows Lyons (1991) and Ward et al. (1991). Abbreviations: DF, drilling frequency.

Second, we reexamined data from Rojas et al. (2015) for 27 drilled specimens of the venerid bivalve I. subrugosa that were preyed upon by the naticid N. unifasciata in a laboratory setting and 89 drilled specimens of I. subrugosa that were target collected (Ottens et al. 2012) at two modern beach localities (Veracruz and El Palmar) in Panama, Central America (Table 1). We used the I. subrugosa case studies, that is, modern beach-collected samples and laboratory-based feeding trials in which numerous cases of edge-drilling predation were identified, to illustrate how the SPPAT approach can be used to quantify small-scale (intravalve) spatial variation in the relative importance of alternative shell-drilling behaviors (i.e., edge and wall drilling).

Quantifying Drilling Location.—Drilled valves were digitally photographed using a standard photogrammetric protocol adapted from Perea et al. (2008). All images were oriented such that the anteroposterior axis was horizontal to facilitate consistency in landmark data collection (Kolbe et al. 2011). The positions of the predatory drillholes were quantified using a two-dimensional morphometric approach proposed by Roopnarine and Beussink (1999) and adapted by Rojas et al. (2015) to assemble the collected drilling data into a standardized prey skeleton. The general procedure was based on the following pseudolandmarks identified on the external view of the valves and placed on the same plane (Supplementary Fig. S1): (1) point of maximum curvature of the ventral edge, (2) anterior end of the valve, (3) the beak on the outline, and (4) posterior end of the valve. A fifth point corresponds to the center of the predatory drillhole. A line joining pseudolandmarks 2 and 4, located on the anterior and posterior margins of the bivalve shell, approximates maximum length of a prey specimen. Casey et al. (2015) adopted a similar approach to quantify drilling locations on the venerid bivalve Mercenaria. These pseudolandmarks do not correspond with specifically internal homologous traits of the prey, as in Roopnarine and Beussink (1999). Instead, they reflect the general external morphology of the prey skeleton, which is an important consideration for drilling predators when manipulating their prey (Kabat 1990; Kingsley-Smith et al. 2003; Rojas et al. 2015). Landmark coordinates obtained for left valves were inverted to compare directly with the right valves following Kowalewski (2004). The Bookstein baseline registration method for two-dimensional data (Bookstein 1986) was used to remove shell size, position information, and remaining orientation from the analyses using pseudolandmarks 1 and 3 as a baseline. This dataset, given the measured spatial locations of drillholes on the standardized prey skeleton (i.e., study area in Bookstein shape units), represents a spatial point pattern (Bivand et al. 2008) (Fig. 1). See Supplementary Table S1 for point-level attributes. The R functions for landmark-based morphometrics of Claude (2008) were used to carry out the morphometric analysis.

Figure1.

Spatial point pattern derived from pooled drilling data on Plio-Pleistocene specimens of the bivalve Lirophora latilirata from the Atlantic coastal plain.

Creating Spatial Point Patterns of Drilling Locations.—A key feature of the SPPAT approach is the creation of a point pattern of drilling predation traces derived via geometric morphometrics. To accomplish this task, all landmark data on drillhole locations for L. latilirata and I. subrugosa were separately pooled and superimposed on a standardized skeleton to create a point pattern, that is, a representation of a given area containing a set of locations (hereafter referred to as “L. latilirata samples” and “beach-collected samples,” respectively; Figs. 2, 3). The pooled point pattern of predation traces on L. latilirata (Fig. 1) was further grouped into sub-patterns (sensu Baddeley et al. 2016) according to geologic age, which are also referred here as point patterns and used to illustrate temporal changes in site-selectivity patterns, including variations in alternative shell-drilling behaviors.

Kernel Density Estimation (KDE) and Hotspot Mapping.—KDE was used to visualize the spatial distribution of predation traces on the prey skeleton for the modern beach-collected and laboratory feeding trial samples of I. subrugosa. This technique fits a curved surface over each drilling location such that the surface is highest above the drillhole and zero at a specified distance (bandwidth, h) from the drillhole. Following Silverman (1986: p. 76, Eq. 4.5), the kernel density was mathematically expressed as:

where ƒ(x) is the density value at drillhole location x, n is the number of drillholes, h is the bandwidth, x – X_i is the distance to each drillhole location i, and K is the quadratic kernel function, which is defined as

Figure2.

Point patterns, kernel density, and hotspot maps of drillholes on fossil samples of Lirophora latilirata from the Atlantic coastal plain grouped by time interval. A, Late Pliocene. B, Early Pleistocene. C, Middle Pleistocene. Kernel density estimated using drilling frequency (DF) as a weighting variable. Units of density maps are number of drillholes per area measured in square Bookstein shape units. Hotspots were detected from the kernel density map using as a threshold the highest 10% estimated values.

Because fossil samples of L. latilirata, unlike our modern beach-collected and laboratory feeding trial samples of I. subrugosa, came from several field surveys, and museum lots examined varied greatly in size, KDE estimates may be biased. Therefore, we implemented a correction for sampling effort using drilling frequency. Drilling frequency (DF) was calculated for each fossil lot l following Kowalewski (2002: p. 19, Eq. 5):

where n_l is the number of valves drilled in the museum lot l, and m_l is the total number of valves in the same lot l. DF was used to reweight the KDE estimated at each drillhole location x and thus to control for the role of different museum lots in the resulting site-selectivity patterns. Drillholes from museum lots with higher DF were assigned greater weight and receive higher KDE estimates. This standard procedure for reweighting the contributions from the individual predatory traces can be used for density estimation when observations are sampled nonuniformly (Fieberg 2007). In our study, for instance, museum lots with large numbers of specimens and an unusually high occurrence of edge-drilling may bias density patterns when pooling drilling data.

Figure3.

Point patterns, kernel density, and hotspot maps derived from drilling data on beach-collected samples of Iliochione subrugosa (Playa Veracruz and Playa El Palmar, eastern Pacific coast of Panama) and feeding trials of Notocochlis unifasciata preying upon I. subrugosa. A, Beach samples. B, Feeding trial samples. Units of density maps are number of drillholes per area measured in square Bookstein shape units. Hotspots were detected from the kernel density map using as a threshold the highest 10% estimated values.

The optimum bandwidth h_opt required for KDE was based on the number of predatory traces and the average distance from each trace to the mean center (e.g., standard distance) (Fotheringham et al. 2000). Descriptive statistics and parameters employed to calculate the optimum bandwidth h_opt used for the KDE are presented in Table 2. Because drillhole configurations are scaled and aligned on the baseline of coordinates (–0.5, 0) and (0.5, 0) derived from Bookstein registration, measurement units (e.g., millimeters) cannot be provided for the density maps. KDE was performed using the function density.ppp (grid spacing eps = 0.01 map units; h = h_opt; weights = DF) in the R package Spatstat v. 1.23-4 (Baddeley and Turner 2005). Based on the estimated density, preferences in drilling location by predators can be reliably characterized as hotspots of drilling traces on the prey skeleton. Operationally, hotspot areas were detected from the kernel estimate map using the highest 10% of kernel estimated values as a threshold (Nelson and Boots 2008).

Distance-based Statistical Analyses.—The spatial distribution of predation traces on a given prey skeleton may be classified as random (CSR), aggregated (individual traces are placed together delineating groups or clusters), or segregated (individual traces are placed farther apart than they would be in CSR) or some combination of these patterns at different distances on the prey skeleton. Here we describe the K- and L-functions, which are used for hypothesis testing, and applied the maximum clustering distance measurement to indicate the distance on the prey skeleton at which traces were clustered.

Table 2.

Descriptive statistics used to calculate the optimum bandwidth (h_opt).

We used Ripley's K-function (Ripley 1976, 1977) to quantify the spatial point patterns of predation traces assembled here. This summary statistic is widely used to quantify deviations from CSR in a statistically consistent framework (Fotheringham et al. 2000; Baddeley et al. 2016). Following Ripley (1988), the estimate of K(r) employed in this study was mathematically expressed as:

where a is the area of the standardized prey skeleton, n is the number of drillholes, and the sum is taken over all ordered pairs of drillholes i and j in the point pattern. d_ij is the distance between the two drillholes and I (d_ij ≤ r) is the indicator that equals 1 if the distance is less than or equal to r. The term e_ij is a standard correction factor to minimize potential bias from edge effects in estimating the K-function (Baddeley et al. 2016). The expected value of K(r) for a random Poisson distribution (K_pois) is πr² and deviations from this expectation indicate either an aggregated or segregated distribution. This function provides the average number of neighboring drillholes that lie within a distance r of a drillhole divided by the density of drillholes per unit area or intensity. Because it is adjusted for intensity, the K-function makes it possible to compare the degree of regularity between point patterns of drillhole locations with different average densities (Baddeley et al. 2016), including both temporal and spatial comparisons. To facilitate the visual assessment of the graphic output, we used Besag's (1977) L-function, which transforms the K-function K_pois(r) = πr2 to the straight-line L_pois (r) = r:

The square-root transformation stabilizes the variance of the empirical L(r), making it easier to assess deviations from CSR and thus to evaluate the following scenarios: (1) null hypothesis [CSR]: L(r) = r; (2) alternative hypothesis 1 [aggregated]: L(r) > r; and (3) alternative hypothesis 2 [segregated]: L(r) < r. The nonparametric test for CSR in drillhole locations was constructed by simulating 999 random spatial patterns of size n from a homogeneous spatial Poisson process and taking the 50^th largest among the simulated values (L_theo) to achieve a desired significance level α = 50/ (999 + 1). The empirical (L_obs) and theoretical (L_theo) L-functions were estimated for a range of r values measured in Bookstein shape coordinates using the Lest function, and critical values were computed using the envelope function. The point of maximum difference between the observed and expected L values (L_obs – L_theo), termed maximum clustering distance (MCD), was used to indicate the distance at which predation traces were aggregated (Mullins et al. 2015). Because the K-function is cumulative and the graphic tool used for statistical inference may be difficult to interpret (Baddeley et al. 2016), we contrasted our results from the L-function against the pair correlation function (PCF). The PCF describes how the (normalized) density of predation traces changes as a function of distance, with the value PCF = 1 indicating CSR, PCF > 1 aggregation, and PCF < 1 segregation (Illian et al. 2008). We contrasted our results against the nearest neighbor index (NNI) (Clark and Evans 1954), a procedure that was previously used by Casey et al. (2015) to quantify deviations from CSR in the location of drilling predation traces. NNI is defined as the ratio of the mean nearest neighbor distance observed divided by the average distance between neighbors under CSR, with NNI < 1 indicating aggregation, and NNI > 1 segregation of traces. We also used a general approach that converts the results calculated at different distances along the prey skeleton into a single summary statistic u, the goodness-of-fit test (GoF) described by Loosmore and Ford (2006). The functions used for spatial statistical analyses are available in the R package Spatstat v. 1.23-4 (Baddeley et al. 2016).

Results

Description and Visualization of Spatial Patterns.—Overall, the spatial point pattern of predation traces for L. latilirata samples showed that drillholes were spread around the shell wall, ranging from the anterior margin to the umbonal area, with some traces located directly on the shell margin near the ventral anterior and posterior regions (Fig. 1). When point patterns for fossil L. latilirata were grouped by geologic age and examined through time, areas of higher concentration of predation traces derived via kernel density were irregularly distributed on the prey skeleton (Fig. 2A–C). Hotspot maps indicate a subtle migration of a single hotspot of predation traces on the prey skeleton of L. latilirata through time, from the dorsal wall region in the late Pliocene (Fig. 2A) toward the central region of the wall area in the early Pleistocene (Fig. 2B) and the anterior region in the middle Pleistocene (Fig. 2C). Only a few cases of edge drilling, a mode of predation in which a predator drills a hole at a point on the commissure between the closed valves (Taylor 1980; Vermeij 1980; Ansell and Morton 1985; Mondal et al. 2014; Rojas et al. 2015), were observed in the fossil L. latilirata samples, and thus no hotspots of predation traces were evident (Fig. 2).

In the point pattern of predation traces derived from beach-collected samples of I. subrugosa, two areas of higher concentration of drilling traces were found, one located at the umbo and the other a point placed at the ventral commissure, reflecting the occurrence of wall- and edge-drilling behavior of the unknown naticid predator(s). However, only a single hotspot was delimited by the analysis, corresponding to predation traces placed near the umbo, that is, wall drilling (Fig. 3A). In the point pattern of predation traces derived from feeding trials of N. unifasciata preying on I. subrugosa, we found a single area of higher concentration of drillholes as well as a single hotspot placed at the umbo (Fig. 3B) but slightly closer to the shell edge than in the beach-collected samples. In both beach-collected samples and feeding trial samples, the areas of higher concentration of drillholes covered a relatively small area of the prey skeleton when compared with those from fossil L. latilirata prey (Fig. 2).

Distance-based Statistics and Hypothesis Testing.—The graphical output of the L-function indicated a significant aggregated pattern of predation traces across a range of distances on the skeleton of L. latilirata in both late Pliocene and early Pleistocene point patterns (Fig. 4A, B) (GoF: u = 0.001; p << 0.01). The average number of predation traces within a distance r shorter than 0.5 Bookstein shape units of another location (derived from Bookstein registration) was statistically greater than that expected for a complete random distribution. In addition, there was a reduction in the point of MCD from 0.25 to 0.19 Bookstein shape units for the late Pliocene to the early Pleistocene, respectively. The empirical PCF and L-function calculated for the middle Pleistocene point pattern remained entirely on the simulated envelopes, and thus there was insufficient evidence against the null hypothesis of CSR at any distance on the prey skeleton, a result supported by the GoF test (u = 0.001, p = 0.03).

The graphical output of the L-function calculated for the point pattern derived from the beach-collected samples of I. subrugosa indicates the occurrence of two segregated clusters, that is, a mixed pattern that includes aggregation and segregation at different distances on the prey skeleton (Fig. 4D) (GoF test for CSR: u = 0.005; p = 0.01). Predation traces located at short distances from other traces in the prey skeleton, either pairs of traces located on the umbo or the ventral margin, define the aggregated pattern. Pairs of predation traces located at large distances from each other, specifically those pairs that include one trace located on the umbo and another on the ventral margin, were located farther apart than expected under CSR, supporting a segregated pattern. In this case, the distance between traces approximates the length of the dorsoventral axis (Fig. 4D). In contrast, only an aggregated pattern of traces was observed for feeding trial data of N. unifasciata preying on I. subrugosa (GoF: u = 0.012; p << 0.01); thus, the segregated pattern observed in the beach-collected samples of the same prey species was not supported (Fig. 4E). Remarkably, the MCD in the point patterns derived from the beach-collected samples (0.15) and feeding trial samples (0.17) was similar.

Results from the L-function and PCF for the assembled point patterns are related but appear to capture different aspects of the distribution of predation traces over a range of distances on the prey skeleton (Fig. 4A–J). The L-function was more conservative, indicating a significant segregated pattern of traces at large distances (>0.45 Bookstein shape distance) exclusively in the beach-collected samples (Fig. 4D), whereas the PCF showed additional minor excursions below the confidence envelope between 0.5 and 0.6 Bookstein shape distances, suggesting a weak aggregated pattern in the late Pliocene (Fig. 4F), early Pleistocene (Fig. 4G), and feeding trial (Fig. 4J) point patterns. The PCF also exhibited multiple peaks indicating aggregation of traces across the prey skeleton that could be attributed to multiple predators preying upon the examined prey (Fig. 4F–G,I–J), but a satisfactory interpretation of such peaks is not straightforward. In contrast, the L-function clearly identified two segregated clusters of traces on the beach-collected samples separated by more than 0.5 Bookstein shape distance, roughly corresponding to the distance separating the two areas with high density of traces in the KDE map (Fig. 2A) and corresponding to edge- and wall-drilling modes of predation (Rojas et al. 2015). The PCF indicated that such a pattern emerges at shorter distances (Bookstein shape distances > 0.3), which is consistent with segregated clusters located relatively closer to each other, either on the shell wall or shell edge, an observation that is not supported by KDE and hotspot maps. The NNI (Fig. 4F–J), previously used by Casey et al. (2015) to study site selectivity in drilling predation, provided little information on the spatial distribution of predation traces across different distances on the prey skeleton. The NNI did not capture the segregated pattern of traces in beach-collected samples (Fig. 4N) and even contradicted results that are consistent across the different distance-based functions implemented in our study, for example, the NNI indicated significant segregation of traces (NNI > 1; p < 0.01) for the middle Pleistocene point pattern, which was found to be random by both L-function and PCF.

Figure4.

Graphical output from the distance-based statistics estimated on the fossil samples, modern beach-collected samples, and laboratory feeding trials. A–E, L-function. Black arrow indicates the point of maximum clustering distance (MCD). F–J, Pair correlation function (PCF). K–O, Histograms of nearest neighbor distances and estimated nearest neighbor index (NNI) for the actual data. Legend: Empirical (“Data”) and expected (complete spatial randomness, “CSR”) functions. The dark gray area is the simulation envelope for 999 Monte Carlo simulations of CSR. z, nearest neighbor z-score; a negative z-value indicates aggregation; a positive z-value indicates segregation.

Discussion

Spatial Patterns in the Distribution of Predation Traces on the Prey Skeleton.—Kernel density maps show significant variation in the concentration of predation traces on the skeleton of fossil L. latilirata. The estimated L-function indicates that naticid predation traces on late Pliocene and early Pleistocene samples, which likely included the abundant Naticarius canrena (Linnaeus, 1758), Naticarius plicatella (Conrad, 1863), and Neverita duplicata (Say, 1822) found in the same surveyed bulk samples (Supplementary Table S2), are significantly aggregated. Hotspots of predation traces in the late Pliocene and early Pleistocene spatial point patterns are unique, restricted to wall-drilling cases, and range in location from the middle region of the shell to the umbo. Predation traces in the middle Pleistocene point pattern were randomly placed (Figs. 4C,H). However, because naticid drilling predation tends to be a stereotyped behavior (Kabat 1990; Kingsley-Smith et al. 2003; Rojas et al. 2015), this pattern is likely due to small sample size. The observed reduction in the MCD from late Pliocene to early Pleistocene point patterns for L. latilirata may represent either temporal increase in site selectivity or changes in predator assemblages. The diversity of naticid species recorded in samples from the surveyed geologic units decreases from 19 to 8 between the late Pliocene and early Pleistocene samples, with at least five species in common (Supplementary Table S2). Although we cannot distinguish the individual contribution of each predator to the observed pattern in concentration of predation traces, we observed a temporal shift in hotspot locations migrating toward the ventral area of the shell. Demonstrating temporal shifts in site-selectivity preferences of drilling predators requires replicated samples through several stratigraphic sections, a task that is beyond the scope of this paper. Instead, this case study illustrates how MCD may provide information on temporal changes in site-selectivity preferences. Our results also indicate that L. latilirata was not subject to significant edge-drilling attacks, pointing to the rarity of such predatory behavior in the examined fossil samples.

Although only a single hotspot characterizes the spatial distribution of predation traces from unknown naticid predators preying on I. subrugosa in the beach-collected samples, distance-based statistical analysis indicates cluster segregation, that is, predation traces are significantly aggregated at small scales or short distances and yet simultaneously significantly distributed in a segregated pattern at large distances on the prey skeleton (Fig. 4D). The segregated pattern results from the high concentration of predation traces located at short distances, either on the umbo or at the ventral edge, which gives an MCD of 0.25 Bookstein shape units. In contrast, the segregated pattern results from the concentration of predation traces at a larger distance on the prey skeleton, which approximates the distance (0.45 Bookstein shape units) separating the umbo and the ventral edge (i.e., dorsoventral axis) (Fig. 4D). The distance on the prey skeleton at which the PCF indicates clustered aggregation (Fig. 4I) is significantly shorter than the dorsoventral axis, pointing to the advantage of the L-function over the PCF in characterizing distribution of predation traces. This case study also highlights the need for combining density estimation, hotspot mapping, and distance-based functions to best characterize the distribution of predation traces on prey.

A previous qualitative analysis of site-selectivity patterns between the studied beach-collected and laboratory feeding trial samples used in this study suggested the occurrence of similar drilling patterns (Rojas et al. 2015). Despite similarities (e.g., occurrence of edge drilling, similar MCD), distance-based statistical analysis indicated that the number of predation traces located on the ventral edge of I. subrugosa in the feeding trials is not enough to form a significant segregated pattern as was the case for the beach-collected samples. This evidence may point to biased site-selectivity patterns either due to artifacts under laboratory conditions (e.g., lack of competition may reduce incidence of edge-drilling behavior; Dietl 2004) or the influence of sample size (e.g., small sample sizes may not resolve minor peaks in density of predation traces), which should be further explored. This case study also illustrates the value of the SPPAT framework over qualitative assessments of site selectivity in evaluating alternative modes of drilling predation (i.e., edge and wall drilling).

Site-Selectivity Assessment Using SPPAT.—Current methods to characterize site-selectivity patterns of drilling predators exclude some spatial information inherent to drillhole location, employing an oversimplified description of the trace position, for example, distance of the drillhole to a selected feature in the skeleton (Kowalewski 2002). The SPPAT approach implemented here exploits the spatially explicit nature of the drilling traces to accurately map predation traces on both modern and fossilized shelled invertebrate prey to provide visual and quantitative information on subtle differences in site-selectivity patterns that are not apparent to the naked eye (see Illian et al. 2008; Mitchell and Butterfield 2018; Mitchell et al. 2018), including recognition of the slight migration of hotspots along particular regions of the prey skeleton that otherwise would be described ambiguously (e.g., shell wall and shell edge) or that depend critically on the arbitrary grid system used by sector-based approaches (Kowalewski 1990, 2002). Limitations of sector-based methods for spatial analysis have been recognized previously in other research fields, such as archaeology (Hodder and Okell 1978; Johnson 1984) and ecology (Culley et al. 1933; Elliott 1971). Because distance-based statistics for hypothesis testing are included (Ripley's K-function, L-function, and PCF), the SPPAT approach quantifies complex patterns that emerge at different distances on the prey skeleton (e.g., mixing of aggregated and segregated patterns) that cannot be captured using either sector-based methods (Hodder and Okell 1978) or nearest neighbor statistics (Li and Zhang 2007). The nearest neighbor approach implemented by Casey et al. (2015) is useful for testing the CSR hypothesis and identifying local associations but not for quantifying scale-dependent patterns. Our results support the notion that the NNI is unreliable in the recognition of complex spatial patterns in the distribution of predation traces on shelled invertebrate prey. Li and Zhang (2007) and Mitchell and Butterfield (2018) discussed other examples for which nearest neighbor statistics were also found to be unreliable in the study of spatial patterns.

Although we illustrated the SPPAT approach through case studies of naticid gastropod predation on bivalve prey, it could be extended to examine the composite pattern of predatory traces produced by co-occurring predators (e.g., naticid and muricid gastropods) relative to each group. In addition, drillhole diameter or any other attribute (predator taxonomy, length, etc.) of the predation traces could be used to create a marked point pattern and thus to explore its effect on the spatial patterns.

Conclusions

Site-selectivity analysis of shell-drilling predators is a spatially explicit question. Taking advantage of the spatial information inherent to drillhole location, the SPPAT approach developed here captures diverse aspects in the spatial distribution of predation traces on modern and fossilized prey remains and introduces a novel set of graphic (i.e., kernel density and hotspot maps) and quantitative tools (i.e., K-function, L-function, MCD), largely used in research fields other than paleobiology, for evaluating changes in site-selectivity patterns across taxa, space, and time. Despite limitations affecting data used here—both stratigraphic (i.e., long and unequal time intervals) and geographic (i.e., sampling localities relatively far apart from each other and potential variation in sampled environments)—we showed how SPPAT can be used to visualize areas of high concentration of predation traces on the prey skeleton and to quantify aggregated and segregated patterns occurring at different spatial scales (or distances on the prey skeleton). The visual information on density of predation traces and hotspots, as well the identification of mixed spatial patterns, could not have been obtained using current approaches that lack spatial information and that are dependent on the scale of analysis. We also showed that the point of MCD is an informative measure to examine large-scale temporal and small (intravalve) spatial-scale variation in site-selectivity patterns. Future studies could investigate variation in hotspot locations and MCD in terms of large-scale geographic variation (e.g., latitudinal trends), co-occurring predators (e.g., naticids and muricids), and prey variables, such as taxonomy, thickness, size, mode of life, and ornamentation that are important to this predator–prey interaction. Our study provides a general conceptual framework to orient future research on site-selectivity patterns in drilling predation. However, the SPPAT approach may be generalized to a wide spectrum of ecologic and taphonomic data preserved on skeletal remains of Phanerozoic marine invertebrates, including encrustation, bioerosion traces, and repair scars (Kowalewski 2002) that are spatially explicit and can be mapped.