Environmental correlates of molluscan predator–prey body size in the northern Gulf of Mexico

Abstract

The Mississippi River watershed drains 40% of the continental United States, and the tremendous primary productivity in the adjacent north-central Gulf of Mexico has created one of the most extensive dead zones on Earth. In contrast, smaller watersheds deliver fewer nutrients to the northeastern gulf, and consequently, productivity is limited and hypoxia is uncommon. How has variation in primary productivity, oxygen availability, and sea-surface temperature affected coastal food webs? Here, we investigate environmental controls on the size of molluscan predators and prey in the northern Gulf of Mexico using Holocene death assemblages. Linear mixed models indicate that bivalve size and the frequency of drilling predation are affected by dissolved oxygen concentrations; drilling frequency declines with declining oxygen, whereas bivalve size increases. In contrast, sea-surface temperature is positively associated with the size of molluscan predators and prey. Net primary productivity contributes relatively little to predator or prey size, and predator-to-prey size ratios do not vary consistently with environmental conditions across the northern gulf. Larger bivalves in areas of oxygen limitation may be due to decreased predation pressure and, consequently, greater prey longevity. The larger size of bivalves and predatory gastropods in warmer waters may reflect enhanced growth under these conditions, provided dissolved oxygen concentrations exceed a minimum threshold. Holocene death assemblages can be used to test long-standing hypotheses regarding environmental controls on predator−prey body-size distributions through geologic time and provide baselines for assessing the ongoing effects of anthropogenic eutrophication and warming on coastal food webs.

Non-technical Summary

The Mississippi River delivers tremendous amounts of freshwater and nutrients to the northern Gulf of Mexico, which results in the explosive growth of phytoplankton populations that are typically nutrient limited. Decomposition of phytoplankton blooms by aerobic bacteria can deplete oxygen concentrations in coastal systems, leading to the establishment of oxygen-limited “dead zones.” Changes in the availability of food and dissolved oxygen, as well as changes in temperature, can have wide-reaching effects on coastal food webs. Here, we investigate how primary productivity, dissolved oxygen, and sea-surface temperature affect the sizes of molluscan predators and prey in the northern Gulf of Mexico using collections of shells preserved in seafloor sediment on the continental shelf. We find that the size of bivalves, and the frequency of predatory drilling by snails, are most affected by dissolved oxygen: prey size increases and drilling predation decreases with decreasing concentrations of dissolved oxygen. Sea-surface temperature is positively associated with the size of both molluscan predators and prey. In contrast, net primary productivity has little direct association with size, and the predator-to-prey size ratio also does not vary consistently with environmental conditions in the northern gulf. Larger bivalves in areas with lower oxygen could reflect reduced pressure from predators and, consequently, greater life spans. Larger predator and prey sizes in warmer waters may reflect more optimal conditions for growth. The shells of recently deceased bivalves, and the associated traces of drilling snails on those shells, can be used to investigate long-standing hypotheses about the roles of environmental variation in body-size evolution through geologic time. Furthermore, future studies comparing these historical data with data from present-day communities may help us understand how coastal food webs are changing in response to various human activities.

Introduction

Photosynthetic algae and other marine primary producers transform light into chemical energy that consumers, like suspension-feeding bivalves, use for growth, reproduction, and other metabolic processes. Transfer efficiency, or the percentage of energy that flows from one trophic level to the next, is estimated to range from 5% to 15%, depending on the ecosystem and estimation approach (Pauly and Christensen 1995; Coll et al. 2008; Libralato et al. 2008). Energy availability and transfer efficiency can strongly influence food web structure and the characteristics of predator and prey populations (Ward and Follows 2016; Eddy et al. 2021) and have been hypothesized to be major drivers of macroevolutionary patterns in the fossil record (e.g., Van Valen 1976).

Changes in the structure of marine ecosystems through geologic time have been attributed by a number of authors to increases in primary productivity, shifts in the composition of the dominant primary producers, and changes in transfer efficiency (e.g., Vermeij 1977, 1987; Bambach 1993; Katz et al. 2004; Knoll and Follows 2016). For example, as primary production increased through the Phanerozoic, as evident by a concurrent increase in phytoplankton diversity and abundance, the fossil records of marine consumers also show increased taxonomic and functional diversity (e.g., Bush et al. 2007; Villéger et al. 2011; Martin and Servais 2020), greater biomass (e.g., Bambach 1993; Martin 1996; Payne and Finnegan 2006; Finnegan et al. 2011; Allmon and Martin 2014), and a shift toward more active life habits (e.g., Vermeij 1977; Signor and Brett 1984; Bush et al. 2007). In a recent study of predator and prey body-size trends, Klompmaker et al. (2017) measured the sizes of fossil shells in two clades of basal consumers (brachiopods and mollusks) and associated traces of drilling predation and documented an increase in predator-to-prey size ratios through the Phanerozoic. This gradual change was caused by an increase in the size of predators, which the authors suggest was triggered by the evolution of more nutritious or more abundant prey or predation among predators (Kowalewski et al. 2005; Klompmaker et al. 2017).

Analyzing size structure in present-day marine ecosystems, for which direct measurement and/or remote sensing environmental data are available, can provide quantitative estimates for the relationships between environmental conditions and body size that are rarely possible in the geologic past. However, studies of body-size data that postdate the onset of human exploitation of marine resources must also consider the direct and indirect effects of human activities on size distributions (e.g., Erlandson et al. 2008; Espinosa et al. 2009; McClenachan 2009; Kusnerik et al. 2018). All told, relatively few studies have examined how the body sizes of marine invertebrates vary along present-day primary productivity gradients. Those that have focused primarily on zooplankton (e.g., Stemberger and Gilbert 1985; Gliwicz 1990), though see Roy et al. (2000), McClain et al. (2012), and Berke et al. (2013) for several studies that examined size variation in marine mollusks. McClain et al. (2012) reported that carbon flux—a measure of available chemical energy—was a significant positive predictor of mollusk size across the Northeast Pacific and Northwest Atlantic. Berke et al. (2013), however, found no association between bivalve body size and environmental factors such as net primary productivity and sea-surface temperature in either the Atlantic or Pacific Ocean. These contrasting results warrant further investigation into the relationships between the body size of marine consumers and environmental variation in present-day marine ecosystems.

