The influence of environmental parameters on spatial variation in zoobenthic density and stable isotopes (δ13C, δ15N, and δ34S) within a large lake

The use of baselines in stable isotope studies to interpret food web structure is essential, but baseline isotope values are often assumed to be spatially homogeneous, even in large aquatic ecosystems. To test this assumption in large lakes, we quantified spatial variation in δ13C, δ15N (deposit‐feeding Oligochaeta and filter‐feeding Dreissena spp.), and δ34S (Dreissena spp. only) and density in Lake Erie between 2014 and 2016. Lake Erie's three distinct basins differ in size, bathymetry, and nutrient loading, making it an excellent system for exploring spatial variation in stable isotopes of baseline organisms. Dreissena spp. densities were highest in the western and lowest in the seasonally hypoxic central basin, while Oligochaeta densities were relatively consistent throughout Lake Erie. Values of δ13C, δ15N, and δ34S exhibited distinct spatial trends that were not related to population densities but followed the west to east direction of water flow within the lake. For both taxa, δ13C was lower in the deeper, oligotrophic east basin than the shallow, mesotrophic west basin, and δ15N and δ34S increased from west to east. Spatial patterns of low δ34S in Dreissena spp. in the western and central basins were likely related to hypoxia, whereas patterns of δ15N in both taxa were probably related to the greater influence of agricultural land uses in the western basin. Spatial trends of stable isotopes in large lake zoobenthos are driven by complex interactions of environmental gradients, which could introduce bias in evaluations of trophic structures within aquatic ecosystems that use stable isotopes.

Anthropogenic stressors are changing the fabric of aquatic ecosystems around the world, and there is an increasing need to make predictions about how they will change, particularly under different climate scenarios (Schmitz and Barton 2014;Pecl et al. 2017).To this end, scientists need to understand the processes and mechanisms that drive aquatic ecosystem structure (Schmitz and Barton 2014;Pecl et al. 2017).However, spatial variability in environmental factors, including nutrient concentrations, depth, light, and others, affect the distribution and abundance of organisms (Jackson et al. 2001;Nye et al. 2009), generating a layer of complexity in understanding aquatic ecosystem processes and functions.Spatial heterogeneity, while not a novel concept, is difficult to quantify and account for in ecosystem studies and is not usually adequately expressed in food web studies, particularly within lower trophic levels (Woodland et al. 2012).
Stable isotopes present an opportunity to understand energy and nutrient flow within aquatic ecosystems (Peterson and Fry 1987;Boecklen et al. 2011;Glibert et al. 2019) and offer an avenue through which spatial heterogeneity can be quantified (Finlay and Kendall 2007;Radabaugh et al. 2013).This is due to the predictable way in which stable isotopes move throughout ecosystem compartments (Fry 2007).Nitrogen (δ 15 N) and carbon (δ 13 C) stable isotopes are the two elements most frequently analyzed to assess trophic relationships in freshwater ecosystems, but isotopes of oxygen (δ 18 O), hydrogen (δ 2 H), and sulfur (δ 34 S) are emerging as useful tracers in ecological research as well (Croisetière et al. 2009;Layman et al. 2012;Ofukany et al. 2014).Significant spatial heterogeneity in stable isotopes has been documented within aquatic ecosystems of various geographic sizes (Graham et al. 2010;Guzzo et al. 2011), and accounting for spatial variability in analyses is necessary, especially when species forage in multiple habitats that have differing stable isotopes (e.g., littoral δ 13 C vs. pelagic δ 13 C; Hobson 1999Hobson , 2005)).
Isoscapes, or spatially explicit maps of isotope values in a resource (e.g., seston, zoobenthos), are a unique and powerful tool to examine spatial trends in ecosystem processes (Bowen 2010;Cheesman and Cernusak 2016;Glibert et al. 2019).Since spatial heterogeneity in isotope values is driven by geochemical nutrient cycling within a system, isoscapes provide a coherent representation of spatially distributed processes occurring in ecosystems (e.g., geological, biological, etc.;Bowen 2010;Kurle and McWhorter 2017), and can be used as isotopic baselines to help interpret trophic position (δ 15 N), and quantify resource use (δ 13 C, δ 34 S) of organisms occupying higher trophic levels (Matthews and Mazumder 2003;Woodland et al. 2012).The use of isoscapes as isotopic baselines helps address the challenges represented by spatial and temporal heterogeneity in lower trophic level organisms typically used as baselines (Matthews and Mazumder 2005;Smyntek et al. 2012;Woodland et al. 2012), since isoscapes are able to integrate this variation and use it to help explain spatial patterns observed in trophic structure and food web dynamics (Kurle and McWhorter 2017).However, isoscapes must be derived from large datasets of the spatially explicit isotopic composition of organisms or other materials within an ecosystem, which often constrains their use due to limited data availability (Glibert et al. 2019).Isoscapes have, thus far, mostly been incorporated into terrestrial and marine research (Bowen 2010;Glibert et al. 2019), but limited work has been done in freshwater systems (Ofukany et al. 2014;Winter et al. 2021).
The relative abundance (or density) of an organism used to create an isoscape is important because rare or uncommon species may not contribute meaningfully to the diet of higher-level consumers and/or be difficult to obtain sufficient samples for spatial analysis.Also, small-bodied organisms (e.g., zooplankton, Oligochaeta) must occur at high enough frequencies to provide sufficient sample mass for stable isotope analysis and be widely distributed enough to interpret spatial and temporal variability.Population density/biomass of an aquatic invertebrate population can be related to the resources available to support it (i.e., carrying capacity), and the presence of environmental conditions it can tolerate (Calizza et al. 2017).Consequently, one could postulate that aquatic invertebrate population densities are related to the water quality within an ecosystem (e.g., dissolved oxygen, temperature, specific conductance, etc.), and will vary according to the spatial variation of these characteristics (Quinn and Hickey 1990;Chambers et al. 2006).
In freshwater systems, isoscapes can be difficult to establish and are rarely detailed enough to quantify the degree of spatial variation within habitats.Background spatial variation in stable isotopes within freshwater ecosystems is often due to spatial heterogeneity in primary productivity, redox conditions (e.g., normoxia vs. anoxia), specific conductance, water temperature, nutrient sources (e.g., sewage treatment plant input, tributary inflow), and surrounding land use (e.g., agricultural vs. urban).For example, δ 13 C values in freshwater ecosystems can be influenced by adjacent land use since both agricultural and urban run-off can lower δ 13 C (Van Dover et al. 1992;Forsberg et al. 1993;DeBruyn and Rasmussen 2002).Sulfur (δ 34 S) isotope values are influenced by redox conditions and are typically lower in hypoxic conditions than normoxic ones (Fry 1986;Croisetière et al. 2009).Values of δ 15 N can be elevated in proximity to sewage treatment plants (Savage and Elmgren 2004) or urban areas (Harvey and Kitchell 2000), and lower in areas of agricultural run-off (Diebel and Vander Zanden 2009), and are therefore also sensitive to nearby land use (Peterson et al. 2007).
While the spatial variation of stable isotopes in lower trophic levels (primary producers, primary consumers) is known to occur and has been shown to impact isotope values of upper trophic levels (secondary and tertiary consumers; Guzzo et al. 2011), less research has been done to investigate reasons for this variation.Here, we investigate the spatial variation of organism density, δ 13 C, δ 34 S, and δ 15 N of pelagic (Dreissena spp.) and benthic (Oligochaeta spp.) zoobenthos within Lake Erie using specimens collected during a comprehensive benthic survey conducted between 2014 and 2016 and relate these to biotic and abiotic parameters.We also quantified population density and determined if this influenced stable isotope values.We expected to measure lower densities of Dreissena spp.and higher densities of Oligochaeta in the central basin of Lake Erie, where there is a history of hypoxia (Burns et al. 2005).For δ 13 C, we anticipated higher values in the western basin than the central and eastern basins due to its shallow depth, and high productivity (i.e., extensive littoral area).Values of δ 34 S were expected to be lowest in regions subject to extended periods of hypoxia due to the metabolism of sulfate by sulfur-reducing bacteria, specifically within the central basin and western basin, where anoxic conditions are common (Olapade 2018).Finally, δ 15 N was expected to increase on a west to east gradient due to greater urban and agricultural influences on the lake in the western basin.We quantified these stable isotope trends to identify spatially constrained baselines that could contribute to better characterizing Lake Erie food webs through stable isotopes and contribute to improving the accuracy of baseline stable isotopes in freshwater ecosystems in general.

