Bidirectional energy subsidies for fish in river mouths of a large lake as revealed by stable isotopes and fatty acids

Freshwater river mouths may facilitate cross‐habitat resource use and bidirectional subsidization through (a) active cross‐habitat movement by consumers like fish and (b) water and material movement both downstream, from tributary to nearshore lake, and in the opposite direction due to backflow processes. Despite potential importance of freshwater river mouths, relative contributions of lentic and tributary energy sources to fishes along these ecotones remain understudied, in part because measuring cross‐habitat subsidization is not straightforward. We assessed cross‐habitat subsidization of three fish consumers round goby (Neogobius malanostomus), yellow perch (Perca flavescens), and alewife (Alosa psuedoharengus) in three river mouths of Lake Michigan, USA. Specifically, we (a) measured stable isotope ratios (δ13C, δ15N, δ2H, δ18O) in ambient water, bivalve shells, potential invertebrate prey, and fish soft tissue, and we assessed fish nutritional status based on essential fatty acid composition, and (b) estimated Bayesian ellipses and mixing models including either two (δ13C‐δ15N) or three (δ13C‐δ15N‐δ2H) isotope ratios. Results revealed evidence of bidirectional habitat subsidies for all fishes, but energy subsidies varied by species and across river mouths; likely related to differences in species‐specific behavior and hydrologic differences among river mouths. Relative to nearshore oligotrophic Lake Michigan, fish in productive tributary environments contained relatively low essential fatty acid content, suggesting a possible trade‐off between prey quantity and quality across habitats. This study advances our understanding of resource connectivity in lake ecosystems and the combined use of multiple stable isotopes (including mixing models) and fatty acids to document cross‐habitat subsidies along freshwater ecotones.

