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 Reference Pauly and Christensen1995; Coll et al. Reference Coll, Libralato, Tudela, Palomera and Pranovi2008; Libralato et al. Reference Libralato, Coll, Tudela, Palomera and Pranovi2008). Energy availability and transfer efficiency can strongly influence food web structure and the characteristics of predator and prey populations (Ward and Follows Reference Ward and Follows2016; Eddy et al. Reference Eddy, Bernhardt, Blanchard, Cheung, Colléter, Du Pontavice, Fulton, Gascuel, Kearney and Petrik2021) and have been hypothesized to be major drivers of macroevolutionary patterns in the fossil record (e.g., Van Valen Reference Van Valen1976).
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 Reference Vermeij1977, Reference Vermeij1987; Bambach Reference Bambach1993; Katz et al. Reference Katz, Finkel, Grzebyk, Knoll and Falkowski2004; Knoll and Follows Reference Knoll and Follows2016). 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. Reference Bush, Bambach and Daley2007; Villéger et al. Reference Villéger, Novack-Gottshall and Mouillot2011; Martin and Servais Reference Martin and Servais2020), greater biomass (e.g., Bambach Reference Bambach1993; Martin Reference Martin1996; Payne and Finnegan Reference Payne and Finnegan2006; Finnegan et al. Reference Finnegan, McClain, Kosnik and Payne2011; Allmon and Martin Reference Allmon and Martin2014), and a shift toward more active life habits (e.g., Vermeij Reference Vermeij1977; Signor and Brett Reference Signor and Brett1984; Bush et al. Reference Bush, Bambach and Daley2007). In a recent study of predator and prey body-size trends, Klompmaker et al. (Reference Klompmaker, Kowalewski, Huntley and Finnegan2017) 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. Reference Kowalewski, Hoffmeister, Baumiller and Bambach2005; Klompmaker et al. Reference Klompmaker, Kowalewski, Huntley and Finnegan2017).
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. Reference Erlandson, Rick, Braje, Steinberg and Vellanoweth2008; Espinosa et al. Reference Espinosa, Rivera-Ingraham, Fa and García-Gómez2009; McClenachan Reference McClenachan2009; Kusnerik et al. Reference Kusnerik, Lockwood, Grant, Tyler and Schneider2018). 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 Reference Stemberger and Gilbert1985; Gliwicz Reference Gliwicz1990), though see Roy et al. (Reference Roy, Jablonski and Martien2000), McClain et al. (Reference McClain, Gullett, Jackson-Ricketts and Unmack2012), and Berke et al. (Reference Berke, Jablonski, Krug, Roy and Tomašových2013) for several studies that examined size variation in marine mollusks. McClain et al. (Reference McClain, Gullett, Jackson-Ricketts and Unmack2012) 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. (Reference Berke, Jablonski, Krug, Roy and Tomašových2013), 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. Reference Lohrenz, Fahnenstiel, Redalje, Lang, Chen and Dagg1997; Caffrey et al. Reference Caffrey, Murrell, Amacker, Harper, Phipps and Woodrey2014; Darnell Reference Darnell2015). Geographic variation in primary productivity is heavily influenced by freshwater input from the continental United States (Elliott et al. Reference Elliott, Pierson and Roman2012). 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 × 1011 m3 of freshwater per year to the northern gulf (Moody Reference Moody1967), and freshwater discharge is projected to increase over the coming century (Tao et al. Reference Tao, Tian, Ren, Yang, Yang, He, Cai and Lohrenz2014; Sinha et al. Reference Sinha, Michalak and Balaji2017). 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 Reference Turner and Rabalais1991; Rabalais et al. Reference Rabalais, Turner, Gupta, Platon and Parsons2007; Steckbauer et al. Reference Steckbauer, Duarte, Carstensen, Vaquer-Sunyer and Conley2011; Breitburg et al. Reference Breitburg, Levin, Oschlies, Grégoire, Chavez, Conley, Garçon, Gilbert, Gutiérrez and Isensee2018). 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. Reference Elliott, Pierson and Roman2012); 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 Reference Turner and Rabalais1991), and their use in the United States has risen dramatically since the 1950s (Cao et al. Reference Cao, Lu and Yu2018). 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 Reference Turner and Rabalais1994; Rabalais et al. Reference Rabalais, Turner, Gupta, Platon and Parsons2007). 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 Reference Rabalais and Turner2001; Osterman et al. Reference Osterman, Poore, Swarzenski and Turner2005; Brunner et al. Reference Brunner, Beall, Bentley and Furukawa2006). Decomposition of algal blooms is mediated by aerobic bacteria, which can result in oxygen depletion in bottom waters, termed “dead zones” (Diaz and Rosenberg Reference Diaz and Rosenberg2008; Breitburg et al. Reference Breitburg, Levin, Oschlies, Grégoire, Chavez, Conley, Garçon, Gilbert, Gutiérrez and Isensee2018). 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. Reference Rabalais, Turner, Gupta, Platon and Parsons2007; Santema et al. Reference Santema, Clarke, Speer and Huettel2015). 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. Reference Osterman, Poore, Swarzenski, Hollander, Turner, Buster and Holmes2011). Recent studies suggest that oxygen depletion on large scales can adversely affect marine organisms’ behavior, physiology, and ecology (e.g., Breitburg Reference Breitburg2002; Levin et al. Reference Levin, Ekau, Gooday, Jorissen, Middelburg, Naqvi, Neira, Rabalais and Zhang2009). 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. Reference Rabalais, Harper, Turner, Rabalais and Turner2001; Baustian and Rabalais Reference Baustian and Rabalais2009; Rabalais and Baustian Reference Rabalais and Baustian2020). 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 Reference Martindale and Aberhan2017; Tomašových et al. Reference Tomašových, Albano, Fuksi, Gallmetzer, Haselmair, Kowalewski, Nawrot, Nerlović, Scarponi and Zuschin2020). 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. Reference Lombard, Labeyrie, Michel, Spero and Lea2009; Saulsbury et al. Reference Saulsbury, Moss, Ivany, Kowalewski, Lindberg, Gillooly, Heim, McClain, Payne and Roopnarine2019; Rillo et al. Reference Rillo, Miller, Kucera and Ezard2020; Reed et al. Reference Reed, Godbold, Grange and Solan2021). However, when temperatures exceed an organism's optimal range, growth rates and reproduction can decline (Pörtner and Farrell Reference Pörtner and Farrell2008; Piazza et al. Reference Piazza, Ullmann and Aberhan2020). 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. Reference Penn, Deutsch, Payne and Sperling2018; Penn and Deutsch Reference Penn and Deutsch2022). 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. Reference Abraham, Baringer, Bindoff, Boyer, Cheng, Church, Conroy, Domingues, Fasullo and Gilson2013; Sinha et al. Reference Sinha, Michalak and Balaji2017; Breitburg et al. Reference Breitburg, Levin, Oschlies, Grégoire, Chavez, Conley, Garçon, Gilbert, Gutiérrez and Isensee2018; Penn and Deutsch Reference Penn and Deutsch2022).
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 Reference Kidwell2007; Kidwell and Tomašových Reference Kidwell and Tomašových2013). 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. Reference Harnik, Torstenson and Williams2017). 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. Reference Kitchell, Boggs, Kitchell and Rice1981; Palmer Reference Palmer1988; Kelley and Hansen Reference Kelley, Hansen, Kelley, Kowalewski and Hansen2003; Kowalewski Reference Kowalewski2004; Klompmaker et al. Reference Klompmaker, Kowalewski, Huntley and Finnegan2017). 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. (Reference Klompmaker, Kowalewski, Huntley and Finnegan2017), 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. Reference Martinelli, Kosnik and Madin2015); 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 (Reference Mikkelsen and Bieler2008) and Tunnell et al. (Reference Tunnell, Andrews, Barrera and Moretzsohn2010) 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 (Reference Mikkelsen and Bieler2008), 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 Reference Stanley1970; Mikkelsen and Bieler Reference Mikkelsen and Bieler2008; Leonard-Pingel and Jackson Reference Leonard-Pingel and Jackson2013). 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 Reference Mikkelsen and Bieler2008; Tunnell et al. Reference Tunnell, Andrews, Barrera and Moretzsohn2010). 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. Reference Crouch, Edie, Collins, Bieler and Jablonski2021).
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 Reference Kowalewski1993; Kelley and Hansen Reference Kelley, Hansen, Kelley, Kowalewski and Hansen2003; Klompmaker et al. Reference Klompmaker, Kelley, Chattopadhyay, Clements, Huntley and Kowalewski2019). 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. (Reference Klompmaker, Kowalewski, Huntley and Finnegan2017). 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. Reference Klompmaker, Kelley, Chattopadhyay, Clements, Huntley and Kowalewski2019 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 Reference Best and Kidwell2000; Lockwood and Work Reference Lockwood and Work2006): 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/m2/day) were retrieved from the Ocean Productivity project at Oregon State University (Oregon State University 2022), using the Standard VGPM algorithm (Behrenfeld and Falkowski Reference Behrenfeld and Falkowski1997) 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. Reference Hijmans, Cameron, Parra, Jones and Jarvis2005, Reference Hijmans, Karney, Williams and Vennes2021) (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. Reference Garcia, Weathers, Paver, Smolyar, Boyer, Locarnini, Zweng, Mishonov, Baranova and Seidov2018; Seidov et al. Reference Seidov, Mishonov, Boyer, Baranova, Nyadjro, Cross, Parsons and Weathers2020), 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. Reference Bates, Mächler, Bolker and Walker2015; Harmon et al. Reference Harmon, Pennell, Brock, Brown, Challenger, Eastman, FitzJohn, Glor, Hunt, Revell, Slater, Uyeda and Weir2020; Fox et al. Reference Fox, Weisberg, Price, Adler, Bates, Baud-Bovy and Bolker2022; Oksanen et al. Reference Oksanen, Simpson, Blanchet, Kindt, Legendre, Minchin and O'Hara2022; Singmann et al. Reference Singmann, Bolker, Westfall, Aust, Ben-Shacher, Hojsgaard, Fox, Lawrence, Mertens, Love, Length and Christensen2022). 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 Reference Peng and Lu2012; Marescaux et al. Reference Marescaux, Falisse, Lorquet, Van Doninck, Beisel and Descy2016). 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 Reference Burnham and Anderson2002). 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 Reference Burnham and Anderson2002). 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 Reference Hocking and Pendleton1983; Craney and Surles Reference Craney and Surles2002; O'Brien Reference O'Brien2007). 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 mm2, 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 2 = 0.89; p < 0.01), and positively associated with dissolved oxygen (R 2 = 0.59; p < 0.01) and sea-surface temperature (R 2 = 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. Reference Osterman, Poore, Swarzenski, Hollander, Turner, Buster and Holmes2011). 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. Reference Brunner, Beall, Bentley and Furukawa2006; Rabalais et al. Reference Rabalais, Turner, Gupta, Platon and Parsons2007; Osterman et al. Reference Osterman, Poore, Swarzenski, Hollander, Turner, Buster and Holmes2011; Rabalais and Baustian Reference Rabalais and Baustian2020).
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. Reference Klompmaker, Kowalewski, Huntley and Finnegan2017). 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. Reference Lombard, Labeyrie, Michel, Spero and Lea2009; Saulsbury et al. Reference Saulsbury, Moss, Ivany, Kowalewski, Lindberg, Gillooly, Heim, McClain, Payne and Roopnarine2019; Rillo et al. Reference Rillo, Miller, Kucera and Ezard2020; Reed et al. Reference Reed, Godbold, Grange and Solan2021), but not others (e.g., Hunt and Roy Reference Hunt and Roy2006), 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. Reference Rabalais, Harper, Turner, Rabalais and Turner2001; Breitburg Reference Breitburg2002; Levin et al. Reference Levin, Ekau, Gooday, Jorissen, Middelburg, Naqvi, Neira, Rabalais and Zhang2009). 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 Reference Mikkelsen and Bieler2008). 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 Reference Lopez and Levinton1987; Mikkelsen and Bieler Reference Mikkelsen and Bieler2008), but depending on prevailing environmental conditions and the taxonomic group, may also suspension feed (Ward and Shumway Reference Ward and Shumway2004). Chemosymbiotic taxa—such as Radiolucina—host sulfide-oxidizing bacteria that are dependent on the decomposition of organic detritus (Taylor and Glover Reference Taylor, Glover and Kiel2010). Studies of gut contents in co-occurring suspension- and deposit-feeding bivalves have documented significant dietary overlap between these groups (Kamermans Reference Kamermans1994), 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 Reference Rabalais and Turner2001; Brunner et al. Reference Brunner, Beall, Bentley and Furukawa2006), 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. Reference Saulsbury, Moss, Ivany, Kowalewski, Lindberg, Gillooly, Heim, McClain, Payne and Roopnarine2019; Reed et al. Reference Reed, Godbold, Grange and Solan2021) and temperature and cell size in some species of foraminifera (Lombard et al. Reference Lombard, Labeyrie, Michel, Spero and Lea2009; Rillo et al. Reference Rillo, Miller, Kucera and Ezard2020). 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. Reference Rabalais, Harper, Turner, Rabalais and Turner2001; Vaquer-Sunyer and Duarte Reference Vaquer-Sunyer and Duarte2008; Kuk-Dzul and Díaz-Castañeda Reference Kuk-Dzul and Díaz-Castañeda2016). 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. Reference Rabalais, Harper, Turner, Rabalais and Turner2001), 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. (Reference Berke, Jablonski, Krug, Roy and Tomašových2013) and Saulsbury et al. (Reference Saulsbury, Moss, Ivany, Kowalewski, Lindberg, Gillooly, Heim, McClain, Payne and Roopnarine2019). 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 Reference Bambach1993; Klompmaker et al. Reference Klompmaker, Kowalewski, Huntley and Finnegan2017). Klompmaker et al. (Reference Klompmaker, Kowalewski, Huntley and Finnegan2017) 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 Reference Alexander and Dietl2001; Mondal et al. Reference Mondal, Chakraborty, Saha, Dey and Sarkar2021); 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. Reference Rabalais, Harper, Turner, Rabalais and Turner2001). 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 Reference Altieri2008; Galligan et al. Reference Galligan S.J., Stuart, McManus and Stoffel2022). Altieri (Reference Altieri2008) 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 Reference Kidwell2013; Kidwell and Tomašových Reference Kidwell and Tomašových2013; Harnik et al. Reference Harnik, Torstenson and Williams2017). 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 Reference Hasler1969; Rabalais and Turner Reference Rabalais and Turner2001; Osterman et al. Reference Osterman, Poore, Swarzenski and Turner2005; Breitburg et al. Reference Breitburg, Levin, Oschlies, Grégoire, Chavez, Conley, Garçon, Gilbert, Gutiérrez and Isensee2018). Agricultural and urban runoff delivers excess nutrients to coastal ecosystems like those of the northern Gulf of Mexico (Hasler Reference Hasler1969). 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 Reference Darnell2015). 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. Reference Steckbauer, Duarte, Carstensen, Vaquer-Sunyer and Conley2011; Breitburg et al. Reference Breitburg, Levin, Oschlies, Grégoire, Chavez, Conley, Garçon, Gilbert, Gutiérrez and Isensee2018). 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 Reference Kidwell2007; Harnik et al. Reference Harnik, Torstenson and Williams2017; Tomašových et al. Reference Tomašových, Albano, Fuksi, Gallmetzer, Haselmair, Kowalewski, Nawrot, Nerlović, Scarponi and Zuschin2020). Considering the significant increase in nutrient inputs to the northern Gulf of Mexico beginning in the 1950s (Turner and Rabalais Reference Turner and Rabalais1994; Brunner et al. Reference Brunner, Beall, Bentley and Furukawa2006), 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.
Acknowledgments
For assistance in the field, we thank the staff of the Dauphin Island Sea Lab, Florida State University Coastal and Marine Laboratory, and the Louisiana Universities Marine Consortium, in particular the vessel captains, crew, and marine technical support staff that we have worked with over the years. We are grateful for the contributions of many undergraduate research students at Colgate University and Franklin & Marshall College who have assisted with sample collection in the Gulf of Mexico and processing in the Paleo Lab. We thank J. Stoller for collecting the shell fragmentation data included in this study. For discussion at various points, we thank S. Finnegan and S. Wang. We thank E. Saupe, J. Huntley, A. Klompmaker, and S. Mondal for their constructive reviews of this article. The work described in this paper was supported by a grant from the National Science Foundation (NSF EAR-2041667 to P.G.H.) and an early-career research fellowship from the Gulf Research Program of the National Academies of Sciences, Engineering, and Medicine (GRP 2000008410 to P.G.H.).
Competing Interests
The authors declare no competing interests.
Data Availability Statement
Data available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.n5tb2rc1h.
Open access