Study organisms and system
Lake Erie, the smallest of the Laurentian Great Lakes (hereafter Great Lakes) by volume, is also the productive lake in the system and can be divided into three bathymetrically delineated basins, each with distinct environmental characteristics (Mortimer 1987;Reutter 2019).Most water enters from two main tributaries (Detroit and Maumee rivers) into the west basin, travels the length of the lake, and exits via the Niagara River in the east.Gradients of decreasing nutrient concentrations and increasing depth are evident from west to east (Mortimer 1987;Watson et al. 2016).
Dreissena spp.are an abundant bivalve that invaded the Great Lakes in the late 1980s (Hebert et al. 1989) and have since been used as pelagic isotopic baselines within these systems (Griffiths et al. 1991;Yuille et al. 2012;Karatayev and Burlakova 2022).As filter feeders, Dreissena spp.typically feed on phytoplankton, but they are flexible in their diets and consume a significant amount of particulate detritus (Garton et al. 2005;Jaeger 2006;Kumar et al. 2016); which can mask a "true" pelagic δ 13 C signature (Heuvel et al. 2021).Benthic/ littoral baselines within the Great Lakes are characterized by examining common benthic or littoral invertebrate species such as Oligochaeta worms.Oligochaeta within the Great Lakes are characterized by two main groups within the family Naididae: Tubificinae ("tubificids"), which are deposit feeders (Wavre and Brinkhurst 1971;Brinkhurst et al. 1972), and Naidinae ("naidids"), which are sediment dwellers but consume a wide variety of diet items depending on the species (McElhone 1980;Bowker et al. 1983;Smith and Kaster 1986).

