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 Reference Bromley1981) (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 (Reference Lyons1991) and Ward et al. (Reference Ward, Bailey, Carter, Horton and Zullo1991). Abbreviations: DF, drilling frequency.

Second, we reexamined data from Rojas et al. (Reference Rojas, Hendy and Dietl2015) 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. Reference Ottens, Dietl, Kelley and Stanford2012) 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. (Reference Perea, García, Acero, Valerio and Goméz2008). All images were oriented such that the anteroposterior axis was horizontal to facilitate consistency in landmark data collection (Kolbe et al. Reference Kolbe, Lockwood and Hunt2011). The positions of the predatory drillholes were quantified using a two-dimensional morphometric approach proposed by Roopnarine and Beussink (Reference Roopnarine and Beussink1999) and adapted by Rojas et al. (Reference Rojas, Hendy and Dietl2015) 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. (Reference Casey, Farrell, Dietl and Veilleux2015) 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 (Reference Roopnarine and Beussink1999). Instead, they reflect the general external morphology of the prey skeleton, which is an important consideration for drilling predators when manipulating their prey (Kabat Reference Kabat1990; Kingsley-Smith et al. Reference Kingsley-Smith, Richardson and Seed2003; Rojas et al. Reference Rojas, Hendy and Dietl2015). Landmark coordinates obtained for left valves were inverted to compare directly with the right valves following Kowalewski (Reference Kowalewski2004). The Bookstein baseline registration method for two-dimensional data (Bookstein Reference Bookstein1986) 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. Reference Bivand, Pebesma and Gómez-Rubio2008) (Fig. 1). See Supplementary Table S1 for point-level attributes. The R functions for landmark-based morphometrics of Claude (Reference Claude2008) were used to carry out the morphometric analysis.

Figure 1.

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 subpatterns (sensu Baddeley et al. Reference Baddeley, Rubak and Turner2016) 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.

Figure 2.

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.

Figure 3.

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.

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 (Reference Silverman1986: p. 76, Eq. 4.5), the kernel density was mathematically expressed as:

(1)$$f\lpar x \rpar = \displaystyle{1 \over {nh}}.\sum\nolimits_{i = 1}^n {} \displaystyle{{K\lpar {x-X_i} \rpar } \over h}$$

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

$$K\lpar x \rpar = \displaystyle{3 \over 4}\lpar {1-x^2} \rpar \comma \;\;\vert x \vert \le 1$$

(2)$$K\lpar x \rpar = 0\comma \;x \gt 1.$$

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 (Reference Kowalewski, Kowalewski and Kelley2002: p. 19, Eq. 5):

(3)$${\rm DF}\lpar l \rpar = \displaystyle{{\;n_l} \over {0.5m_l}}\comma \;\;m_l \ge 1$$

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 Reference Fieberg2007). 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.

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. Reference Fotheringham, Brunsdon and Charlton2000). 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 Reference Baddeley and Turner2005). 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 Reference Nelson and Boots2008).

Table 2.

Descriptive statistics used to calculate the optimum bandwidth (h _opt).

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.

We used Ripley's K-function (Ripley Reference Ripley1976, Reference Ripley1977) 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. Reference Fotheringham, Brunsdon and Charlton2000; Baddeley et al. Reference Baddeley, Rubak and Turner2016). Following Ripley (Reference Ripley1988), the estimate of K(r) employed in this study was mathematically expressed as:

(4)$$K\lpar r \rpar = \displaystyle{a \over {n\lpar {n-1} \rpar }}\;\sum\nolimits_i {} \sum\nolimits_j {} I\lpar {d_{ij} \le r} \rpar e_{ij}$$

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. Reference Baddeley, Rubak and Turner2016). 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. Reference Baddeley, Rubak and Turner2016), including both temporal and spatial comparisons. To facilitate the visual assessment of the graphic output, we used Besag's (Reference Besag1977) L-function, which transforms the K-function K _pois(r) = πr2 to the straight-line L _pois (r) = r:

(5)$$L\lpar r \rpar = \sqrt {\displaystyle{{k\lpar r \rpar } \over \pi }} .$$

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. Reference Mullins, Van Ert, Hadfield, Nikolich, Hugh-Jones and Blackburn2015). Because the K-function is cumulative and the graphic tool used for statistical inference may be difficult to interpret (Baddeley et al. Reference Baddeley, Rubak and Turner2016), 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. Reference Illian, Penttinen, Stoyan and Stoyan2008). We contrasted our results against the nearest neighbor index (NNI) (Clark and Evans Reference Clark and Evans1954), a procedure that was previously used by Casey et al. (Reference Casey, Farrell, Dietl and Veilleux2015) 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 (Reference Loosmore and Ford2006). The functions used for spatial statistical analyses are available in the R package Spatstat v. 1.23-4 (Baddeley et al. Reference Baddeley, Rubak and Turner2016).