The Gulf of Mexico is well suited for investigating these themes due to the pronounced primary productivity gradient that occurs in the northern gulf (Lohrenz et al. 1997; Caffrey et al. 2014; Darnell 2015). Geographic variation in primary productivity is heavily influenced by freshwater input from the continental United States (Elliott et al. 2012). Several large rivers flow into the north-central Gulf of Mexico, including the Mobile and Tensaw Rivers, which drain into Mobile Bay in Alabama, and the leading driver of the gulf's energy gradient, the Mississippi River, which flows into the gulf in southern Louisiana. The Mississippi River delivers approximately 1.03 × 10¹¹ m³ of freshwater per year to the northern gulf (Moody 1967), and freshwater discharge is projected to increase over the coming century (Tao et al. 2014; Sinha et al. 2017). Elements essential for plant growth, such as nitrogen and phosphorus, are transported to the northern gulf by these river systems and contribute to the growth of marine phytoplankton and organic carbon accumulation in marine sediments (e.g., Turner and Rabalais 1991; Rabalais et al. 2007; Steckbauer et al. 2011; Breitburg et al. 2018). Primary productivity in the northern gulf generally decreases with distance from the Mississippi River delta, such that the coast of Louisiana has more available chemical energy than coastal environments offshore Florida (Elliott et al. 2012); data averaged from 2003 through 2018 reveal that mean annual primary productivity offshore Louisiana is 2.2 times greater than offshore north Florida (Table 1). While freshwater delivery of nutrients to the northern gulf occurs naturally, human activities have resulted in further nutrient enrichment, which has triggered a cascading array of ecological changes.

Table 1. Geographic coordinates of each sampling station and associated environmental data. NPP, net primary productivity; DO, dissolved oxygen concentration; SST, sea-surface temperature.

Synthetic, chemical fertilizers were first introduced to the United States in the 1930s (Turner and Rabalais 1991), and their use in the United States has risen dramatically since the 1950s (Cao et al. 2018). Increases in marine primary productivity coincide with trends in fertilizer use, as evidenced by elevated concentrations of biologically bound silica and organic carbon since the 1950s in sedimentary cores from regions of the continental shelf adjacent to the Mississippi River delta (Turner and Rabalais 1994; Rabalais et al. 2007). Nutrient enrichment in the northern gulf has led to seasonal bouts of hypoxia that can be traced back to at least the 1970s (Rabalais and Turner 2001; Osterman et al. 2005; Brunner et al. 2006). Decomposition of algal blooms is mediated by aerobic bacteria, which can result in oxygen depletion in bottom waters, termed “dead zones” (Diaz and Rosenberg 2008; Breitburg et al. 2018). The occurrence and spatial extent of hypoxia varies with river discharge and is also affected by thermal stratification of the water column and mixing by currents and storms. Consequently, benthic habitats closer to the Mississippi River have substantially lower dissolved oxygen concentrations than areas of coastal Florida that are adjacent to much smaller watersheds (Rabalais et al. 2007; Santema et al. 2015). Studies of sedimentary cores indicate that hypoxic conditions occurred in coastal Louisiana before the twentieth century but have increased in frequency and extent toward the present-day (Osterman et al. 2011). Recent studies suggest that oxygen depletion on large scales can adversely affect marine organisms’ behavior, physiology, and ecology (e.g., Breitburg 2002; Levin et al. 2009). Studies of benthic response to seasonal hypoxia in the northern Gulf of Mexico report decreased species richness, abundance, and biomass due to increased mortality and migration of organisms to more oxygenated waters (Rabalais et al. 2001; Baustian and Rabalais 2009; Rabalais and Baustian 2020). Fossil assemblages from periods of past oxygen depletion show similar reductions in the taxonomic richness and abundance of marine species, albeit over much longer timescales (e.g., Martindale and Aberhan 2017; Tomašových et al. 2020). Analyses of environmental gradients in the northern Gulf of Mexico, therefore, provide a system that can offer insights into the drivers structuring the body sizes of predators and prey in coastal food webs.

Temperature may also play a critical role in shaping body-size distributions. Increasing temperatures can result in elevated growth rates and, consequently, larger shell sizes due to enhanced enzymatic activity and metabolic rates (Lombard et al. 2009; Saulsbury et al. 2019; Rillo et al. 2020; Reed et al. 2021). However, when temperatures exceed an organism's optimal range, growth rates and reproduction can decline (Pörtner and Farrell 2008; Piazza et al. 2020). Environmental limits on organismal growth are shaped not only by temperature, but by the interactions between temperature and oxygen availability; increasing temperature results in greater metabolic demand, which is limited by dissolved oxygen concentrations (Penn et al. 2018; Penn and Deutsch 2022). As such, oxygen availability and sea-surface temperature may both contribute to predator and prey body-size variation across the northern Gulf of Mexico. Determining how trophic size structure varies spatially with temperature, primary productivity, and oxygen availability can inform predictions for future ecological reorganization in the face of ongoing anthropogenic climate change and eutrophication (Abraham et al. 2013; Sinha et al. 2017; Breitburg et al. 2018; Penn and Deutsch 2022).