Sample collections and processing
Between 2014 and 2016, samples of benthic invertebrates were collected from 106 locations throughout Lake Erie in support of the US EPA and Environment and Climate Change Canada-sponsored Cooperative Science and Monitoring Initiative (Richardson et al. 2012; Lake Erie Partnership Working Group 2018).Sampling sites were selected across Lake Erie using a probabilistically based stratified random design to estimate lakewide Dreissena spp.and benthic invertebrate population abundance and biomass.Approximately equal sampling effort was apportioned to each of the three basins and to the lake margin vs. offshore habitats.The number of samples allocated among offshore substrate types and depths within each basin was weighted by the area of depth classes, the relative area of hard vs. soft substrate classes, and the relative variability in the density of dreissenids estimated from previous surveys.Sampling was conducted through a collaborative arrangement with Environment and Climate Change Canada, providing cruise time on Canadian Coast Guard Ships Limnos and Shark.
For each site, sample substrate was recorded based on the average size of particles (e.g., clay, silt, or sand), as this can have an impact on zoobenthos density and production (Bonsdorff et al. 1995).Triplicate samples of substrate were collected using a Petite (15 Â 15 cm) or standard (25 Â 25 cm) Ponar grab, rinsed in a sieve bucket to remove silt and clay, and preserved in an ethanol-formalin solution (15 : 6 : 1 parts water : 95% ethanol : 100% formalin).Specific conductance, dissolved oxygen concentration, and redox potential were measured in situ 50 cm above the substrate with a Hydrolab multisonde meter.Average annual water temperature of the water column from 0 to 20 m was estimated from the Great Lakes Aquatic Habitat Framework (Wang et al. 2015).Geographic coordinates and water depth at each site were provided from ship logs at the conclusion of each cruise.
In the laboratory, benthic samples were elutriated to separate organic from inorganic fractions.The organic matter sample portion was rinsed through a nested sieve series (4.0, 1.0, 0.5, and 0.25 mm) to separate the material into size fractions.Each fraction was sorted beneath a dissecting microscope at Â6.4 and/or Â12 magnification.Sample fractions were subsampled (typically 25% or 50%) if they contained large amounts of organic material.All invertebrates were retained, identified to the lowest practical taxonomic level, and archived in 70% ethanol.
Dreissena spp.and oligochaete densities were calculated as the number of individuals/sampler area (i.e., Petite Ponar Grab = 0.152 m Â 0.152 m; Standard Ponar Grab = 0.229 m Â 0.229 m) where number of individuals was obtained from the number of specimens counted within an entire sample (or inferred from subsample estimates).Finally, geometric means of densities were calculated for all replicates/sites.Only Dreissena spp.retained on the 4-and 1-mm sieves were enumerated and measured for stable isotope analysis because smaller dreissenids make up a negligible proportion of the biomass of a sample.Shells were hand-measured to the nearest 0.1 mm using digital calipers.
Up to 20 oligochaete worms from each size fraction (or 80 individuals total) were slide-mounted for identification using CMC 10 aqueous mounting medium (Masters Company).The remainder of the specimens in a sample were used for stable isotope analysis.To minimize the potential for sizerelated bias in identification or isotopic analysis, the worms from a size fraction were emptied into a Petri plate and arranged by length from smallest to largest.The 20 worms to be designated for identification from a sample of N worms were selected by slide-mounting every nth specimen when n = 1/(20/N).For stable isotope analysis, oligochaetes were further grouped into two size fractions to increase the sample mass available for stable isotope analysis ("large" = 4.0-and 1.0-mm size fractions; and "small" = 0.5-and 0.25-mm size fractions), on the basis that some taxa of Oligochaeta are larger or smaller than others.The relative proportions of Tubificinae and Naidinae comprising a stable isotope sample were later inferred from the identification of the specimens in the slide-mounted subset.