versa (Richardson and Sato 2015), and where delivery of nutrients and prey may contribute to coastal productivity (Carassou et al. 2011).
Although studies of habitat utilization and resource subsidies in freshwater river mouths are rare, some studies have attempted to use stable isotope analyses to describe consumers' resource use in these transition zones (Dufour et al. 2008;Senegal et al. 2020;O'Reilly et al. 2023).In fish, stable isotopes in soft tissues are integrative indices of resource use, reflecting food consumed and habitats occupied over the past 2-6 months, depending upon tissue turnover (Weidel et al. 2011).Stable carbon (δ 13 C) and nitrogen (δ 15 N) isotope ratios are among the most commonly used markers in foraging ecology studies including those of estuarine transition zones (Post 2002;Hoffman et al. 2010).Both δ 13 C and δ 15 N of primary producers can vary widely within and among systems depending on the form of assimilated carbon and nitrogen (Post 2002).For example, agricultural nitrogen inputs may influence the 15 N abundance in receiving water bodies, and variation in δ 15 N from allochthonous inputs may be reflected in consumer soft tissues (Larson et al. 2013).While most trophic studies in freshwater river mouths have relied upon δ 13 C and δ 15 N, some studies have incorporated ratios of additional isotopes (e.g., Senegal et al. 2020).Isotope ratios of hydrogen (δ 2 H) and oxygen (δ 18 O) in a consumer can differentiate between allochthonous and autochthonous resources and index fish habitat use across hydrological gradients (Doucett et al. 2007;Soto et al. 2013;Vander Zanden et al. 2016).The simultaneous consideration of isotope ratios of four elements represents a potentially effective approach to elucidate relative lentic and lotic resource contributions in food webs and may provide additional resource discrimination over two element measurements in the vicinity of freshwater river mouths (Vander Zanden et al. 2016;Senegal et al. 2020).
Fatty acids play a role in a variety of physiological processes over the ontogeny of organisms and may serve as additional food web tracers.Long-chained polyunsaturated fatty acids (PUFAs), such as alpha-linoleic acid (ALA, 18:3n-3), arachidonic acid (ARA, 20:4n-6), eicosapentaenoic acid (EPA, 20:5n-3), and docosahexaenoic acid (DHA, 22:6n-3) are particularly important diet components for the growth and reproduction of fish (Tocher 2010).Because fish cannot synthesize several essential fatty acids, the composition of these fatty acid compounds in fish tissues will partially reflect their prey (Müller-Navarra et al. 2004;Jardine et al. 2015).For example, lake eutrophication and brownification alter phytoplankton structure, potentially reducing the availability of DHA and long-chain PUFA, which may be reflected in the fatty acid composition of fish (Taipale et al. 2016).Simultaneously, such responses may have important implications for the nutritional condition of individual consumers (Jardine et al. 2015;Taipale et al. 2016).
Herein, we describe a study to elucidate resource contributions supporting three species of fish in river mouth systems connected to Lake Michigan, USA.In the Laurentian Great Lakes, including Lake Michigan, river mouth habitats are known to be relatively productive (Höök et al. 2007;Schoen et al. 2016), provide diverse foraging opportunities, and serve as vital spawning and nursery grounds for various fish species (Dufour et al. 2005;Höök et al. 2007;Janetski and Ruetz 2015).Lake Michigan is an oligotrophic system that has recently experienced further oligotrophication due to decreased nutrient loading and the filtering activities of invasive Dreissena spp.mussels (Barbiero et al. 2019).It is plausible that due to decreased lake productivity, the relative ecological importance of tributaries and river mouths in supporting the nearshore Lake Michigan fish assemblage has increased (Höök et al. 2007;O'Reilly et al. 2023).However, river mouths connected to Lake Michigan differ broadly in terms of their discharge characteristics and morphology, which may influence the relative bidirectional exchange of material and resource subsidization.Specifically, we studied (a) invasive round goby (Neogobius melanostomus) that primarily feeds demersally in nearshore habitats and tends to occupy small territorial ranges (Ray and Corkum 2001), (b) yellow perch (Perca flavescens), a native species in Lake Michigan that occupies diverse habitats across nearshore and river mouth habitats and feeds on a broad spectrum of prey from zooplankton to benthic macroinvertebrates and fish (Happel et al. 2015;Foley et al. 2017), and (c) introduced alewife (Alosa pseudoharengus), a primarily planktivorous, anadromous fish species in its native range, which retains some seasonal migratory behavior in Lake Michigan (Höök et al. 2007).
Our primary objective was to quantify the proportional subsidization of resources ("Lake" vs. "River") in fish collected in the vicinity of three river mouth systems (encompassing the lower river reach to the local zone of potential fluvial influence in nearshore Lake Michigan).The three river mouths differed in key attributes such as watershed size, upstream land use, discharge, and morphology.We hypothesized (a) that differences in river mouths attributes will promote variation in habitat subsidization, and (b) that the degree of cross-habitat resource subsidization will vary among the three study species, reflecting differences in feeding habits and migration behaviors.We relied on stable isotope analyses of δ 13 C and δ 15 N, which directly reflect assimilated diets in tissues, as well as measurements of δ 2 H and δ 18 O that additionally reflect the composition of local ambient water.We evaluated isotopic overlap among fish collected in different habitats using isotopic ellipses and quantified resource use and energy subsides using dual-isotope and tripleisotope Bayesian mixing models, with δ 13 C and δ 15 N, and including δ 2 H, respectively.Further, we quantified fatty acids of study species and focused on essential fatty acids (e.g., DHA and EPA) as a supplemental index of resource reliance and individual nutritional condition.

Study site
This study was conducted in three southeast Lake Michigan river mouth areas (Trail Creek, St.Joseph River, and Muskegon River; Fig. 1).We define a river mouth as the zone encompassing a tributary (upstream to the extent that lake influences, such as backflow from the recipient lake, are apparent) and downstream as far as the river plume is potentially evident in the recipient lake (Larson et al. 2013).Although Lake Michigan is generally oligotrophic in relation to its tributaries, the rivers selected for this study differed in drainage areas, hydrology, geomorphology, and watershed land use, resulting in differences in productivity relative to Lake Michigan.In brief, Trail Creek, the smallest of the rivers (mean flowrate of 2 m 3 s À1 ; watershed area of 181 km 2 ), flows through an urban and industrially developed watershed (19% agricultural land) into a small dredged harbor, with evidence of river backflow in summer and fall (Grimm 2015;Essig 2016).In contrast, St.Joseph River flows through a large, mostly agricultural watershed (watershed area of 7540 km 2 ; ≈ 70% agricultural land; DeGraves 2006) and has a higher discharge (mean flowrate 97 m 3 s À1 ; Essig 2016).St.Joseph River flows directly into Lake Michigan and can be characterized by a relatively large river plume during spring and summer (Essig 2016;Jameel et al. 2018).Muskegon River (348 km long, watershed of 4382 km 2 ) drains a largely forested area (49% forested; Ray et al. 2010) and flows through a drowned river mouth lake, Muskegon Lake (up to 24 m depth), and an artificially reinforced Muskegon Channel to Lake Michigan (Marko et al. 2013).Muskegon Lake is a highly productive habitat (Höök et al. 2007), where riverine and watershed-derived materials settle and are processed before reaching Lake Michigan (Marko et al. 2013).
All three river mouth sites were sampled monthly during spring, summer, and fall from May to October in 2011 and 2012.Two to three locations were typically sampled along the 10 m depth contour in Lake Michigan, encompassing a range of 0.2-2.0km from the mouth of each river.Due to inconsistent catches of fish, all sampling sites from nearshore Lake Michigan, as well as all sites within tributaries, were pooled together to represent two environmentally distinct habitats: "Lake" and "River" (Fig. 1a,b).In Muskegon, the spatial complexity and additional sampling within the Muskegon river mouth facilitated comparisons among four unique river mouth habitats: Muskegon River, Muskegon Lake, Muskegon Channel, and nearshore Lake Michigan (Fig. 1c).

Fish and resource sampling
Fish species: At each river mouth location, round goby, yellow perch, and alewife were collected with adaptive sampling approaches-gillnets, boat electrofishing (Model SR16H; Smith-Root Inc.), beach seine nets, and hook and line sampling to target round gobies in Trail Creek and Muskegon Lake (for methodology details, see Stein 2018).
Age 0 individuals were excluded because young fish are likely to differ in habitat occupancy and resource use from older fish (Roswell et al. 2014).Length thresholds were used to identify fish Age 1 and older.For round goby, individuals > 45 and > 55 mm total length (TL) collected before or after August, respectively, were considered ≥ Age 1+ (Corkum et al. 2004).For yellow perch, individuals > 80 and > 90 TL collected before or after August, respectively, were considered ≥ Age 1+ (Fitzgerald et al. 2004).Alewife > 80 mm TL collected through August were considered ≥ Age 1+ (Höök et al. 2007).Upon capture, fish were transported on the dry ice to the laboratory.
For a subset of fish, $ 0.2 g of white muscle tissue was removed from the dorsal mid-section and stored at À80 C for subsequent fatty acid analysis.The remaining whole fish soft tissues, after removing fish guts, were homogenized (Schielke and Post 2010) and stored at À20 C for stable isotope analysis.Stein (2018) additionally analyzed a subset of yellow perch to evaluate isotopic consistencies between homogenized whole fish and muscle tissue and found small variations in stable isotope ratios and qualitatively similar patterns among individuals.
Ambient water, shells, and prey: Ambient water was collected weekly $ 0.3 m below the surface in 1-dram glass vials, sealed, and stored at 4 C.In addition, we analyzed δ 13 C and δ 18 O ratios in whole crushed sessile bivalve shells (Dreissena spp.) to described local water dynamics integrated over a longer time period (McConnaughey and Gillikin 2008).Shells were collected in petite PONARs (see below), and the δ 13 C and δ 18 O of these shells should primarily reflect the isotopic composition of dissolved inorganic carbon and oxygen in ambient water.
Potential prey samples, including bulk zooplankton (henceforth, "seston") and benthic invertebrates, were collected across all sites.Seston was collected with a 64-μm vertical plankton net (2.0 m long with 0.5 m diameter), towed from 0.5 m above the bottom to the surface with a speed of 0.2 m s À1 .Benthic invertebrates were collected with a petite PONAR grab (0.023 m 2 , 500 μm), and each grab's contents were sorted to separate chironomid larvae as a benthic source, and other benthic invertebrates as a "mixed" benthic source (including amphipods, oligochaetes, and gastropods).Given that they feed on seston and represent a pelagic source, soft tissues from Dreissena spp.were not included as benthic isotope sources.The inclusion of potential prey items as isotopic sources was also informed by stomach contents of collected round goby, yellow perch, and alewife as summarized by Stein (2018).All prey samples were stored at À20 C until further analyses.

Stable isotope and fatty acid analyses
Stable isotopes: Analyses of δ 13 C, δ 15 N, δ 2 H, and δ 18 O were conducted at the Stable Isotope Ratio Facility for Environmental Research, University of Utah, USA.All biological samples were dried at 60 C for 24 h and lipid washed prior to stable isotope analyses, air-dried, and homogenized samples to a powder (Parrish 1999).All isotopic values of samples are reported using δ notation as follows: δX = ([R sample /R standard ] -1), in units of per mil (‰), where R is a ratio of 13 C : 12 C, 15 N : 14 N, Delta Plus or Delta V mass spectrometer, and calibrated to the Vienna Pee Dee Belemnite (VPDB) and atmospheric N 2 (AIR) reference scales using two-point linear calibration based on two organic standards (UU-CN-1: δ 13 C = 23.9‰,δ 15 N = 49.6‰;UU-CN-2: À28.2‰, À4.6‰; average precision for both δ 13 C and δ 15 N based on replicate samples was AE 0.1‰).Samples for δ 2 H and δ 18 O ratios were analyzed in duplicate on a Thermochemical Elemental Analyzer coupled to a Delta V mass spectrometer.Samples were stored alongside two in-house keratin standards (KHS: non-exchangeable ) for a period of 2 weeks to allow exchangeable hydrogen in all materials to equilibrate with hydrogen in the ambient laboratory water vapor.Non-exchangeable Η and O isotope data were calibrated to the Vienna Standard Mean Ocean Water-Standard Light Antarctic Precipitation (VSMOW-SLAP) scale, using a two-point linear calibration based on measured data for the KHS and CBS standards (average precision based on replicate samples was 4.1‰ for δ 2 H and 0.4‰ for δ 18 O).
Ambient water samples were analyzed using Cavity ring-down spectroscopy (Picarro L2130-i) according to Jameel et al. (2016).Accuracy and precision were checked throughout the period of analysis using internal laboratory reference waters; precision of the instruments used was within AE 0.20‰ and AE 0.04‰ for δ 2 H and δ 18 O, respectively (Jameel et al. 2016).
Whole, crushed bivalve shells were converted to CO 2 in an orthophosphoric acid bath, and the purified gas was reacted and analyzed on a Thermo Finnigan Gas Bench II coupled with a Delta V Isotope Ratio Mass Spectrometer and calibrated to the VPDB reference scale using two-point linear calibration to measured values for the Carrara Marble standard (δ 13 C = 2.1‰, δ 18 O = À1.8‰) and Lithium Carbonate (LSVEC) standard (δ 13 C = À46.6‰,δ 18 O = À26.7‰;average precision for both δ 13 C and δ 18 O based on replicate analyses of a third calcite standard was AE 0.1 ‰).
Fatty acids: To quantify fatty acids, white muscle tissues samples were dried for 48 h, weighed, and then lipids were extracted with a modified method using chloroform : methanol (2:1, v/v) with 0.01% butylated hydroxytoluene (Parrish 1999;Budge et al. 2006).Muscle lipids were transesterified to fatty acid methyl esters (FAMEs) with methanol and a 3% concentrated sulfuric acid solution (Budge et al. 2006).Overall, 43 individual FAMEs were detected using a Shimadzu gas chromatograph model GC-2010 with flame ionization detector and autoinjector with EZStart software version 7.4 (Shimadzu Corp.;Feiner et al. 2018).The proportions of individual fatty acids were expressed as a mass percentage of the total fatty acid content.

Data analysis
Biplots and isotopic standard ellipses: Within each river mouth, we considered intraspecific resource overlap across paired habitats (nearshore Lake Michigan vs. tributary; differences among the four Muskegon habitats) by fitting ellipses that comprise a central trend of isotopic data and encompass approximately 95% of data for δ 13 C vs. δ 15 N and δ 2 H vs. δ 18 O biplots (Jackson et al. 2011).For all species, the isotopic standard ellipse area was corrected for a small sample size (SEA C ), and proportional ellipse overlaps among habitats were estimated using the function "maxLikOverlap" in the SIBER package (version 2.1.7;Jackson et al. 2011) in R version 4.0.3(R Core Team 2020).The percentage overlap of two habitats was calculated as the isotopic area of overlap divided by the total isotopic area covered by the two ellipses (non-overlap SEA C area).Muskegon tributary overlaps were estimated between each habitat pair including across nearshore Lake Michigan, Muskegon Channel, Muskegon Lake, and Muskegon River habitats.Finally, potential differences among isotopic niche widths (‰ 2 ) were obtained from Bayesian standard ellipse area (SEA B ) including credibility intervals for statistical comparisons (Supporting Information Table S1; see Jackson et al. 2011).Ellipse areas and metrics were calculated using default model settings in SIBER (see Jackson et al. 2011).
Mixing model including δ 13 C and δ 15 N: Resource subsidies for fish across habitats were estimated using Bayesian mixing models in MixSIAR package (version 3.1.12;Stock et al. 2018) in R version 4.0.3(R Core Team 2020).MixSIAR models were modified from the standard rjags package to operate using the runjags package in R which facilitated parallel processing (i.e., four cores) and extension of Markov Chain Monte Carlo (MCMC) chain length (Turschak et al. 2022).Models for each fish species collected in each habitat were structured with one random factor-"Habitat" (e.g., "Lake" or "River").Trophic discrimination factors (TDFs) (mean AE SD) of 0.4‰ AE 1.3‰ for δ 13 C and 3.4‰ AE 1.0‰ for δ 15 N were provided as model inputs (Post 2002).Models were run with four parallel MCMC chains of length 50,000, burn-in length 20,000, and chains were thinned by 10, with uninformed priors.Model convergence was evaluated using trace plots, and Gelman-Rubin Diagnostic not exceeding 1.05 for more than 5% of parameter estimates as described in Turschak et al. (2022).Model endmembers from nearshore Lake Michigan and river mouth habitats were weighted by assumed equal mean contribution of benthic and pelagic prey (50 : 50).In the Muskegon tributary, we estimated Lake Michigan vs. tributary subsidies for fish in each of the four habitats only including two end-members in the mixing models, that is, nearshore Lake Michigan and Muskegon River sources.Note that resulting model estimates of habitat subsidies in fishes may be conservative if seston and invertebrate habitat source end-members were partially subsidized by reciprocal habitats.For Bayesian mixing model applications, assumptions, error structures, and best practices see Stock et al. (2018).
Mixing models including δ 13 C, δ 15 N, and δ 2 H: We included a third isotopic tracer, δ 2 H, to potentially improve the model accuracy and to assess the effect of δ 2 H tracer on the estimated source proportions.δ 2 H integrates dietary and environmental sources, with contributions of approximately 33% ambient (environmental) water and 67% diet for fish (Soto et al. 2013;Vander Zanden et al. 2016).MixSIAR models with three isotopes (δ 13 C-δ 15 N-δ 2 H) had the same model parameters as the dual-isotope models (δ 13 C-δ 15 N).TDF of δ 2 H was assumed to equal zero due to the reported ambient water accumulation effect or "compounding effect" at each trophic level (Solomon et al. 2011;Soto et al. 2013;Coulter et al. 2017).In this case, we can use the mass-balance model for fish tissue δ 2 H and we assumed a contribution of ambient water on fish protein isotope values of 33% (ω; Solomon et al. 2009;Soto et al. 2013;Vander Zanden et al. 2016).Accordingly, we corrected δ 2 H values for ambient water effect δ 2 H ambient-w in fish (same as δ 2 H diet ) as derived from the equation by Vander Zanden et al. ( 2016): where τ is a trophic difference between resource and consumer (= 1).Ultimately, we executed two versions of δ 13 C-δ 15 N-δ 2 H mixing models.One version included δ 2 H values accounting for ambient water contribution (δ 2 H ambient-w ) using sitespecific ambient water δ 2 H values measured in this study, and the second version only included originally measured δ 2 H fish values (not accounting for ambient water effect [ω]).
Fatty acid statistical analyses: We analyzed indices of essential fatty acids that fish cannot typically synthesize in freshwater ecosystems (Tocher 2010;Happel et al. 2015;Jardine et al. 2015).For each fish species, we evaluated cross-habitat differences between total PUFAs and relative contents of selected essential fatty acids (DHA, EPA, ALA, and ARA) using the Welch test of one-way ANOVA.Prior to statistical tests, we transformed (arcsine of square root) all fatty acid relative content data to approximate normality.Furthermore, we applied the post hoc Dunnett-T3 test for pairwise multiple comparisons to determine which groups were significantly different.This statistical analysis was conducted using IBM SPSS software version 28.0.Evaluation of statistical significance in all tests was set to p < 0.05.

Biplots and isotopic standard ellipse analyses
Ambient water in all river mouth sites showed a trend of increasing δ 2 H and δ 18 O values from tributaries to Lake Michigan, and data generally grouped into three isotopically distinct regions: "Nearshore Lake Michigan," "Rivers," and "Muskegon tributary" (Fig. 2).In nearshore Lake Michigan δ 18 O mean values ranged from À6.2‰ to À5.8‰ and δ 2 H mean values ranged from À46.0‰ to À44.7‰.Water isotope ratios in St.Joseph River and Trail Creek were distinct from Lake Michigan, with Trail Creek being more enriched and hence more similar to Lake Michigan than St.Joseph River.The lowest mean water δ 2 H and δ 18 O values were observed in Muskegon tributary sites, particularly in Muskegon River (À8.6‰ δ 18 O and À60.0‰ δ 2 H).However, even in the Muskegon tributary, a gradual increase in isotopic values was observed from Muskegon River to Muskegon Lake (À7.9‰ δ 18 O and À57.0‰ δ 2 H) to Muskegon Channel (À7.6‰ δ 18 O and À54.2‰ δ 2 H) to nearshore Lake Michigan (Fig. 2).
The δ 13 C and δ 18 O values of Dreissena spp.bivalve shells followed a pattern similar to ambient water with isotopic enrichment in nearshore Lake Michigan, relative to tributaries (Fig. 3).The largest within habitat variability of δ 13 C and δ 18 O among individual bivalve shells was observed in Trail Creek and Muskegon Channel.In contrast, bivalve shell δ 13 C and δ 18 O values were the least variable in St.Joseph River.
There were clear isotopic differences in fish between nearshore Lake Michigan and tributary habitats, with differences generally tracking isotopic gradients of water and potential prey (Fig. 5).In biplots δ 13 C-δ 15 N ellipses indicated habitat overlaps which were largely supported by complementary δ 2 H-δ 18 O ellipses (Fig. 5; Supporting Information Table S1).In comparing nearshore Lake Michigan vs. tributary habitats, round goby exhibited relatively low habitat overlaps (Fig. 5; Supporting Information Table S1).In contrast, yellow perch and alewife displayed greater cross-habitat overlaps.For example, using δ 13 C-δ 15 N ellipses in Trail Creek, no between-habitat overlap was detected for round goby, while yellow perch showed 8%, and alewife up to 47% overlap between tributary and nearshore Lake Michigan ellipses (Fig. 5a).Furthermore, we noticed that δ 2 Hδ 18 O ellipses tended to show slightly higher cross-habitat overlaps in comparison to estimated overlaps based on δ 13 C-δ 15 N ellipses.Relative to St.Joseph and Trail Creek, the Muskegon river mouth showed more complex overlap patterns across the four habitats (Fig. 5c; Supporting Information Table S1).Finally, we also compared the isotopic niche widths using SEA B between sites and species and found that in most cases alewife, followed by yellow perch, had larger SEA B niche widths than round gobies (Supporting Information Table S1).
Mixing models with δ 13 C, δ 15 N, and including δ 2 H Cross-habitat subsidization of fish species was reflected by MixSIAR model estimates of resource use varied among species and river mouths (Fig. 6; Supporting Information Table S2).In general, round goby and yellow perch had a strong reliance on local resources (except for round goby in nearshore Lake Michigan near St. Joseph), whereas alewife was estimated to have received the most subsidization across habitats (Fig. 6).For example, in nearshore Lake Michigan near Trail Creek, round goby, and yellow perch incorporated 95-98% of the local lake resource, whereas alewives were nearly equally supported by Lake Michigan and Trail Creek river resources (55% of "Michigan" resource, model δ 13 C-δ 15 N, Fig. 6a; Supporting Information Table S2).Interestingly, while round goby ellipses did not indicate habitat overlap for individuals collected near the St.Joseph river mouth (Fig. 5b; Supporting Information Table S1), the mixing model estimated a large contribution of riverine subsidies to these nearshore Lake Michigan individuals (44% of "River" resource, model δ 13 C-δ 15 N; Fig. 6b; Supporting Information Table S2).In contrast, round goby collected in the tributary relied entirely on local resources in the St.Joseph River (Fig. 6b).Finally, model results in Muskegon implied broad resource subsidization among multiple habitats (Fig. 6c).
We compared MixSIAR models using δ 13 C and δ 15 N with two additional modeling approaches employing a third isotope tracer, δ 2 H.Both three-isotope models demonstrated that including a third tracer, δ 2 H or δ 2 H ambient-w , only slightly affects source proportion estimates across habitats and species (relative to two isotope ratio model estimates; Fig. 6; Supporting Information Table S2).

Fatty acids
Total PUFA levels (% PUFA relative to all fatty acids detected) among fish species were similar (round goby 32- 48%, yellow perch 34-55%, alewife 36-56% of total fatty acids; Fig. 7).Within the fatty acid pool, DHA displayed particularly high variability, and DHA content in round goby was much lower (5.2%AE 1.3%; Fig. 7a) than in yellow perch (18.0%AE 1.8%; Fig. 7b) and alewife (15.1% AE 2.2%; Fig. 7; Supporting Information Table S3).Fatty acid levels of all three species varied significantly among habitats, but alewife showed less consistent patterns across habitats than the other  Lake Michigan (in a legend as "Michigan") and Trail Creek (a), St.Joseph (b), and Muskegon (c) river mouths (in a legend as "River," "Channel," and "Lake" habitats).Each data point represents a single fish that was collected between May andOctober (2011-2012), and the number of individuals sampled varies depending on the habitat.Ellipses incorporate 95% of the data and are corrected for a small sample size (SEA C ). Color coding and habitats of each tributary are shown in a legend.Text indicates the relative isotope overlap area (Overlap %) comparing "Michigan" and "River" habitats in all cases.Note that overlap area may be somewhat underestimated in some habitats due to the limited sample size (*).For additional Overlap % estimates across "Channel" and "Lake" habitats in Muskegon tributary, see Supporting Information Table S1.
two species.In Muskegon, PUFA and DHA levels of all fishes showed a clear spatial difference across habitats, and for some fish species significant changes were also detected in EPA, ARA, and ALA levels.DHA increased from river to nearshore Lake Michigan for round goby (Welch's ANOVA, F 1,27 = 30.27,p < 0.05; Fig. 7a) and yellow perch (Welch's ANOVA, F 3,30 = 6.51, p < 0.05; Fig. 7b), but this gradient was not consistent for alewife (Welch's ANOVA, F 2,19 = 15.85,p < 0.05; Fig. 7c).In the other two river mouths, only yellow perch in Trail Creek exhibited a significant difference in DHA levels between tributary and nearshore habitats (Welch's ANOVA, F 1,9 = 9.2, p < 0.05; Fig. 7b).Limited successful collection of Fig. 6.Dual (δ 13 C-δ 15 N) vs. triple-isotope mixing model (δ 13 C-δ 15 N-δ 2 H ambient-w ) resource proportion estimates for round goby, yellow perch, and alewife across nearshore Lake Michigan (in a legend as "Michigan") and Trail Creek (a), St.Joseph (b), and Muskegon (c) river mouths (in a legend as "River," "Channel," "Lake" habitats).X-axis labels represent double-isotope model estimates, except with the superscript "H" they represent triple-isotope model estimates with δ 2 H ambient-w .Black and white lines represent estimated mean posterior densities and bars indicate 95% credible intervals.Bars illustrating subsidies in fish species across Lake Michigan and tributary habitats from two sources-"Michigan" (color-filled bars in line with the colors of ellipses of "Michigan" habitat in Fig. 5) vs. "River" (transparent bars).n.d., no data available.For individual model mean, median, and credible interval estimates, see Supporting Information Table S2.For site abbreviations, see the legend in Fig. 2.
round goby (Fig. 7a) in Lake Michigan and of yellow perch (Fig. 7b) in the St.Joseph River precluded some analyses of fatty acid composition between lake and river habitats.

Discussion
Study findings revealed evidence of bidirectional resource subsidies for three fish species between tributaries and nearshore Lake Michigan habitats.Consistent with hypotheses, the magnitude and direction of cross-habitat subsidies varied among river mouths and species.Moreover, essential fatty acids revealed that fish occupying relatively oligotrophic nearshore Lake Michigan tended to express greater nutritional status relative to fish in productive tributary environments.Collectively, this study suggested complex resource subsidization in these freshwater ecotones, where the absence of pronounced salinity gradients may facilitate high cross-habitat exchanges.

Isotopic patterns in ambient water, shells, and prey
The δ 13 C, δ 15 Ν, δ 2 H, and δ 18 O ratios of ambient water, shells, and potential prey collectively demonstrated broad differences across nearshore Lake Michigan and tributary habitats, as well as distinct patterns within river mouth systems.Our study measured higher δ 13 C values in nearshore prey Fig. 7. Mean composition (mass % of the total fatty acid content AE SE) of the sum of PUFAs and commonly used essential fatty acid biomarkers (DHA, EPA, ARA, ALA) in round goby (a), yellow perch (b), and alewife (c) in nearshore Lake Michigan (as "Michigan") and Trail Creek, St. Joseph, and Muskegon river mouths (in a legend as "River," "Channel," and "Lake" habitats).Significant differences between each habitat (Welch's ANOVA and post hoc Dunnett-T3 tests, p < 0.05) are indicated with letters (A-B) and connected with gray lines to indicate spatially significant fatty acid differences across habitats.n.d., no data available.Mean contents of 18 detected fatty acids (saturated, monounsaturated, and polyunsaturated) are summarized in Supporting Information Table S3.For site abbreviations, see the legend in Fig. 2.
from Lake Michigan and lower δ 13 C values in prey from tributaries (Larson et al. 2012;Marko et al. 2013).This trend of low-to-high δ 13 C values in prey along the tributary to lake continuum aligns with previous studies on δ 13 C patterns in lake river mouths (Larson et al. 2012;Marko et al. 2013;Turschak et al. 2014;Senegal et al. 2020).In lakes, autochthonous production by primary producers and less reliance on terrestrial carbon can contribute to greater δ 13 C values (e.g., Finlay et al. 2002).In contrast, tributary allochthonous production will potentially result in relatively low δ 13 C because of different terrestrial carbon sources that primarily enter tributary habitats with low δ 13 C values (Larson et al. 2012(Larson et al. , 2013;;Marko et al. 2013).
We also identified notable δ 15 N habitat differences and isotopic enrichment in prey in some tributary habitats.Enriched δ 15 N values in St.Joseph River prey may be due to substantial nitrogen loading from agriculture (≈ 70% agricultural land cover; DeGraves 2006) and domestic sewage inputs rather than natural biological processes (Peterson et al. 2007).In contrast, Trail Creek with a smaller watershed area and lower proportional agricultural land use (Essig 2016) exhibits lower δ 15 N values in potential prey.Finally, in Muskegon Lake, natural processes, hydrology, and microbial processing may have increased δ 15 N values in prey, influencing isotopic differentiation relative to adjacent upstream and downstream habitats (Marko et al. 2013).Drowned river mouth lakes, like Muskegon Lake, with recurring hypolimnetic hypoxia and high organic matter loading from wetlands may favor microbial processes like denitrification and, as a result, notably elevate δ 15 N values in biota (Ray et al. 2010;Marko et al. 2013).Nonetheless, identifying the exact processes contributing to habitat tributary differences in δ 15 N ratios may be challenging due to additional trophic interactions affecting isotopic composition in prey (Post 2002).
Ambient water isotope patterns presented herein are generally consistent with past studies in nearshore Lake Michigan (e.g., Jasechko et al. 2014;Foley et al. 2017Foley et al. , 2022;;Senegal et al. 2020) and with relatively few available studies that have simultaneously considered tributaries (e.g., Dufour et al. 2005;Jameel et al. 2018;Senegal et al. 2020;O'Reilly et al. 2023).Similar to δ 13 C in prey samples, δ 2 H and δ 18 O values in the ambient water followed a general trend of relatively low δ 2 H and δ 18 O values in tributaries and isotopic enrichment in nearshore Lake Michigan waters (e.g., Foley et al. 2022).The rationale is that the lighter isotopes would be preferentially removed during water evaporation, and this typically has a greater impact on long-residing water bodies like large lakes (Bowen et al. 2012).
River mouths experience seasonal water discharge variability, with particularly high river discharge and nutrient loading in spring and lower discharge and even potential backflow in other seasons (Essig 2016;Jameel et al. 2018).δ 13 C and δ 18 O isotopic variation in bivalve shells indicated an occurrence of backflows from Lake Michigan in Trail Creek and Muskegon Channel.For example, Trail Creek flow is highly variable, with spring discharge into Lake Michigan and episodic wind-driven backflow of nearshore waters and other materials in other seasons (Grimm 2015;Essig 2016).In contrast, St.Joseph River and Muskegon Lake showed less evidence of influence from Lake Michigan.This suggests that the fast-flowing, highdischarge St.Joseph River and sheltered Muskegon Lake are less likely to receive materials or experience water backflows from Lake Michigan.Collectively, the δ 13 C, δ 15 N, δ 2 H, and δ 18 O ratios in water, bivalve shells, and lower trophic level consumers demonstrated a consistent spatial distinction among habitats.These isotopic baselines can be used to characterize habitats and potential resource subsidization in large lakes.Moreover, the findings suggest that analysis of δ 2 H and δ 18 O in complex freshwater environments can help detect riverine system impacts on large lakes and vice versa.

Isotopic patterns in fish
Similar to water and potential prey isotopic patterns, all three fish species reflected isotopic differences between nearshore lake and tributary habitats.Fish in tributaries are expected to have lower δ 2 H, δ 13 C, and δ 18 O values than fish in lake waters, while fish in Lake Michigan would have lower δ 15 N.In addition, we found that fish in river mouth regions are subsidized bidirectionally across habitats, likely due to two mechanisms: (a) transport of water, nutrients, and prey across environments and (b) active movement of fish between habitats.Thus, in line with other studies in nearshore Lake Michigan (Höök et al. 2007;Larson et al. 2013;Happel et al. 2015;Foley et al. 2022), and other Great Lakes (Larson et al. 2016), our findings support that a component of nearshore fish production is subsidized by river mouth resources.Importantly, we also identified a subsidization from oligotrophic Lake Michigan to fish inhabiting productive tributaries, though the extent varied by species and river mouth.
For each fish species we observed dominant habitat utilization strategies.Specifically, round goby seemed to use local resources (showed no overlap of isotopic ellipses; Fig. 5b), while yellow perch and especially alewife, showed more evidence of broad foraging.Such habitat utilization strategies were consistent across river mouth systems and were in line with other Great Lake studies (Edsall 1970;Schoen et al. 2016;Senegal et al. 2020).Based on isotopic ratios in bivalve shells, we expected that subsidization of lake sources would be more common in tributary systems with recurring backflow effects (i.e., Trail Creek and Muskegon Channel), but not in tributary fish from fast-flowing tributary systems like St.Joseph River.However, mixing model outputs of yellow perch and round goby collected in Trail Creek had relatively low lacustrine resource subsidization (< 1% lake "Michigan" source), while alewife showed high lake subsidization in both Trail Creek (39% lake "Michigan" source) and St.Joseph tributary habitat (65% lake "Michigan" source; Supporting Information Table S2).The lack of subsidization from Lake Michigan in Trail Creek fish, except alewife, was somewhat surprising.This may partially reflect that Trail Creek tributary end-members were also partially influenced by Lake Michigan inputs, and thus do not represent pure riverine inputs.Furthermore, alewife exhibit strong migratory strategies between habitats (Edsall 1970) and thus may act to integrate resources (Höök et al. 2007;Dufour et al. 2008).Hence, these findings suggest that fish migration behavior may have contributed to differences in estimated resource subsidization across tributaries and migratory fish like alewife can exploit both riverine and lacustrine resources across diverse habitats.
Muskegon Lake is known to be a highly productive habitat (Höök et al. 2007;Janetski and Ruetz 2015), and in comparison to Trail Creek and St. Joseph, the Muskegon tributary system suggested a complex resource habitat network.Mixing model outputs showed that all three fish species, including less mobile round goby, integrate a variety of resources across Muskegon habitats.Resource exchange indicates that Muskegon habitats are likely interconnected due to the physical characteristics of this system (Essig 2016;O'Reilly et al. 2023) and that fish may move between these habitats for feeding and refuge.In addition, based on isotopic variation in bivalve shells, Muskegon Channel may be an important convergence zone for fish between the productive Muskegon Lake and nearshore Lake Michigan.
Round goby in nearshore Lake Michigan near the St.Joseph River displayed interesting resource reliance.On the one hand, isotopic biplots and ellipses revealed no overlap and between round goby in nearshore Lake Michigan and St.Joseph tributary (Fig. 5b).However, based on isotopic mixing models we observed that round goby in nearshore Lake Michigan near St.Joseph were notably subsidized by tributary resources (46% "River" source; Supporting Information Table S2).We suggest that these two conclusions are consistent, and that round goby relied on a combination of tributary and nearshore Lake Michigan resources.If there was no subsidization, we would expect even more distinct ellipses between these two fish habitat groups.In total, these results are consistent with St.Joseph River having a relatively high river discharge and influx of river sources into Lake Michigan.Yellow perch and alewife in nearshore Lake Michigan near St.Joseph were seemingly also subsidized by riverine resources.On their own, isotopic ellipses may not reveal resource subsidization across habitats, while mixing models allow for quantification of such phenomenon.Collectively, these findings confirm that subsidies vary among tributaries and diverse habitats and reveal the potential for bidirectional exchange of energy resources in river mouth systems.

Including δ 2 H in mixing models
By integrating measurements of three stable isotopes with δ 2 H as a third tracer in Bayesian mixing models, we generated resource contribution estimates very similar to dual δ 13 C-δ 15 N model estimates and confirmed patterns of resource subsidization across habitats.Moreover, our study identified only slight differences in estimated habitat contributions when including original δ 2 H values in mixing models and when accounting for ambient water (δ 2 H ambient-w ).The similarity of results in both three-isotope models may be associated with minimal trophic compounding of ambient water between fish and potential prey (one trophic level) in this study, whereas more complex trophic step consideration in a model may clarify some individual differences of resource subsidization (Solomon et al. 2009(Solomon et al. , 2011;;Heikkinen et al. 2022).In dietary ecology, δ 2 H complicates interpretation as the isotopic composition of consumer tissues reflects a mix of ambient water and dietary isotopic compositions (Soto et al. 2013;Vander Zanden et al. 2016).Trophic compounding and the contribution of ambient water can be associated with trophic feeding levels in fish and is an important consideration for δ 2 H ratios in isotope models (Solomon et al. 2009;Soto et al. 2013;Vander Zanden et al. 2016).As long as this uncertainty can be accounted for (Solomon et al. 2011;Heikkinen et al. 2022), the δ 2 H models alone may be a tool for estimating resource proportions in fish (Doucett et al. 2007).Nevertheless, future studies should simultaneously examine δ 13 C, δ 15 Ν, δ 2 H, and δ 18 O patterns in fish tissues (e.g., Hayden et al. 2015) to explore the potential of using multi-isotope approaches in conjunction with multi-stable isotope mixing models.

Fatty acid patterns in fish
The relative concentrations of different fatty acids varied by species and generally agreed with other Lake Michigan studies (Czesny et al. 2011;Feiner et al. 2018).Importantly, across species, fatty acid concentrations, especially PUFAs, and DHA, varied by habitats.Earlier studies have also identified fatty acid composition-based spatial diet heterogeneity in nearshore Lake Michigan (Czesny et al. 2011;Happel et al. 2015Happel et al. , 2019;;Foley et al. 2017).We observed a higher proportion of PUFA and DHA in nearshore lake fish than in tributary sites, indicating differences in fatty acid availability, uptake and/or retention, and generally suggesting better nutritional conditions in Lake Michigan (Müller-Navarra et al. 2004).In contrast, tributaries that are biologically more productive and serve as an important nursery habitat (Dufour et al. 2005), may be associated with low-quality diets.Allochthonous contributions from land use and watershed modifications may reduce the supply of essential fatty acids to tributary consumers by limiting the input of riparian plant materials and by increasing sedimentation, which can contribute to lower PUFA, DHA, and EPA levels in food webs (Müller-Navarra et al. 2004;Taipale et al. 2016).Moreover, these habitats presumably vary in terms of dominant autochthonous primary producers, including different phytoplankton taxa, which synthesize different concentrations of specific PUFAs (Taipale et al. 2020).Differential fatty acid availability may affect fish growth and survival in tributaries (Tocher 2010).In contrast, some tributary habitats, such as Muskegon Lake, contain greater prey densities (e.g., Höök et al. 2007).Thus, future studies might consider the trade-off between prey quantity and quality across these habitats.

Conclusions
This study provides empirical evidence of bidirectional resource subsidies for fish between productive tributary and nearshore oligotrophic Lake Michigan habitats.The magnitude and direction of subsidies varied across river mouths and fish species.In general, migratory alewife displayed the greatest integration of resources across habitats.While less mobile round goby displayed distinct isotopic profiles among habitats, there was also evidence of cross-habitat resource subsidization by this species.For example, round goby collected in nearshore Lake Michigan near the large St.Joseph River discharge showed clear subsidization by riverine resources.By assessing levels of essential fatty acids (e.g., PUFA and DHA), we found that fish in nearshore oligotrophic Lake Michigan were generally characterized by greater nutritional conditions relative to fish in tributaries.In total, study findings contribute to our understanding of the complex connectivity of habitats and resources within a large lake system, and specifically point to the need for continued study and potential conservation of river mouth transition zones.

2H
Fig. 1.Map of Lake Michigan, USA (light gray) and catchment areas (dark gray) for focal river systems.Maps of river mouth areas are displayed in panels (a) Trail Creek, (b) St. Joseph, and (c) Muskegon.Sampling locations in nearshore Lake Michigan and tributary habitats are indicated by black dots.

Fig. 2 .
Fig. 2. Isotopic composition of δ 18 O and δ 2 H of ambient water in nearshore Lake Michigan and tributary sites.Biplot shows the mean isotopic variation (AE SE) of water samples collected from May throughout October(2011)(2012). Dashed ellipses are intended to demonstrate spatial isotopic variations that group together in distinct regions.Abbreviations for tributary names and sites are listed in the legend.Note that the same legend applies to Figs.3, 4, 6, and 7.

Fig. 3 .
Fig. 3. Isotopic composition of δ 13 C (left) and δ 18 O (right) of bivalve shells (Dreissena spp.) in nearshore Lake Michigan (bottom bars) and tributary sites (top bars).Horizontal bars in a box plot indicate median proportional values.Edges of the boxes represent the approximate 1 st and the 3 rd quartiles.Horizontal error bars extend to the lowest and highest data value inside a range of 1.5 times the interquartile range (1.5 IQR), respectively.Outliers are represented by the dots outside the 1.5 IQR.Sample size (n) is displayed for each site.For site abbreviations, see the legend in Fig. 2.