The skeletal remains of calcifying organisms, such as bivalve mollusks, can be preserved on the seafloor for millennia, resulting in time-averaged death assemblages (Kidwell 2007; Kidwell and Tomašových 2013). Radiocarbon dates for empty bivalve shells collected from surficial sediments on the inner shelf in the north-central Gulf of Mexico can span the last 3000 years or more, with a median age of approximately 600 years before present that predates the Industrial Revolution (Harnik et al. 2017). Drill holes of predatory gastropods can be preserved in the shells of their bivalve prey, and drill-hole diameter is proportional to the size of the predator (Kitchell et al. 1981; Palmer 1988; Kelley and Hansen 2003; Kowalewski 2004; Klompmaker et al. 2017). Therefore, molluscan death assemblages from the northern Gulf of Mexico can be used to investigate variation in the body sizes of predators and prey along broadscale environmental gradients in coastal ecosystems. These data can provide a historical baseline for body-size variation that predates particular anthropogenic stressors, such as the onset of industrial fishing and synthetic fertilizer use. Here we use Holocene death assemblages from across the northern Gulf of Mexico to investigate the environmental drivers of size variation in bivalve prey and their drilling predators. We hypothesized that the sizes of predators and prey would increase with available chemical energy (primary productivity) and sea-surface temperature, provided there was sufficient oxygen available to meet organisms’ physiological demands. Similar to Klompmaker et al. (2017), we hypothesized that elevated primary productivity off the coast of Louisiana would be associated with larger prey, larger predators, and greater predator-to-prey size ratios and that these size measures would positively correlate with primary productivity across the northern gulf.

Data and Methods

Sample Collection of Death Assemblages

Samples of molluscan death assemblages were collected at 15 stations on the continental shelf offshore Louisiana, Alabama, and Florida during the summers of 2018 and 2019 (Table 1, Fig. 1). Five stations were sampled in each region, each spaced 10 km apart along the −20 m isobath. All stations are located in U.S. federal waters, 9 or more nautical miles from shore. This depth on the continental shelf was selected in order to investigate regional-scale variation in environmental conditions while minimizing transit time on research vessels; environmental conditions at more inshore localities can be heavily influenced by point source inputs, which were not the focus of the present study. Stations in each region were numbered from west to east, 21 to 25 (e.g., LA21 is the westernmost station in Louisiana; AL25 is the easternmost station in Alabama). The top 10–30 cm of seafloor sediment was collected using box corers (Louisiana and Alabama) and a Van Veen sediment grab (Florida). Samples were sieved through a 2-mm mesh to remove smaller individuals, including some juveniles, and all bivalve specimens—both whole shells and shell fragments—greater than 2 mm were retained. Bivalve specimens were then separated into live and dead individuals, with dead individuals used for this study.

Figure 1. Location map. White hatched rectangle on the inset map of the contiguous United States denotes the location of the study area in the northern Gulf of Mexico. Detailed map includes 15 stations (open circles) located on the continental shelf offshore Louisiana, Alabama, and Florida. In each region, five stations spaced approximately 10 km apart were sampled along the −20 m isobath; stations are labeled from west to east, 21 to 25 (e.g., LA21 is the westernmost, and LA25 the easternmost, of the five stations sampled in Louisiana); see Table 1 for station geographic coordinates and associated environmental data. Ten-meter bathymetric contours (dotted gray lines) are plotted out to the continental shelf break (−200 m). Major river systems in the detailed land map are plotted as dashed gray lines.

Body-Size Data