Stable isotope analysis
Dreissena spp.and Oligochaeta samples preserved in ethanol were placed in Petri dishes to allow the ethanol to evaporate.Shells were removed from Dreissena spp.samples to prevent a bias in δ 13 C from the valves' CaCO 3 (Jaschinski et al. 2008).Each sample for stable isotope analysis was a composite of multiple individuals from a single Ponar replicate for each site since individual masses were often less than the minimum (400 μg for δ 13 C and δ 15 N; 5000 μg for δ 34 S) required to be analyzed on the mass spectrometer.Samples were frozen at À20 C for 24 h and then lyophilized for 48 h to ensure samples were anhydrous prior to processing for stable isotope analysis.
Dreissena spp.tissues were homogenized into a fine powder using dissection scissors, and a 400-600 μg aliquot was measured out for δ 13 C and δ 15 N analysis on a Delta V Advantage Continuous Flow Mass Spectrometer (Thermo Fisher Scientific) coupled to a 4010 Elemental Combustion System without further processing.Aliquots used for δ 34 S analysis weighed 5000-5500 μg and 300-600 μg of vanadium pentoxide was added to each capsule as an accelerant before being analyzed on a Delta V Plus Continuous Flow Mass Spectrometer (Thermo Fisher Scientific) coupled to a 4010 Elemental Combustion System (Costech Instruments).Instrument precision, as assessed by four internal lab standards (NIST1577c, USGS 42, NIST 8555, and tilapia muscle) run in every 10 samples, was ≤ 0.3 (SD) for δ 34 S. Instrument accuracy for δ 34 S was +0.3‰ (mean difference from certified value) for NIST 8554 and NIST 8555, and À0.3 for NIST 8529 (n = 30 for all).
The mass of most of the Oligochaeta samples was so small (< 400 μg) that δ 34 S could not be determined for them.Accordingly, tests were conducted on the mass spectrometer used for δ 13 C and δ 15 N to determine the smallest sample mass that could be weighed and still provide an accurate measurement of δ 13 C and δ 15 N.For this, a certified standard (bovine liver, NIST 1577c) and a large Oligochaeta sample (mass > 3000 μg) were weighed out in triplicate to masses of 100, 200, 300, 400, 500, and 600 μg and analyzed on a Delta V Advantage Continuous Flow Mass Spectrometer (Thermo Fisher Scientific) coupled to a 4010 Elemental Combustion System to create a calibration curve that showed machine accuracy at different weights.All standards for each test weight were within AE 0.1‰ (mean difference of expected vs. measured values) for δ 13 C, and AE 0.3‰ for δ 15 N of certified values, so it was concluded that oligochaete samples as small as 100-200 μg would yield consistently accurate measurements of δ 13 C and δ 15 N. Instrument precision, assessed by the standard deviation of standards (NIST1577c, USGS 40, urea, and tilapia muscle) run in every 10-12 samples, was ≤ 0.1‰ (SD) for δ 13 C, and ≤ 0.2‰ for δ 15 N. Instrument accuracy of δ 13 C and δ 15 N checked throughout the period of time samples were run was +0.1‰ (mean difference from certified value) for δ 13 C and À0.1‰ for δ 15 N, and was assessed by replicates of NIST standards 8542 (δ 13 C only), 8547 (δ 15 N only), 8573, and 8574 for both (n = 50 for all).
To analyze spatial patterns of stable isotopes in Dreissena spp.and Oligochaeta, sites that had data from multiple Ponar grabs (i.e., two or more samples for Dreissena spp. or Oligochaeta were associated with a site) were averaged to minimize bias from resampling in geospatial statistics.In addition, within-site variation of stable isotopes of δ 13 C, δ 15 N, and δ 34 S was between AE 0.2‰ and 0.3‰ SE, which is consistent with that expected of sample replicates, and insignificant in comparison to the spatial patterns observed.This resulted in a total of 43 and 80 sites for each taxon, respectively, used to document spatial distribution, and permit spatial interpolation, and geospatial regression.Spatial distribution of the density and stable isotope data were assessed with Moran's I to determine the initial dispersion before spatial modeling was conducted (Haining 2001).The distribution of Dreissena spp.and Oligochaeta density and stable isotope values was mapped to represent spatial variation using an inverse-weighting function in ArcMap (version 10.8.1).
Multiple linear regressions were performed on the density and stable isotope data for both Oligochaeta and Dreissena spp. to assess the spatial variability of each.Variables included in multiple linear regression analysis for zoobenthos densities and stable isotopes were water quality data (annual average water temperature, DO, and specific conductance), site depth, site distance to the mouth of the Detroit River, site sediment quality percent Tubificinae (proportion of Tubificinae in sample; Oligochaeta models only), and shell length (Dreissena spp.models only).For the most part, these variables were measured and recorded in the field except for the distance to the Detroit River, which was calculated using the "Point Distance" tool in ArcMap.Sediment quality was input into models as a dummy variable, which ranked the classifications of sediment type by particle size (i.e., clay = 0, silt = 1, sand = 2).Densities of dreissenids and oligochaetes were log-transformed for statistical analysis since the density distributions were skewed.A model produced by the multiple linear regression was considered valid and "passing" if it met the criteria of an R 2 > 0.2, all variables contributed to the model significantly (p < 0.05), residuals were normally distributed (Jarque-Bera test; p > 0.05), residuals were not spatially autocorrelated (Moran's I p > 0.05), and VIF was less than 5.0.The best-fit model for each species was identified as the model that met all the above criteria and had the highest R 2 with the lowest Akaike's information criterion.Ordinary least squares linear models were fitted to the best model for each isotope and taxon.Multiple linear regression analyses and ordinary least squares regressions were performed in ArcMap (10.8.1) using the exploratory regressions and ordinary least squares tools.

Taxonomic composition of Oligochaeta samples
The relative proportions of Oligochaeta subfamilies varied little between the size fractions collected or among sample locations.Overall, most oligochaetes collected were Tubificinae (> 90% throughout the entire lake), and most of the remaining specimens were either Naidinae (7.6 AE 1.3% in the small size fraction and 0.4% in the large size fraction) or individuals belonging to other taxa within the family (e.g., Lumbriculidae; < 1.1% on average in either size fraction; Table 2).The relative taxonomic composition of Oligochaeta did not vary much among basins at the subfamily level, and many samples from the small and large size fractions in all basins consisted mainly of Tubificinae (> 89% on average; Table 2).Furthermore, values of δ 13 C and δ 15 N did not differ significantly between size fractions in any basin (Kruskal-Wallis tests; χ 2 ≤ 1.6, df = 1, p ≥ 0.2 for all; Table 2).

Spatial distribution
There was little evidence of spatial autocorrelation of sample densities of Dreissena spp.and Oligochaeta across Lake Erie based on Moran's I analysis (Moran's I; Oligochaeta: Moran's I = 0.0, z = À0.1, p = 0.9; Dreissena spp.: Moran's I = 0.0, z = À0.2, p = 0.9).Spatial interpolation of the population density of both taxa showed that the density of dreissenids was highest around the mouths of the Maumee and Sandusky rivers and lowest in the central basin (Fig. 1A.Oligochaeta were more uniformly distributed throughout the lake, with several regions of high and low densities within each basin (Fig. 1B).
The distribution of δ 13 C, δ 15 N, and δ 34 S (Dreissena spp.only) values was spatially autocorrelated for both Oligochaeta and Dreissena sp.based on Moran's I analysis (Moran's I; Oligochaeta: Moran's I δ13C = 0.4, z = 9.1, p < 0.001; Moran's I δ15N = 0.5, z = 9.8, p < 0.001; Dreissena spp.: Moran's I δ13C = 0.2, z = 5.0, p < 0.001; Moran's I δ15N = 0.4, z = 8.1, p ≤ 0.001; Moran's I δ34S = 0.2, z = 5.1, p < 0.001).The distribution of all three isotopes for both taxa exhibited spatial heterogeneity and gradients among basins.Values of δ 13 C were generally higher for both Oligochaeta and Dreissena spp. in the western basin and lower in more easterly locations of Lake Erie (Fig. 2).In addition, δ 13 C values were higher in samples collected from shallower water (e.g., western basin, Long Point Bay, and the Pennsylvania Ridge; Fig. 2) than elsewhere.Values of δ 15 N were lowest in the western basin for both Oligochaeta and Dreissena sp. and were greater in samples collected in more eastern locations of the lake (Fig. 3).Similar to δ 15 N, δ 34 S values, which were only measured in Dreissena spp., were greater at eastern than at western locations.The lowest values occurred in the central basin and in the plume of the Maumee River in the western basin (Fig. 4).

Geospatial regression
No models were identified as significant for any combination of variables that could predict the density of either Oligochaeta or Dreissena spp.In both taxa, none of the variables were significant in the multiple regression analysis, and consistent with this, the R 2 values for models were low.This indicates that the environmental variables used in the analysis were likely not the main variables affecting the distribution of zoobenthos density throughout Lake Erie and the use of time-averaged values might be better suited to estimating spatial variability in their densities.
No models passed the criteria set in the multiple linear regression analysis for values of δ 13 C or δ 15 N in Oligochaeta.The lack of any models identified as passing in the multiple linear regression for δ 13 C, and δ 15 N values suggests that we lacked important variables explaining the spatial variation of stable isotopes of Oligochaeta.Dreissena sp.δ 13 C was best explained by a model that included one variable, shell length (R 2 = 0.28, F 1,41 = 24.1,p < 0.001; Table 3).Three other models also passed the criteria set for multiple regression analysis of δ 13 C values in Dreissena sp., which included one 1-variable model and two 2-variable models (Supporting Information Table S1).A total of five models passed the criteria set for exploratory regression for Dreissena sp.δ 15 N values.The best-fit model for Dreissena sp.δ 15 N explained 70% (R 2 = 0.68, F 3,38 = 41.4,p < 0.001; Table 3) of the variation observed and included two variables: shell length (partial R 2 = 0.33), and water depth (partial R 2 = 0.35; Table 3).The other four models that met the criteria set for the multiple regression analysis all included shell length as a variable in a combination of two-and threevariable regression models (Supporting Information Table S1).Exploratory regression of Dreissena spp.δ 34 S found no passing models (Supporting Information Table S1).The absence of any other passing models was mostly due to variables not loading significantly onto tested models in the multiple regression analysis.

Discussion
This study highlights that spatial variation of stable isotopes can be significant in large lakes, which could impact the interpretation of stable isotopes in upper trophic levels in food web studies.Values of δ 13 C and δ 15 N had ranges of 6-12‰ for both Oligochaeta and Dreissena spp.across Lake Erie, and ranges of 5-6‰ within geographically distinct basins, highlighting the spatial heterogeneity of stable isotope values in large lake ecosystems.Values of δ 34 S, measured only in Dreissena spp., spanned close to 9‰ across the lake and 4-6‰ Table 3. Summary of the best-fit models for Dreissena spp.δ 13 C and δ 15 N values in Lake Erie collected between 2014 and 2016 were identified through multiple regression analysis.Best-fit models were selected based on the combination of factors that met the following criteria: high R 2 , low AIC c , and residuals were normally distributed (Jarque-Bera statistic p > 0.05) and not spatially autocorrelated (Moran's I p > 0.05).No models passed the criteria for δ 34 S values.within basins.Diet tissue discrimination factors are usually assumed to be $ 0.5-1.0‰for δ 13 C, $ 3.4‰ for δ 15 N, and $ 0‰ for δ 34 S (Peterson and Howarth 1987;Post 2002), spatial variation on a magnitude of 6-12‰ in organisms used as baselines to interpret trophic structure (e.g., Dreissena spp.) introduces significant bias into trophic ecology studies if not properly addressed in statistical analyses.For example, Guzzo et al. (2011) found that mismatches between the location of capture for white perch and yellow perch and the corresponding baseline organism (Dreissena spp.) in the western basin of Lake Erie biased their trophic position calculations.
Given that many fish species in Lake Erie, such as the commercially and recreationally valuable walleye (Sander vitreus), make regular and large annual movements, food web studies in large aquatic ecosystems need to assess and incorporate spatial variation in stable isotope analyses.Significant spatial variation of δ 13 C, δ 15 N, and δ 34 S was observed in both Oligochaeta and Dreissena spp.collected from across Lake Erie.Although population densities showed no spatial analysis trends, there were regions of higher and lower densities throughout Lake Erie.Since zoobenthos communities are likely more influenced by long-term trends in water quality, it is possible zoobenthos densities would be related to variables such as duration of hypoxia.Trends of δ 13 C decreased, and δ 15 N and δ 34 S (Dreissena spp.only) increased moving east through the lake in both taxa, which is the same direction as the flow of water and gradients of increasing depth and decreasing nutrient concentrations in the system.Episodic high-nutrient run-off from agricultural land in the Maumee River basin has significant impacts on ecosystem dynamics, which contrasts with low nutrient concentrations from the Detroit River and likely influences the spatial patterns observed in all isotopes and taxa (Bridgeman et al. 2012).While no significant models were identified for values of δ 13 C or δ 15 N in Oligochaeta, they were correlated with values of both isotopes in Dreissenidae, indicating ecosystem processes probably drove some of the spatial heterogeneity observed.Specific research that examines stable isotope distribution in freshwater systems with water and lake characteristics is limited, and more are needed.More specific analysis, for example, time-integrated sampling, might help clarify these results.Although local inputs, that is, high agricultural inputs in the western basin, and water flow were likely important determinants of the spatial distribution of stable isotopes in benthic invertebrates in Lake Erie, results also highlight the importance of environmental conditions in large aquatic systems.Based on these results, food web studies using stable isotopes in such systems need to consider not only inputs and water flows but also variations in environmental parameters such as water quality indices (e.g., DO), productivity (nutrient availability) and climate (e.g., temperature regimes, seasonal influxes of nutrients).
Oligochaeta populations throughout Lake Erie were predominantly Tubificinae spp.and taxonomic composition had no effect on δ 13 C or δ 15 N in our study.This is likely because Tubificinae spp.were the dominant taxa in both size fractions analyzed for stable isotopes and Naidinae spp.did not make up more than 8.1% of a sample on average, suggesting that the isotope values reflect mostly that of Tubificinae spp.The taxonomic composition of Oligochaeta samples is consistent with previous benthic community analyses of Lake Erie's western and eastern basins, which observed an Oligochaeta community dominated by Tubificinae spp.(Carr and Hiltunen 1965;Dermott and Kerec 1997;Karatayev et al. 2022).
Ontogeny likely had an influence on isotope values in Dreissena spp., as indicated by the inclusion of shell length in the best-fit models for δ 13 C and δ 15 N.While a relationship between shell length and δ 13 C and δ 15 N values has not previously been observed in Dreissena spp., it has been observed in Unionid mussels (Yasuno et al. 2014), which was associated with ontogenetic changes in diet.Here, the relationship between δ 13 C and δ 15 N values with length could be associated with ontogeny as well, since larger Dreissena spp.could consume larger particles and potentially zooplankton.
The spatial distribution of density of both Dreissena spp.and Oligochaeta were random according to Moran's I, and were not significant in stable isotope models, demonstrating population density did not influence stable isotope values and, by convention, their feeding ecology.Dreissena spp.densities were highest around the mouths of the Maumee and Sandusky rivers, two tributaries that discharge nutrient-rich water into Lake Erie.Oligochaeta were more evenly distributed throughout the lake, which is likely related to their dependence on sediment-derived resources rather than plankton.The lowest densities of Oligochaeta and Dreissena spp.occurred in areas of the central basin.For dreissenids, this is probably due to the occurrence of episodic hypoxia in that basin in most years, especially during the late summer and early fall (Charlton et al. 1993;Burns et al. 2005).Dreissena spp.densities were much lower at most locations (< 100 individuals m À2 ) in the central basin, excluding two sample sites, which had densities of 324 and 515 individuals m À2 .Dreissena spp.are sensitive to hypoxia, which has been shown to limit their population densities in the central basin of Lake Erie (Karatayev et al. 2018), whereas Oligochaeta are tolerant to hypoxia (Karatayev et al. 2022).This was not supported by multiple regression analysis, which did not find any significant models to explain the spatial distribution of either taxon.However, the data used for the multiple regression analysis was from in situ measurements taken at the time of sampling and not based on time-averaged values, which could confound our results.Indeed, while many benthic invertebrate communities within Lake Erie have recovered from poor water quality conditions in the 1950s and 1960s, those present in the central basin still resemble the communities present during the 1960s and are dominated by hypoxia-tolerant taxa like oligochaetes (Karatayev et al. 2022).Overall, the densities of Dreissena spp.measured in this study were similar to those measured in 2000 (Patterson et al. 2005), 2004, and 2010 (J. Ciborowski et al., University of Windsor, unpubl.), and 2019 (Karatayev et al. 2022).Patterson et al. (2005) reported that lakewide average densities were around 600 individuals m À2 .Oligochaete densities are also similar to past reports in 2004 (J.Ciborowski et al., University of Windsor, unpubl.), and 2019 (Karatayev et al. 2022).
Dreissenids had lower δ 13 C and δ 15 N than oligochaetes based on estimates of a lakewide mean value for each basin.These differences are indicative of their feeding ecology and the structure of pelagic and benthic food web pathways.Dreissena spp.are pelagic filter feeders (e.g., phytoplankton; Garton et al. 2005), whereas oligochaetes are benthic deposit-or epipelic-feeding organisms that depend on detritus and bacteria (Tubificinae; Brinkhurst et al. 1972) or epibenthic biofilm and algae (Naidinae; Brinkhurst et al. 1972;Timm 2012).Pelagic food webs are mainly supported by the production of primary producers such as phytoplankton, whereas profundal food webs are supported by microbial decomposition of detritus (Berglund et al. 2007).Bacteria-based food chains such as those formed in benthic ecosystems are typically longer than their pelagic counterparts (Berglund et al. 2007), which explains the higher δ 15 N observed in oligochaetes than dreissenids.In addition, δ 13 C tends to be lower in organisms dependent on planktonic algae than those dependent on benthic algae (Croisetière et al. 2009), which explains the lower δ 13 C observed in dreissenids.Similar patterns of higher δ 13 C and δ 15 N in Oligochaeta (mean AE SD, δ 13 C = À18.4AE 1.3‰, δ 15 N = 8.2 AE 0.6‰) than Dreissena spp.(δ 13 C = À25.8AE 0.4‰, δ 15 N = 7.3 AE 0.3‰) were observed in Lake Simcoe (Ozersky et al. 2012).Values of between 6.0 and 9.8‰ for δ 15 N were reported for Lake Erie Dreissena spp.and around 11‰ in sediment-associated benthic organisms (e.g., Chironomidae;Garton et al. 2005;Campbell et al. 2009;Guzzo et al. 2011), consistent with the values reported in this study.Campbell et al. (2009) and Guzzo et al. (2011) both measured δ 13 C in Dreissena spp.(approximately À22-28‰) that were consistent with values in this study.No records of Oligochaeta δ 13 C were found for Lake Erie, but values measured in this study were consistent with sediment-associated benthic invertebrates (e.g., chironomids) reported in previous research (Campbell et al. 2009).
Among-basin and depth-related spatial patterns of δ 13 C and δ 15 N in Lake Erie were similar in both Oligochaeta and Dreissena spp., indicating that spatial variation observed in both isotopes likely reflect ecosystem-level processes (i.e., nutrient inputs, cycling, and water flow) rather than differences in either species feeding ecology.This was corroborated by the results of multiple regression analyses, which identified depth as one of the top variables correlating with the spatial variation observed in δ 15 N in Dreissenidae.Research reported in the Great Lakes and elsewhere has also noted that depth appears to be an important correlate of trends in δ 13 C and δ 15 N (Johannsson et al. 2001;Brock et al. 2006;Hu et al. 2012).Trends of lower δ 13 C in macrophytes and sediments found at greater depth have been observed in a combination of marine and freshwater systems (Brock et al. 2006;Hu et al. 2012), which is not consistent with the patterns observed in δ 13 C in Dreissena spp.and Oligochaeta in this study.As a result, spatial patterns of δ 13 C values in zoobenthos in this study are likely an effect of the net autotrophic/high productivity status of the western basin transitioning to the net heterotrophic/low productivity status of the eastern basin.Johannsson et al. (2001) reported that Lake Ontario zooplankton collected from deeper water had greater δ 15 N than zooplankton collected from shallow water, although data limitations restricted the interpretation of the relationship, and match the patterns observed in Lake Erie Dreissena spp. in this study.In addition, increasing δ 15 N values with depth were observed in zoobenthos in Lake Superior and were thought to reflect the decomposition of sinking phytoplankton (Sierszen et al. 2006).In other aquatic ecosystems, trends of higher δ 15 N in deeper water have been associated with the depletion of nearshore 15 N by N 2 fixing prokaryotes, and the influence of depth on the remineralization of organic matter (Brock et al. 2006;Puccinelli et al. 2018).
In the western basin, spatial heterogeneity of stable isotopes likely reflects influxes of nutrients from Lake Erie's two main tributaries, the Detroit and Maumee Rivers, which supply 95% of water to the lake annually (Bolsenga and Herdendorf 1993;Niu and Xia 2021).Distinct plumes of δ 13 C were observed in sediments of Alaskan coastal lakes adjacent to stream discharge locations (Brock et al. 2006), similar to the patterns of δ 13 C and δ 15 N in Oligochaeta and Dreissena spp.observed at the mouths of the Detroit and Maumee rivers in this study.The Detroit River is the main water inflow into Lake Erie ($ 90% total water input; mean river flux 5000 m 3 s À1 ) and carries water predominantly from the upper Great Lakes (Superior, Michigan, and Huron), which is colder, clearer, and more oligotrophic than Lake Erie (Binding et al. 2012;Niu et al. 2018;Niu and Xia 2021).This influx of lower nutrient water is likely responsible for lower δ 15 N values observed in both taxa on the north side of the western basin where the Detroit River plume discharges.The warmer Maumee River has the largest drainage area of any Great Lakes tributary (mean river flux $ 200 m 3 s À1 ), and its watershed is dominated by agricultural and industrial land uses, which release a large quantity of nutrients into the western basin of Lake Erie (Cipoletti et al. 2019).In this study, lower δ 13 C in Oligochaeta, and higher δ 15 N in both Dreissena spp.and Oligochaeta around the mouth of the Maumee River could be explained by the agricultural and urban influence on nutrientstable isotopes being deposited there.This is also consistent with observations of higher δ 15 N values in zoobenthos in wetlands at the mouths of rivers in the Great Lakes (Peterson et al. 2007).In addition, patterns of spatial variation of both δ 13 C and δ 15 N in Lake Erie's western basin observed in yellow perch (Perca flavescens), white perch (Morone americana), zooplankton, and dreissenids were consistent with the findings of this study (Guzzo et al. 2011).
Relatively few Dreissena spp.samples from the central basin of Lake Erie were available for analysis, limiting our ability to interpret spatial trends of δ 34 S in this region.Dreissena spp.are relatively intolerant of hypoxic conditions (Karatayev et al. 2018), which is reflected in their low densities in the central basin reflects.Values of δ 34 S generally increased from west to east, in the direction of water flow within the system, which is probably related to the large benthic influence in the western basin.No adequately fit model was identified for δ 34 S, likely due to a combination of the absence of key variables, lack of samples in the central basin, and the fact that the water quality measurements used were single time points and not time-integrated averages.Hypoxia in the western and central basins might be responsible for the lowest values of δ 34 S observed, considering that δ 34 S is lower when sulfur-reducing bacteria are present (Finlay and Kendall 2007).The bathymetry of the central basin of Lake Erie renders it vulnerable to hypoxia since it is deep enough to stratify, but the hypolimnion is typically thin ($ 5-7 m thick) in comparison to the epilimnion (15-20 m thick; Charlton et al. 1993;Burns et al. 2005;Watson et al. 2016).As a result, the volume of water in the hypolimnion is usually not sufficient for the oxygen demand of the benthic processes (Burns et al. 2005;Watson et al. 2016).
Despite the sparsity of data points for Dreissena spp.(n = 5) in the central basin, the pattern of δ 34 S matches what is known about the average extent of the hypoxic zone there (Burns et al. 2005).Episodic hypoxia is common in the south-western part of the western basin due to periodic thermal stratification on calm days during summer months (Bridgeman et al. 2006), which is also consistent with small areas of lower δ 34 S in Dreissena spp. in that region.The understanding of δ 34 S dynamics in freshwater systems is less developed than in marine or estuarine ecosystems, but our results suggest that its distribution in lower trophic levels in freshwater aquatic food webs is partially related to oxygen availability and hypoxia events.Additional research investigating the effect of hypoxia on δ 34 S in lower trophic levels used as isotopic baselines is needed in Lake Erie and elsewhere.
In conclusion, large spatial heterogeneity was observed in δ 13 C, δ 15 N, and δ 34 S (Dreissena spp.only) of Lake Erie zoobenthos among Lake Erie's three basins, possibly reflecting the influence of nutrient inputs from agricultural, urban, and industrial land use from the Maumee River watershed.Zoobenthos population density varied randomly throughout the lake and was unrelated to the spatial patterns of stable isotopes throughout the lake.Distinct plumes of δ 13 C and δ 15 N were observed in both taxa collected near the mouth of the two main tributaries of Lake Erie, and isotope data were further correlated with water quality parameters and water depth, suggesting that environmental conditions play a role in influencing isotope ratios in freshwater ecosystems.Dynamics of δ 34 S in Dreissena spp.highlighted the large variability of stable isotopes in natural systems and corroborated recent findings highlighting the use of δ 34 S in trophic ecology studies as a complementary tracer to the typically used δ 13 C and δ 15 N. Finally, care should be taken when designing trophic ecology studies using isotopes to make sure that the spatial extent of the study area is appropriate to the spatial use of the targeted organism(s).

Fig. 2 .
Fig. 2. Spatial distribution of δ 13 C in (A) Dreissena spp.and (B) Oligochaeta within Lake Erie between 2014 and 2016 as determined by inverse-distance weighted interpolation.Darker values indicate higher (less negative) values of δ 13 C, and sites where sufficiently large samples were collected are indicated by dark red triangles and black dots for Dreissena spp.and Oligochaeta, respectively.

Fig. 3 .
Fig. 3. Spatial of δ 15 N in (A) Dreissena spp.and (B) Oligochaeta within Lake Erie between 2014 and 2016 as determined by inverse-distance weighted interpolation.Lighter values indicate lower values of δ 15 N, and sites where sufficiently large samples were collected are indicated by dark red triangles and black dots for Dreissena spp.and Oligochaeta, respectively.

Fig. 4 .
Fig. 4. Spatial of δ 34 S in Dreissena spp.within Lake Erie between 2014 and 2016 as determined by inverse-distance weighted interpolation.Lighter values indicate lower values of δ 34 S, and sites where sufficiently large samples were collected are indicated by dark red triangles.

Table 1 .
Density, δ 13 C, δ 15 N, and δ 34 S of Oligochaeta and Dreissena spp.collected in Lake Erie between 2014 and 2016.Density is reported as the number of individuals m À2 , all stable isotope values are in ‰ and Oligochaeta does not have any δ 34 S data.Bolded lines are lakewide averages.

Table 2 .
Summary of taxonomic composition (mean AE SE) of Tubificinae and Naidinae in Oligochaeta collected in different size frac-