Body-size data were collected from whole bivalve shells in sampled death assemblages. Data were collected from one sample of the death assemblage per station, with the exception of LA24, where two samples were analyzed due to limited numbers of specimens. Six bivalve genera were selected for study based on their presence in samples from across the northern gulf, their range of feeding ecologies, and family-level diversity. Drilling intensity can vary among congeneric species (e.g., Martinelli et al. 2015); however, we chose to focus on the genus level so that body-size data for individual genera extended across the pronounced environmental gradients that exist in the northern gulf, along which species-level turnover occurs. As discussed later, variation among genera was accounted for using linear mixed-effects models. The six genera we considered include a corbulid (Caryocorbula), two venerids (Chione, Lirophora), two protobranchs (Nucula and Nuculana, in the families Nuculidae and Nuculanidae, respectively), and a lucinid (Radiolucina). In total, these genera comprise roughly 16–45% of bivalve shells with intact umbos in our samples; the median and interquartile range (IQR) for Louisiana, Alabama, and Florida were 45% (36–46%), 37% (15–38%), and 16% (13–17%), respectively. These values are approximate, as the total number of bivalves in each sample includes specimens with intact umbos that may not be identifiable to the genus level. Genera were identified using taxonomic descriptions in Mikkelsen and Bieler (2008) and Tunnell et al. (2010) and were cross-checked against the World Register of Marine Species (WoRMS) Database (https://www.marinespecies.org). We note that some of these taxa in the northern gulf (e.g., corbulids) have not been the subject of recent taxonomic work, that differences in taxonomic frameworks exist among researchers, and that uncertainty can also occur when identifying small or juvenile specimens. For example, in Florida, some corbulids with very similar morphologies have been interpreted as either mature Caryocorbula contracta or juvenile Corbula dietziana; following Mikkelsen and Bieler (2008), we placed those specimens in Caryocorbula. Similarly, juvenile Chionopsis intapurpurea can resemble juvenile Chione elevata, and some may have been identified as Chione in our dataset. All six genera are shallow infaunal organisms, based on burrowing depths reported in the literature for these or confamilial genera and based on their body size (Stanley 1970; Mikkelsen and Bieler 2008; Leonard-Pingel and Jackson 2013). We have also observed live individuals of each of these genera in the upper 10–30 cm of sediment at one or more of our stations in the northern gulf. There is considerable dietary breadth among these six genera, including suspension feeders (Caryocorbula, Chione, Lirophora), mixed deposit–suspension feeders (Nucula, Nuculana), and chemosymbiotic feeders (Radiolucina) (Mikkelsen and Bieler 2008; Tunnell et al. 2010). These genera come from across the bivalve tree of life; the last common ancestor for these five extant families dates to the early Paleozoic (Crouch et al. 2021).

Shell length and width were measured using a Leica M125 stereomicroscope and the Leica Application Suite X software, resulting in a total sample size of 3568 shells (Table 2). Shell area, or prey size, was approximated as an ellipse: area = shell width/2 * length/2 * π; shell width was the anterior–posterior dimension of the shell parallel to the hinge, and length (sometimes referred to as “shell height”) was the dorsal–ventral dimension perpendicular to our width measurement. If a predatory drill hole was present, the drill-hole diameter (predator size) was measured, resulting in a total sample size of 445 drill holes (Table 2). Different taxonomic groups of predators produce morphologically distinct drill holes (Kowalewski 1993; Kelley and Hansen 2003; Klompmaker et al. 2019). For example, naticid gastropods create beveled drill holes in which the outer diameter exceeds the inner diameter, whereas muricid gastropods produce vertical-walled drill holes. To determine the taxonomic affinity of drilling predators in our samples, inner and outer drill-hole diameters were measured and used to distinguish naticid versus muricid predators; <1% of observed drill holes could not be assigned to one of these groups. Instances of failed predation (i.e., incomplete drill holes) were also measured, although these were relatively uncommon (~10% of drill holes). Eight shells contained more than one drill hole, and only the largest diameter drill hole was measured in these instances. The predator-to-prey size ratio was determined by dividing drill-hole area by shell area, following the procedure of Klompmaker et al. (2017). Drilling frequency was calculated for each station by dividing the number of drilled shells by the total number of individuals of the six bivalve genera (i.e., the total number of measured shells divided by two). All body-size data and associated metadata generated in this study are available from the Dryad Digital Repository (https://doi.org/10.5061/dryad.n5tb2rc1h). All specimens are stored in the Paleobiology Lab collections in the Department of Earth and Environmental Geosciences at Colgate University (Hamilton, N.Y.).

Table 2. Summary data for individual stations, ordered from west to east across the northern Gulf of Mexico. Sample size is the total number of bivalve shells of the six focal genera. Drilled shells is the total number of drilled shells, including both complete and incomplete drill holes; instances of failed predation were uncommon (~10% of all drill holes). Drilling frequency is the frequency of shells with predatory drill holes, relative to the total number of bivalve individuals (i.e., the number of bivalve shells divided by two). Median shell area and drill-hole diameter represent median prey and predator size, respectively.

Taphonomic Data

Taphonomic processes have the potential to bias measures of bivalve size and frequencies of drilling predation (Klompmaker et al. 2019 and references therein). To assess variation in taphonomic conditions among stations, we measured the extent of shell fragmentation by randomly subsampling approximately 200 bivalve skeletal elements (whole and fragmentary shells) from death assemblage samples from each station; median sample size per station was 226 elements, with an interquartile range of 211 to 241. Each skeletal element was scored using a three-part scale (Best and Kidwell 2000; Lockwood and Work 2006): 0 = no fragmentation (100% of valve); 1 = minimal fragmentation (>20% of valve remaining); 2 = substantial fragmentation (<20% of valve remaining). The mean fragmentation index was calculated for each station, with higher values indicating more extensive taphonomic damage. To evaluate whether prey size tended to be larger at stations with more extensive taphonomic damage, we used Spearman rank-order correlation to evaluate the relationship between mean fragmentation and median prey size. Using this same approach, we also assessed the association between mean fragmentation and the median predator-to-prey size ratio.

Environmental Data

Values of net primary production (NPP, in mg C/m²/day) were retrieved from the Ocean Productivity project at Oregon State University (Oregon State University 2022), using the Standard VGPM algorithm (Behrenfeld and Falkowski 1997) and the 2160 × 4320 grid size yielding a spatial resolution of 0.16°. Only full years of NPP data collection were considered, averaging the data from 2003 until 2018 for each month, then calculating the mean and standard deviation across months. For each of the 15 sampled stations, we matched the spatial coordinates to the nearest grid point of the NPP data, based on the shortest geographic distance using the World Geodetic System of 1984 (WGS 84) and the R package geosphere v. 1.5-18 (Hijmans et al. 2005, 2021) (Table 1).

Sea-surface temperature (SST, in °C) and dissolved oxygen (DO, in μmol/kg) values were downloaded from the Gulf of Mexico Regional Climatology (v. 2), derived from the World Ocean Database (Garcia et al. 2018; Seidov et al. 2020), at 15 m water depth and 0.1° and 1° spatial resolution, respectively. We matched the 15 stations to these environmental data as described for NPP. The mean and standard deviation of SST and DO values were calculated for the summer months of July, August, and September, across the years of 1955 to 2017 (Table 1). These months had the best geographic coverage for the 15 sampled stations (i.e., minimized the farthest distance between these stations and grid points with environmental data). The mean distance between our stations and grid points with NPP, DO, and SST data was NPP = 5.5 km, DO = 77.9 km, and SST = 6.74 km, respectively. The month of June had larger distances between the environmental data and our studied stations, reaching a maximum distance of 198.4 km between the nearest DO value and station AL21. This month was, therefore, not included in the analysis. As such, we analyze mean JAS (July–August–September) SST and DO data, rather than JJA (June–July–August) climatology used in some other studies to look at summer conditions in the Northern Hemisphere. It is also worth noting that, due to the low spatial resolution of the DO data in the World Ocean Database, a total of four unique DO values were retrieved for the 15 studied stations (Table 1).

Statistical Analyses

Raw data were processed and analyzed using RStudio and the lme4, vegan, geiger, car, and afex packages (Bates et al. 2015; Harmon et al. 2020; Fox et al. 2022; Oksanen et al. 2022; Singmann et al. 2022). Histograms of all size measurements show nonnormal distributions with right skews; therefore, the values for prey size, predator size, and predator-to-prey size ratios were log transformed before analysis. Linear mixed-effects models were used to evaluate relationships between body size (i.e., prey size, predator size, and predator-to-prey size ratio) and environmental gradients in the northern gulf (i.e., longitude, NPP, DO, and SST). Environmental predictors were treated as fixed effects, and genus identity as a categorical, random effect. A mixed-effects model is appropriate here, as we are trying to estimate the overall effects of environmental conditions on the size of bivalve prey and their drilling predators, using data for an ecologically and phylogenetically diverse subset of bivalve taxa (e.g., Peng and Lu 2012; Marescaux et al. 2016). Mixed-effects models are widely used when data are hierarchically structured and exhibit non-independence (e.g., body-size measures for different bivalve genera). In addition to modeling the additive effects of each environmental variable on body size, we also modeled the interactive effect of NPP and DO (NPP*DO) on body size, because these two variables are causally linked by the aerobic decomposition of phytoplankton biomass that results in depleted oxygen concentrations in bottom waters. All environmental predictors were scaled to zero mean and unit variance to allow their relative effects on body size to be assessed on a comparable scale. Environmental predictors of drilling frequency were also assessed using logistic regression mixed-effects models. Model selection was conducted using the Akaike information criterion (AIC), where AIC balances the fit of a given model to the data with the number of parameters of that model (Burnham and Anderson 2002). If the best-supported model had an Akaike weight <0.90, parameter estimates were averaged across the ensemble of models considered using their respective Akaike weights (Burnham and Anderson 2002). Variance inflation factors (VIFs) were calculated for each model to determine the extent of collinearity between predictors included as fixed effects. Parameter estimates can be unstable when predictor variables are highly correlated, and in those instances, it is generally recommended that researchers eliminate one or more predictors (Hocking and Pendleton 1983; Craney and Surles 2002; O'Brien 2007). High VIFs (values > 6) necessitated the removal of longitude as a predictor from all models that included one or more other fixed effects. The set of models that were considered included NPP, DO, and SST as potential environmental predictors of body size.

Results

Individual stations may be affected by local to regional-scale variation in sedimentation and mixing, as well as differences in biological production, and consequently, stations contain variable numbers of whole and fragmentary bivalve shells. Most bivalve skeletal material at the 15 stations in the northern gulf is fragmented; shell fragments comprise 81% of the bivalve skeletal elements that were assigned a fragmentation index score. Some variation in fragmentation was observed among regions; the median fragmentation index at the five stations in Louisiana was 1.46 (IQR = 1.45–1.49), in Alabama it was 1.3 (IQR = 1.3–1.31), and in Florida it was 1.35 (IQR = 1.28–1.43). In total, 3568 whole shells of the focal bivalve genera were identified and measured across the 15 stations (Table 2). The greatest number of specimens occurred in the sample from AL21, representing 19.2% of the total sample size. In contrast, the sample from FL21 yielded the fewest specimens, making up only 2.0%. The relative abundance of genera varied among regions (Fig. 2), but Caryocorbula was the most common genus overall in the dataset. LA24 had the largest median prey size (Fig. 3A), predator size (Fig. 3B), and predator-to-prey size ratio (Fig. 3C). In contrast, LA21, AL25, and FL25 had the lowest median prey size (Fig. 3A), predator size (Fig. 3B), and predator-to-prey size ratio (Fig. 3C), respectively. Across the three regions, drilling frequency was highest offshore Alabama and lowest offshore Louisiana (Fig. 4). Of the 445 drill holes that were observed, 94% exhibited the beveled morphologies characteristic of naticid predators. Driller identity varied little among regions, with beveled morphologies comprising 91%, 95%, and 95% of drill holes in Louisiana, Alabama, and Florida, respectively. Across stations, median prey size and median predator-to-prey size ratios were not correlated with mean fragmentation index (p > 0.05 for both Spearman rank order correlation tests, with rho values of 0.16 and −0.33, respectively), indicating that variation in predator and prey body size among our stations is not readily attributable to taphonomic processes.

Figure 2. Proportional abundance of the six focal bivalve genera in Louisiana, Alabama, and Florida. Proportions were calculated by pooling the abundance data for the five stations in each region.

Figure 3. Geographic variation in prey size (shell area, in mm², A), predator size (drill-hole diameter, in mm, B), and predator-to-prey size ratio (% of shell area drilled, C); all size measures are plotted on a logarithmic scale. Size measurements greater or less than 1.5 times the interquartile range are shown as open circles. Solid gray vertical lines separate the three study regions (Louisiana, Alabama, and Florida).

Figure 4. Regional variation in drilling frequency, measured as the frequency of bivalve shells with predatory drill holes at each of the five sampling stations in each region, relative to the total number of bivalve individuals (i.e., the number of bivalve shells divided by two). Median drilling frequency varies markedly across the northern Gulf of Mexico, with the lowest frequency in coastal Louisiana.

Environmental Gradients

Pronounced environmental gradients exist across the northern gulf (Table 1). Bivariate linear regression models fit to the environmental data for the 15 stations indicate that station longitude is negatively associated with net primary productivity (R ² = 0.89; p < 0.01), and positively associated with dissolved oxygen (R ² = 0.59; p < 0.01) and sea-surface temperature (R ² = 0.26; p = 0.05). In other words, stations in the north-central gulf tend to be characterized by higher primary productivity, lower dissolved oxygen, and lower sea-surface temperatures than stations in the northeastern gulf.

Model Selection Framework

Linear mixed-effects models were used to evaluate hypothesized relationships between environmental gradients and predator and prey body size in the northern Gulf of Mexico, while accounting for genus identity as a random effect. Observed variations in prey size were best described by a mixed model that includes dissolved oxygen (p < 0.01) and sea-surface temperature (p < 0.01) as fixed effects and genus identity as a random effect (AIC = 5042; Akaike weight = 0.787; Table 3). Variations in predator size across the northern gulf were best described by a model that included sea-surface temperature as a fixed effect (p < 0.01) and genus identity as a random effect (AIC = 502; Akaike weight = 0.743; Table 4). Variation in predator-to-prey size ratios was equally well described by two models: one that included dissolved oxygen (p = 0.497) and genus identity, and another that included sea-surface temperature (p = 0.328) and genus identity, although the estimated coefficients for the fixed effects in both of these models were indistinguishable from zero (AIC = 1023; Akaike weight = 0.338; Table 5). The best-supported model for drilling frequency included dissolved oxygen (p < 0.01) and net primary productivity (p = 0.08) as fixed effects, and genus identity as a random effect (AIC = 2471; Akaike weight = 0.292; Table 6).

Table 3. Results for a set of linear mixed-effects models that consider the associations between prey size and environmental variables across the northern Gulf of Mexico. Models are ordered according to their relative support using Akaike weights; the first model is the best-supported model in the set. In all models, genus was included as a random effect. Environmental variables were scaled to zero mean and unit variance to allow their fixed effects on size to be assessed on a comparable scale. Confidence intervals are provided for each fixed effect along with the associated p-value; significance at p = 0.05. AIC, Akaike information criterion; AW, Akaike weight; VIF, variance inflation factor; NPP, net primary productivity; DO, dissolved oxygen concentration; SST, sea-surface temperature; NPP*DO is the multiplicative effect of net primary productivity and dissolved oxygen concentration.

Table 4. Results for a set of linear mixed-effects models that consider the associations between predator size and environmental variables across the northern Gulf of Mexico. Models are ordered according to their relative support using Akaike weights; the first model is the best-supported model in the set. In all models, the taxonomic identity of the prey (genus) was included as a random effect. Environmental variables were scaled to zero mean and unit variance to allow their fixed effects on size to be assessed on a comparable scale. Confidence intervals are provided for each fixed effect along with the associated p-value; significance at p = 0.05. AIC, Akaike information criterion; AW, Akaike weight; VIF, variance inflation factor; NPP, net primary productivity; DO, dissolved oxygen concentration; SST, sea-surface temperature; NPP*DO is the multiplicative effect of net primary productivity and dissolved oxygen concentration.

Table 5. Results for a set of linear mixed-effects models that consider the associations between predator-to-prey size ratio and environmental variables across the northern Gulf of Mexico. Models are ordered according to their relative support using Akaike weights; the first model is the best-supported model in the set. In all models, the taxonomic identity of the prey (genus) was included as a random effect. Environmental variables were scaled to zero mean and unit variance to allow their fixed effects on the predator-to-prey size ratio to be assessed on a comparable scale. Confidence intervals are provided for each fixed effect along with the associated p-value; significance at p = 0.05. AIC, Akaike information criterion; AW, Akaike weight; VIF, variance inflation factor; NPP, net primary productivity; DO, dissolved oxygen concentration; SST, sea-surface temperature; NPP*DO is the multiplicative effect of net primary productivity and dissolved oxygen concentration.

Table 6. Results for a set of binomial linear mixed-effects models that consider the associations between drilling frequency (drilled vs. non-drilled) and environmental variables across the northern Gulf of Mexico. Models are ordered according to their relative support using Akaike weights; the first model is the best-supported model in the set. In all models, the taxonomic identity of the prey (genus) was included as a random effect. Environmental variables were scaled to zero mean and unit variance to allow their fixed effects on drilling frequency to be assessed on a comparable scale. Confidence intervals are provided for each fixed effect along with the associated p-value; significance at p = 0.05. AIC, Akaike information criterion; AW, Akaike weight; VIF, variance inflation factor; NPP, net primary productivity; DO, dissolved oxygen concentration; SST, sea-surface temperature; NPP*DO is the multiplicative effect of net primary productivity and dissolved oxygen concentration.

None of the models relating environmental conditions to the size of bivalve prey, their predators, predator-to-prey size ratios, or drilling frequency, had an Akaike weight greater than 0.90 (Tables 3–6). Therefore, model averaging was used to estimate the fixed effects of each environmental predictor across the ensemble of models considered for each response variable (e.g., prey size; Table 7). Dissolved oxygen had the greatest effect on prey size, with larger bivalves found in lower-oxygen conditions (Table 7, Fig. 5A). Prey size was also positively associated with sea-surface temperature (Fig. 5B), with a relatively similar effect size as dissolved oxygen (Table 7). Net primary productivity was also significantly associated with prey size but had the smallest effect size (Table 7, Fig. 5C). Sea-surface temperature was the only environmental factor that was significantly associated with the size of drilling predators, with larger predators observed in warmer waters (Table 7, Fig. 5D). Dissolved oxygen concentrations had the strongest effect on the frequency of drilling predation, with greater drilling in more oxygenated waters (Table 7). Net primary productivity was also positively associated with drilling frequency, with greater drilling in more productive environments, provided there was sufficient oxygen (Table 6). The effect of net primary productivity on drilling frequency was substantially weaker than that estimated for dissolved oxygen (Table 7). Predator-to-prey size ratios were not significantly associated with any of the environmental predictors.

Table 7. Associations between environmental variables and prey size, predator size, predator-to-prey size ratio, and drilling frequency across the northern Gulf of Mexico. Estimated fixed effects were calculated by model averaging of parameter estimates using the Akaike weight of each model in the set that was considered. Because the environmental variables were scaled to zero mean and unit variance, their relative effects on a given response variable can be assessed on a comparable scale. NPP, net primary productivity; DO, dissolved oxygen concentration; SST, sea-surface temperature; NPP*DO is the multiplicative effect of net primary productivity and dissolved oxygen concentration.

Figure 5. Variation in the size of bivalve prey and their drilling predators along different environmental gradients in the northern Gulf of Mexico, plotted on a logarithmic scale. Statistically significant size–environment associations are plotted in each panel; see Tables 3 and 4 for further detail. Prey size is inversely associated with dissolved oxygen concentration (A), and positively associated with sea-surface temperature (B) and net primary productivity (C). Predator size (D) is positively associated with sea-surface temperature. Solid black dots are the median sizes, and vertical lines are the interquartile ranges. The low spatial resolution of available dissolved oxygen data meant there were a total of four sites in the World Ocean Atlas with dissolved oxygen data proximal enough to be matched to our stations.

Discussion

We used Holocene death assemblages to investigate the environmental drivers of predator and prey body-size distributions along present-day environmental gradients in the northern Gulf of Mexico. By focusing on this system, we were able to use environmental data generated through direct measurement and remote sensing methods to evaluate a series of factors long hypothesized to be important drivers of body-size trends through geologic time (e.g., primary productivity, sea-surface temperature, and dissolved oxygen concentration). In doing so, we assume that environmental data for the northern gulf collected in the late twentieth and early twenty-first centuries provide, at least to a first order, an accurate characterization of geographic differences in environmental conditions faced by predator and prey populations over past centuries. This assumption is reasonable, as watersheds have delivered nutrients to the northern gulf for millennia, and variation in watershed size can serve as a rough proxy for nutrient inputs and resulting coastal primary productivity (Van der Zwaan 2000). The net primary productivity gradient that exists today across the northern gulf also existed historically due to the pronounced differences in watershed size in the region, from the extensive Mississippi River watershed that affects coastal Louisiana and surrounding areas of the north-central gulf, to the smaller Mobile Bay watershed that discharges in coastal Alabama, to the smaller-still Apalachicola watershed in north Florida. Sedimentary cores also document the occurrence of hypoxic conditions in coastal Louisiana over the past 300 years (Osterman et al. 2011). Anthropogenic eutrophication and associated hypoxia have steepened these preexisting environmental gradients, leading to greater primary productivity and more extensive coastal hypoxia over the past century in the north-central gulf (Brunner et al. 2006; Rabalais et al. 2007; Osterman et al. 2011; Rabalais and Baustian 2020).

We hypothesized that basal consumers in coastal food webs would be larger in areas of greater net primary productivity, as greater chemical energy would allow bivalves to invest more energy into growth and reproduction, potentially leading to greater body sizes and/or larger populations. Consequently, increased energy availability would flow up to higher trophic levels, resulting in larger-bodied predators and, potentially, an increase in the ratio of predator-to-prey body size (Klompmaker et al. 2017). We also hypothesized that basal consumers and their predators would be larger in warmer waters due to faster growth rates. This predicted relationship is consistent with the findings of some previous studies (e.g., Lombard et al. 2009; Saulsbury et al. 2019; Rillo et al. 2020; Reed et al. 2021), but not others (e.g., Hunt and Roy 2006), suggesting that further work is needed to understand size–temperature relationships in marine ecosystems. Dissolved oxygen may also play a role in body-size distributions through its effect on food web structure. Due to their higher metabolic demand, predators may be excluded from regions affected by hypoxia (e.g., Rabalais et al. 2001; Breitburg 2002; Levin et al. 2009). Provided there is sufficient oxygen available for basal consumers, predator exclusion has the potential to lead to greater prey longevity and, consequently, larger prey size.

The six bivalve genera that we studied span a range of feeding strategies including suspension feeding, mixed deposit–suspension feeding, and chemosymbiosis with sulfide-oxidizing bacteria. To what extent is net primary productivity likely an important control on food supply, and by extension growth, across these three ecological groups? Suspension-feeding bivalves—such as Caryocorbula, Chione, and Lirophora—filter organic matter from the water adjacent to the sediment–water interface, which may include phytoplankton but also resuspended organic matter (Mikkelsen and Bieler 2008). Mixed feeders (i.e., those capable of deposit feeding and suspension feeding)—such as Nucula and Nuculana—consume a mixture of nonliving organic detritus and sediment-associated microbes (Lopez and Levinton 1987; Mikkelsen and Bieler 2008), but depending on prevailing environmental conditions and the taxonomic group, may also suspension feed (Ward and Shumway 2004). Chemosymbiotic taxa—such as Radiolucina—host sulfide-oxidizing bacteria that are dependent on the decomposition of organic detritus (Taylor and Glover 2010). Studies of gut contents in co-occurring suspension- and deposit-feeding bivalves have documented significant dietary overlap between these groups (Kamermans 1994), and phytoplankton and benthic microalgae both contribute to the organic detritus that chemosymbiotic bivalves depend upon. Given the information currently available, we suggest that net primary productivity, as measured by chlorophyll a concentrations, can serve as a general measure of food availability for this ecologically diverse set of basal consumers in the northern gulf.

We used a mixed-effects modeling framework to disentangle the unique fixed effects of each of these environmental factors on predator and prey size while accounting for genus identity as a random effect. When the environmental fixed effects on predator and prey size were examined independently, each environmental variable seemed to contribute significantly to trophic size structure. However, when coefficients were estimated using multivariate models that accounted for the covariation between predictors, and then averaged over the full ensemble of models, it became clear that dissolved oxygen and sea-surface temperature were important factors associated with the size of bivalve prey and their drilling predators, whereas net primary productivity contributed primarily through its role in the formation and extent of bottom-water hypoxia. Our analysis of drilling frequencies across the northern gulf provides a possible explanation for why prey size is negatively associated with dissolved oxygen concentrations. Drilling frequencies increase significantly with dissolved oxygen and net primary productivity. Despite high primary productivity in areas of coastal Louisiana that are adjacent to the Mississippi River, development of hypoxic conditions limits the activity of drilling predators. In contrast, in coastal Alabama, where primary productivity is high but hypoxia is more limited and spatially and temporally variable (Rabalais and Turner 2001; Brunner et al. 2006), drilling frequencies are at their greatest. Drilling frequencies decline from Alabama to Florida as productivity declines, but dissolved oxygen levels remain stable. These results suggest that hypoxic bottom waters may act as a predation refuge for at least some bivalve taxa that can tolerate such conditions.

Our finding that the body size of predators and prey is positively associated with sea-surface temperature is consistent with some previous studies that have reported a positive relationship between temperature and growth rate in marine bivalves (Saulsbury et al. 2019; Reed et al. 2021) and temperature and cell size in some species of foraminifera (Lombard et al. 2009; Rillo et al. 2020). The sensitivity of taxa to changes in dissolved oxygen depends on individual tolerances, yet in general, physiological stress tends to occur when oxygen is limited (e.g., Rabalais et al. 2001; Vaquer-Sunyer and Duarte 2008; Kuk-Dzul and Díaz-Castañeda 2016). Our finding that the average size of bivalves increased as oxygen concentrations declined is inconsistent with the expectation that growth and metabolic processes slow under hypoxia (e.g., Rabalais et al. 2001), but may reflect the lower metabolic demand of basal versus secondary consumers in coastal food webs. Our finding that primary productivity has little direct effect on bivalve body size is consistent with the results of Berke et al. (2013) and Saulsbury et al. (2019). This pattern is in contrast to the interpretation of some size-related trends in the fossil record that assert that increasing availability of chemical energy at lower trophic levels shapes the size distributions of higher-level consumers (e.g., Bambach 1993; Klompmaker et al. 2017). Klompmaker et al. (2017) documented an increase in predator size and the predator-to-prey size ratio through the Phanerozoic, which they attributed to greater energy availability for higher-level consumers over time due to the evolution of more nutritious prey, an increase in prey population size, or increased predation among predators. Similar geographic trends in predator-to-prey size ratios were not observed in our study when we compared populations in productive habitats adjacent to the Mississippi River with those in less-productive environments offshore coastal Florida. Variation in the size of drilling predators, and the frequency of drilling predation, could potentially reflect spatial variation in the composition of predators (Alexander and Dietl 2001; Mondal et al. 2021); however, this is unlikely to be driving the observed patterns in the northern gulf, as >90% of observed drill holes exhibited morphologies indicative of naticid predation, and predator identity did not vary geographically in our samples.

Why was net primary productivity not more important in the gulf, given the emphasis of previous studies on the role of food availability in structuring size distributions in modern and ancient food webs? Hypoxia provides one possible explanation for this apparent inconsistency. In the northern gulf, drilling frequency is positively associated with both net primary productivity and dissolved oxygen. Low-oxygen conditions, like those present at our Louisiana stations, can lead to decreased species richness, abundance, and biomass due to increased mortality and migration of organisms to more oxygenated waters (Rabalais et al. 2001). Because metabolic rate increases with activity level, drilling predators will tend to have greater aerobic demand than their bivalve prey, which could lead to predator exclusion in hypoxic settings. Hypoxic zones have been hypothesized to be predation refuges for less metabolically active prey species, like bivalves, which allows these organisms to grow to maturity under diminished predation pressure (e.g., Altieri 2008; Galligan et al. 2022). Altieri (2008) reported that the release of top-down control from predators in hypoxic conditions resulted in longer survivorship and a greater number of larger-bodied prey. This form of predatory release could explain why we observed a negative correlation between prey size and dissolved oxygen concentrations. Future studies of bivalve sclerochronology in the northern gulf could further disentangle the environmental controls on bivalve body size in this region. Do larger bivalves occur in shelf environments adjacent to the Mississippi River delta because nutrient delivery in these settings leads to increased primary productivity and elevated bivalve growth rates? Or does predatory release due to hypoxic conditions in the vicinity of the delta increase survivorship of bivalves, such that their larger size reflects greater longevity (i.e., life span)? These hypotheses are not mutually exclusive, and larger prey size may reflect some combination of increased growth rates and greater survivorship.

Radiocarbon dates for shells in surficial death assemblages collected from continental shelf environments can range up to 3000 or more years before present, with a median age that can predate the Industrial Revolution (Kidwell 2013; Kidwell and Tomašových 2013; Harnik et al. 2017). Therefore, the geographic differences we observed in predator and prey body sizes most likely represent spatial variation in trophic size structure before the pronounced increase in anthropogenic nutrient delivery in the 1950s following the development and widespread application of synthetic, chemical fertilizers (e.g., Hasler 1969; Rabalais and Turner 2001; Osterman et al. 2005; Breitburg et al. 2018). Agricultural and urban runoff delivers excess nutrients to coastal ecosystems like those of the northern Gulf of Mexico (Hasler 1969). Current projections indicate that the Mississippi River transports approximately 32 million tons of nitrogen and 7.6 million tons of phosphorus into the northern gulf every year (Darnell 2015). As a result of this nutrient enrichment, phytoplankton normally limited by the availability of nitrogen and phosphorus experience unhindered growth (e.g., Steckbauer et al. 2011; Breitburg et al. 2018). The subsequent decomposition of these autotrophs has depleted dissolved oxygen concentrations in coastal ecosystems, leading to environmental degradation. Our results provide a historical baseline for future studies investigating ecosystem response to anthropogenic eutrophication and hypoxia. Discordance between live and death assemblages can reveal biotic response to recent ecological and environmental changes resulting from human activities (Kidwell 2007; Harnik et al. 2017; Tomašových et al. 2020). Considering the significant increase in nutrient inputs to the northern Gulf of Mexico beginning in the 1950s (Turner and Rabalais 1994; Brunner et al. 2006), we expect live–dead comparisons to reveal significant shifts in trophic size structure and drilling frequency over recent decades.

Conclusions

Our study used Holocene-aged death assemblages to evaluate the environmental controls of bivalve predator and prey body-size distributions across the northern Gulf of Mexico. We found that dissolved oxygen and sea-surface temperature were associated with the size of molluscan predators and prey, while primary productivity contributed little to geographic differences in body-size distributions. Additionally, we found that drilling frequency was limited under hypoxic conditions and increased with oxygen availability and primary productivity. Coastal areas characterized by low-oxygen conditions may serve as a predation refuge for some bivalves, resulting in greater survivorship and, consequently, larger body size. Late Holocene gradients in trophic size structure and drilling frequency provide a historical baseline for future studies evaluating the impacts of anthropogenic nutrient loading on coastal ecosystems in the northern Gulf of Mexico and beyond.