Phosphorus Dynamics and Availability in the Nearshore of Eastern Lake Erie: Insights From Oxygen Isotope Ratios of Phosphate

Introduction

In the North American Laurentian Great Lakes, blooms of filamentous green algae (mostly Cladophora sp.) were one of the most visible signs of cultural eutrophication in the early 1950s through to the late 1970s, particularly in Lake Erie (Shear and Konasewich, 1975; Higgins et al., 2008). Phosphorus (P) abatement programs under the Great Lakes Water Quality Agreement (GLWQA) resulted in rapid and substantial reductions in external P inputs, that while primarily aimed at reducing extensive phytoplankton blooms, hypoxia, and fish kills, also appeared to reduce Cladophora blooms although confirmational data are limited (Painter and McCabe, 1987). Despite these apparent improvements, reports of expansive blooms of Cladophora increased in frequency in eastern Lake Erie in the mid to late 1990s (Higgins et al., 2005) following P load reductions but also colonization of Lake Erie by exotic zebra (Dreissena polymorpha) and quagga (D. rostriformis bugensis) mussels (Griffiths et al., 1991; Roe and MacIsaac, 1997). It was postulated that the recurrence of these blooms during a period of stable or declining P loads (Dolan and Chapra, 2012) was in part related to structural and functional changes to the ecosystem mediated by large populations of filter feeding dreissenid mussels (Hecky et al., 2004) that increased water clarity (Howell et al., 1996), expanded areas of hard substrate for algal attachment (Haltuch et al., 2000), and enriched the benthic environment through the excretion of metabolic wastes and undigested algae (Arnott and Vanni, 1996; Conroy et al., 2005). At the same time, changing agricultural practices, precipitation patterns, and weather conditions have increased soluble P loadings to western Lake Erie, resulting in a recurrence of cyanobacterial blooms (Michalak et al., 2013). While there is general agreement reductions in P availability are the most practical means with which to reduce the extent and severity of Cladophora blooms (Auer et al., 2010), uncertainty remains as to the relative importance of direct P inputs, P circulating within the lake, and P potentially recycled in the benthic environment by dreissenid mussels (Bootsma et al., 2015). These uncertainties have confounded efforts to develop revised P loading targets for the east basin of Lake Erie to address Cladophora blooms (GLWQA, 2015).

Several recent studies have indicated that the relative abundance of oxygen stable isotopes in phosphate (δ¹⁸O_P) carries information regarding the sources and biological cycling of phosphate (PO₄, hereafter P_i) in aquatic environments (Colman et al., 2005; McLaughlin et al., 2006b, 2013; Elsbury et al., 2009; Joshi et al., 2015; Li et al., 2017). The strong bond between oxygen and phosphorus atoms in phosphate are resistant to oxygen exchange under most environmentally relevant conditions (Paytan and McLaughlin, 2012). While isotopic fractionation effects of abiotic processes such as sorption and desorption appear to be relatively minor (Jaisi et al., 2010, 2011), fractionation during biologically mediated enzymatic reactions can be large and at times, substrate and enzyme dependent (Blake, 2005; Liang and Blake, 2006, 2009). The uptake and internal cycling of P_i by intracellular pyrophosphatases (EC 3.6.1.1) results in complete exchange of all oxygen atoms from P_i with oxygen atoms from water, leading to a rapid temperature-dependent equilibrium between oxygen in P_i and intracellular water (Longinelli and Nuti, 1973; Blake, 2005), which is expected to be equal to that of ambient water (δ¹⁸O_W) (but see Li et al., 2016b). This isotopic equilibration has been observed in many environments to be an indicator of P_i turnover (McLaughlin et al., 2006b, 2013). Hydrolysis of organic phosphorus (OP) by various extra- and intracellular phosphatases and phosphohydrolases (for example, alkaline phosphatase; EC 3.1.3.1, 5' nucleotidase; EC 3.1.3.5, phosphodiesterase EC 3.1.4.1) impart a negative fractionation at the P-O bond site, shifting the δ¹⁸O_P of the product (P_i) further away from equilibrium (Liang and Blake, 2006, 2009). Consequently, knowledge of these reactions and pathways may provide a means to identify sources of P or the importance of specific P cycling pathways governing P_i availability in aquatic ecosystems (Elsbury et al., 2009; McLaughlin et al., 2013; Joshi et al., 2015; Granger et al., 2017b; Li et al., 2017).

Initial application of this approach in central and western Lake Erie indicated the presence of two distinct sources of P_i; riverine inputs (~+12‰) and a second, isotopically heavier source (~+17‰) of unknown origin (Elsbury et al., 2009). Subsequent investigations have suggested that this P_i is derived primarily by redox mediated release from anoxic sediments in the central basin (Paytan et al., 2017). Given the apparent success at identifying different P_i sources in central and western Lake Erie and the residual uncertainty regarding the relevance of external vs. internal sources of P_i associated with Cladophora blooms, we sought to apply this approach to (1) examine the influence of the Grand River on nearshore phosphorus concentrations (2) investigate the utility of stable isotopes of oxygen in phosphate for differentiating P_i sources in the nearshore of the eastern basin and (3) to improve our understanding of P dynamics in a dynamic coastal environment. Our results presented here add to the expanding library of freshwater δ¹⁸O_P data, particularly at the oligotrophic end of the spectrum. In addition, we provide data from a eutrophic river and to our knowledge, the first measurements of δ¹⁸O_P generated by filter feeding bivalves. An improved understanding of the sources of P_i and the manner in which it is processed and recycled within the nearshore are of critical importance to establishing effective nutrient management strategies to ameliorate Cladophora blooms as required by the amended (GLWQA, 2012).

Materials and Methods

Study Region

The eastern basin of Lake Erie is the deepest of the three basins (64 m maximum depth) and is considered oligotrophic based on long term monitoring of phosphorus concentrations (mean spring TP ~0.32 μM; Dove and Chapra, 2015). The eastern basin receives ~40% of its total P load from the Grand River (Figure 1) while the remainder comes via exchange with the central basin and other smaller tributaries around the north and south shores of the basin (Maccoux et al., 2016). With prevailing winds from the southwest direction, the north shore littoral region is subjected to a large fetch and this highly exposed shoreline is dominated by limestone bedrock with smaller regions of glacial till and sand (Rukavina and St. Jacques, 1971). Dreissenid mussels have colonized most of the hard substrate along this shore and remain a dominant feature of the benthos (Patterson et al., 2005; Chomicki et al., 2016).

Water Sample Collection and Nutrient Analysis

Water samples were collected on four to five occasions in 2013 and 2014 between the months of April and September at 12 sites along the northern shore of the east basin of Lake Erie (Figure 1). Vertical profiles of temperature, specific conductivity, pH, dissolved oxygen, the physical structure of the water column were collected at each site with a YSI-6600 or EXO-2 water quality sonde prior to initiation of sampling. All water samples were collected using a submersible pump and stored in acid-cleaned 4 L HDPE jugs, placed on ice in coolers until processing for total phosphorus (TP), total dissolved phosphorus (TDP), soluble reactive phosphorus (SRP), ammonia (NH₃), nitrate/nitrite (NO₃/NO₂), total chlorophyll a (CHLA), particulate carbon (PC), and particulate nitrogen (PN). Water samples were processed for the above parameters following (Dove et al., 2009) and analyzed at the National Laboratory for Environmental Testing, Burlington ON. Particulate P (PP) was calculated as TP—TDP and dissolved organic P (DOP) was calculated as TDP-SRP.

For isotopic analyses of water and P_i, additional samples were collected with a target volume of ~200 L (to attain a minimum of 10 μM P_i) from selected depths using a submersible pump mounted on a tripod (Table 1): surface (S; 1.0 m depth) and ~0.1 m above the lake bottom (LB). To minimize boundary layer disturbance, the pump flow rates were set to ~50–60 mL min⁻¹, resulting in collection times of ~1 h. Pumped water was passed under low pressure first through a 1 μm polypropylene depth cartridge filter (Hytrex^®) and then through a 0.2 μm pleated polypropylene filter (Flotrex^®; FPN921BGS) to remove particles and colloidal material and collected in 10 20 L acid cleaned carboys. To each carboy, ~50 g of MgCl₂ (Fisher Scientific, ACS grade) was added and dissolved followed by an addition of 100 mL of 1M NaOH (Fisher Scientific, ACS grade) to achieve a loading of ~0.3% v:v to induce precipitation of brucite (MgOH₂) (Karl and Tien, 1992; Thomson-Bulldis and Karl, 1998). Fixed sample bottles were kept darkened and cool until further processing upon return to shore (see below).

Phosphate and Water Oxygen Isotope Analysis

To analyze δ¹⁸OP in P_i, the brucite floc from each set of carboys was allowed to settle before siphoning supernatant water away. Samples were combined into a smaller number of carboys and this process was repeated until ~10 L of brucite floc remained. Sufficient 1M HNO₃ (ACS grade) was added to dissolve the brucite floc, but keep the pH of the sample between 5 and 7. 1M NaOH was then added to re-precipitate brucite and reduce the overall sample volume. At least three cycles were required to reduce sample volumes from 200 L to ~50 mL. At each stage, SRP and TDP concentrations in both the sample and supernatant were measured after reagent addition to ensure >99% recovery of SRP (i.e., SRP in supernatant was not detectable). Recovery of SRP was always >99% while recovery of TDP varied between 5 and 53%, depending on the sample and amount of 1M NaOH required.

Initially, we sought to precipitate silver phosphate (Ag₃PO₄) following the protocol of McLaughlin et al. (2004). However, we were unable to generate Ag₃PO₄ that was sufficiently clean of contaminant oxygen using this approach. Consequently samples collected in 2013 were deemed contaminated and are not presented herein. For samples collected in 2014, we adapted the approach of Colman et al. (2005) to remove organic matter prior to undergoing sequential precipitation and crystallization steps. Briefly, concentrated samples (~50 mL) were run through Oasis HLB resin columns (1 mL min⁻¹ flow rate) to remove organic matter. P_i in the now cleaned sample was first precipitated as ammonium phosphomolybdate (APM) and then magnesium ammonium phosphate (MAP) (Kolodny et al., 1983; Tamburini et al., 2010) to separate P_i from remaining impurities. Following dissolution of the MAP precipitate, cations were removed with BioRad AGx50 cation resin, and samples were evaporated to a total volume of ~1.0 mL prior to precipitation of Ag₃PO₄. Ag₃PO₄ crystals were cleaned with H₂O₂ (15%; Tamburini et al., 2010), washed copiously with distilled water, dried at 60°C and stored in a dessicator prior to analysis. When possible, samples were processed in duplicate with ¹⁸O enriched reagents (±59.7‰) to assess potential for hydrolysis of organic P compounds. For the limited number of samples with sufficient phosphate recovered that we could process with enriched reagents (n = 9), the difference between spiked and un-spiked samples was on average < 1.5‰, thus we conclude that minimal hydrolysis occurred during processing and purification.

Between 0.25 and 0.35 mg of sample and nickelized carbon were added to silver capsules (pre-baked @ 500 C, 4 h) and shipped to the National Hydrology Research Centre in Saskatoon, SK. Samples were run in duplicate or triplicate if sufficient Ag₃PO₄ was recovered, but the majority of samples only had sufficient Ag₃PO₄ for a single analysis. Oxygen stable isotope analyses were accomplished by thermal conversion of Ag₃PO₄ to CO by reaction with glassy carbon at 1450°C in a Thermo Finnigan TC/EA device. Reaction gasses were separated by gas chromatography and the resultant CO introduced to a Thermo Finnigan Delta V Isotope Ratio Mass Spectrometer (IRMS). We used crimped silver capillaries (Longinelli and Nuti, 1973) containing amount equivalent aliquots of USGS UC03 (δ¹⁸O = +29.8‰) and IAEA VSMOW (δ¹⁸O = 0.0‰) to normalize raw δ¹⁸O values from the mass spectrometer to the SMOW scale. An internal laboratory standard (0.3 μM; K₂HPO₄; Fisher, ACS Grade) was prepared in 10 L of filtered (< 0.2 μm) Lake Erie water that had been previously stripped of P_i using brucite flocculation as above. Internal standards were run as samples and averaged +9.6 ± 0.3‰ (un-spiked reagents) and + 9.9 ± 0.1‰ (spiked reagents). Certified reference material B2207 (Isomass Scientific) measured 22.1 ± 0.5‰. Analytical precision averaged 0.6‰. All values are reported using standard delta notation relative to Vienna Standard Mean Ocean Water (VSMOW).

Samples for water δ¹⁸O_W were collected whenever samples were taken for δ¹⁸O_P by filtering ~2 mL of lake water through a 0.45 μm cellulose acetate syringe filter into 2 mL glass vials. Vials were stored on ice and refrigerated until analysis by Off-Axis Integrated Output Spectroscopy (OA-ICOS) on a Los Gatos Research DLT-100 Liquid Water Isotope Analyser. We used two internal references waters, ROD3 and INV1 (δ¹⁸O = −1.1 and −28.3‰, respectively), to calibrate the raw values to the SMOW scale. The analytical precision of δ¹⁸O_W values was determined by duplicate analyses of samples and internal references and is ± 0.2‰.

Dreissenid Mussel and Cladophora Sampling

Benthic surveys to collected dreissenid mussels and Cladophora were conducted at similar intervals in 2013 and 2014 at a larger number of stations in the study area as part of a targeted monitoring effort, however we only report data relevant to the stations identified in this study for 2014 as these surveys were conducted within 1 to 4 days of sampling for water chemistry and isotopes. Briefly, mussels and algae were harvested from three 0.15 m² quadrats at each station (with the exception of 1 m depth stations where a 0.125 m² quadrat was used). Mussels were removed by hand or by scraping while applying an airlift device to collect samples in a mesh bag. Samples were stored cool on ice until return to the laboratory in Burlington, ON and subsequently cleaned of debris and frozen prior to analysis. After lyophilization, mussels were counted, weighed, measured, and soft tissue removed to estimate densities of live individuals and biomass (as shell free dry mass, SFDM m⁻²). For sites with algae collected concurrently, algal filaments were first cleaned in the field by washing with lake water. Upon return to the laboratory, algal filaments were cleaned under DI water and subsequently frozen. Algal biomass was determined after lyophilization and expressed as dry mass (g DM m⁻²).

Dreissenid Incubations

Live dreissenid mussels were collected in July of 2013 and May of 2014 to aid in characterizing the δ¹⁸O_P of P_i excreted by mussels (δ¹⁸O_{P_DM}). Divers harvested intact mussels by hand into a fine mesh bag. Mussels were placed in ziplock bags with damp paper towel and stored on ice for transport to the laboratory in Burlington, ON. Upon arrival, mussels were gently separated from rocks, shell material and other debris, brushed clean with a toothbrush under running DI water. Between 75 and 150 individual mussels (10–30 mm valve length) were placed into acid washed 4 L polypropylene jars filled with filtered (< 0.2 μm) Lake Erie water. Jars were aerated with air stones, and mussels were allowed to defecate and excrete for 12–36 h under a 16:8 L:D cycle at ± 2°C of lake temperature. SRP concentrations were monitored in the chambers until sufficient P had accumulated for analysis (~10 μmol). Lakewater was gently decanted and siphoned from the jar, filtered through a 0.45 μm cellulose acetate filter, and processed as described above for δ¹⁸O_P.

Calculation of Isotopic Equilibrium Values for Phosphate

The theoretical isotopic equilibrium values for P_i (δ¹⁸O_{P_EQ}) can be calculated using P_i water fractionation equations derived by Chang and Blake (2015) using Equation (1);

δ¹⁸O_{P_EQ} values were calculated for each sample using the measured δ¹⁸O_w and temperature from the water column profiles. To assess uncertainty, we calculated δ¹⁸O_{P_EQ} using the minimum and maximum δ¹⁸O_w for the entire study period for lake and river sites respectively.

Fraction of River Water and Expected Nutrient Concentrations

To estimate the relative influence of the river on nutrient conditions at the study sites, we used a linear mixing model to estimate the fractional contribution of river water (FRW) at lake sites using the measured values of δ¹⁸Ow

where δ¹⁸O_ws is the water δ¹⁸O_w value at each site, δ¹⁸O_wr is the water δ¹⁸O_w value at the river mouth and δ¹⁸O_wl is the highest water δ¹⁸O_w value observed at the open lake sites. We elected to use the highest δ¹⁸O_w value for each survey period as there was some spatial and temporal variation in δ¹⁸O_w, possibly due to larger scale circulation features. For each site (i), we calculated the expected concentration of PP, TDP, and SRP assuming conservative mixing between the open lake and river;

Where C_P(i) is the expected concentration of P (PP, TDP, or SRP) at each site, C_P(l) is the concentration of P for the open lake end member, and C_P(r) is the concentration of P at the river mouth. We used the average concentrations for sites with FRW < 0.1 to derive the open lake end member values for each survey, and C_P(r) is represented by the phosphorus concentration at the river mouth.

Results

Influence of the Grand River and Patterns of Nutrient Variation in the Nearshore

δ¹⁸O_w values ranged from −12.6 to −8.6‰ at the river mouth but were more constrained at the lake sites, ranging from −7.9 to −6.0‰ (Table 1). The calculated fraction of river water (FRW) was variable over all surveys but was highest in late May—early June, reaching up to 0.8 at shallow stations ~2 km due east of the river mouth (Figure 2). While FRW was highest after the spring runoff period, values ranged from near 0 to ~0.35 between June and August surveys, particularly at stations east of the river mouth (Figure 2). FRW tended to increase with increasing river discharge (Figure 2) although the correlation was weak (r = 0.23, p < 0.1).

Nutrient concentrations were often higher and more variable in proximity to the Grand River mixing area. Total phosphorus, total dissolved phosphorus and soluble reactive phosphorus (TP, TDP, and SRP) were negatively correlated with distance from the river mouth (r = −0.42, r = −0.39, r = −0.34, p < 0.001, p < 0.001, p < 0.01, respectively). Nitrate (NO₃) and dissolved organic carbon (DOC) concentrations were also strongly negatively related to distance from the Grand River mouth (r = −0.48, r = −0.48, p < 0.001, respectively). Ammonia (NH₃) was the only nutrient parameter measured that did not display any spatial gradient in relation to the Grand River (r = −0.17, p > 0.2).

Despite apparently elevated phosphorus concentrations at sites close to the river, deviation from expected concentrations based on a simple two member mixing model indicated strongly non-conservative behavior for both particulate and dissolved fractions in late May—early June. Notably, the deviations were strongest as FRW increased suggesting rapid consumption/settling or uptake processes were active (Figure 3A). From mid-June to August, the same general trend was apparent for particulate P, but was somewhat attenuated for TDP and SRP, largely because the concentration difference of TDP and SRP between river mouth station and open lake sites was small (Figures 3B–D). Some positive excursions for particulate P were evident at 1 m depth sites (Figures 3B,C) suggesting possible resuspension of sediment or sloughed periphyton. During mid-June, July and August surveys, TDP and SRP were generally lower than predicted for LB samples from shallow water (1 and 3 m depth) sites, with the exception of PSN 1371 and 1274 in July (Figure 3C). TDP and SRP in LB samples at 10 m sites were often (but not always) elevated relative to expected concentrations predicted from mixing alone (Figure 3).

Dreissenid Mussel and Cladophora Distribution

Plots of Cladophora and dreissenid mussel abundance and biomass (SFDM m⁻²) are shown in Figure 4. Cladophora biomass ranged from 0 to 189 g DM m⁻², reaching maximum values at 1 m depth in June and July (Figure 4). Cladophora biomass at 3 m depth sites did not exceed 25 g DM m⁻² and at 10 m sites did not exceed 2 g DM m⁻² (Figure 4). Dreissenid abundance and biomass varied widely between surveys and locations (6–2131 m⁻² and 0.2–65 g SFDM m⁻² respectively). Dreissenid abundance and biomass was clearly higher at 10 m depth sites than 3 m sites. No strong spatial patterns in dreissenid abundance or biomass were observed, and although Cladophora biomass was high at the 1 m sites close to the Grand River, we do not have data for similar depths at greater distances to the west side of the river, which obscures the lack of spatial patterns observed in other years (D.Depew, unpubl.data).

Oxygen Isotopic Composition of River and Lake Samples

δ¹⁸O_P and calculated δ¹⁸O_{P_EQ} values for all stations and dates are shown in Table 1; and their spatial/seasonal and depth distributions in Figure 5 through Figure 7. δ¹⁸O_{P_EQ} values at the river mouth were isotopically lighter compared to lake sites (Table 1). The seasonal dynamics of δ¹⁸O_P at the river mouth are displayed in Figure 5. During freshet period (April-late May), despite large variation in river discharge, δ¹⁸O_P varied relatively little, ranging between +11.1 ± 0.6 and +11.9 ± 0.7‰ and overlapped the δ¹⁸O_{P_EQ} window at times (Figure 5). From mid-May to late August, δ¹⁸O_P declined steadily (δ¹⁸O_P in the river was significantly correlated with day of year; r = −0.95, p < 0.001) through the season reaching a low of +5.7‰ in late August (Figure 5).

The spatial and temporal variation of δ¹⁸O_P at lake sites are shown in Figure 6. δ¹⁸O_P at lake sites ranged from +4.7 to +14.9‰ and did not show strong relationships with FRW or with distance from the Grand River (Figure 6). Within-transect variation (i.e., variation across different depths) in δ¹⁸O_P was often greater than differences from west to east (Figure 6). No obvious differences were observed for samples collected near the lake bed at 1 and 3 m sites and surface samples at 10 m sites, but samples collected near the lake bed at 10 m sites (LB10) were consistently isotopically heavier than companion surface samples (paired t-test, t = 11.1, n = 12, p < 0.001) and frequently overlapped the range of δ¹⁸O_{P_DM} (Figure 7)_. Notably, all but one lake sample are below the calculated theoretical equilibrium based on water temperature and δ¹⁸O_W (Figure 6, Table 1).

The degree of isotopic offset of δ¹⁸O_P from the calculated δ¹⁸O_{P_EQ} was strongly correlated to SRP concentration and negatively correlated to the seston N:P ratio, suggesting larger negative excursions for samples associated with low P availability (Figure 8). Such patterns are expected if the hydrolysis of DOP is an important mechanism of P_i regeneration. Indeed, we observed low δ¹⁸O_P values and larger negative excursions from the δ¹⁸O_{P_EQ} values as the TDP pool became increasingly dominated by DOP (Figure 8). The only samples that appeared to not follow this pattern were those collected at the lakebed at 10 m depth sites.

Oxygen Isotopic Composition of Phosphorus Excreted by Dreissenid Mussels

Oxygen isotopic compositions of phosphate collected during dreissenid excretion experiments (δ¹⁸O_{P_DM}) ranged from +10.2 to +13.5‰ and averaged +11.9 ± 1.1‰ (Table 2). δ¹⁸O_P at LB10 and LB3 sites were positively related to dreissenid abundance and shell free biomass (Figure 9).

Discussion

Variation in water quality and nutrient conditions along the north shore of the eastern basin are shaped to a large degree by the discharge of nutrient rich water from the Grand River (Chomicki et al., 2016) as well as interaction with hydrodynamic features. Wind forcing (He et al., 2006) and large scale circulation features (Rao and Schwab, 2007) often trap river inputs alongshore, creating a band of enriched water that extends to approximately the 5 m depth contour (Chomicki et al., 2016). Estimates of FRW at our study sites are in general agreement with prior studies and existing knowledge and observations of nearshore circulation and river plume dynamics. For example, we observed FRW estimates approaching 80% at sites to the east of the river mouth, with larger values typically encountered at 1 and 3 m depth sites which are < 1 km from the shoreline. It is unclear if the high FRW observed ~20 km to the east represents input from the Grand River as there are a number of smaller tributaries that discharge along the shore between the Grand River and the eastern-most study sites. Previous estimates of river influence extent range from 8 to 10 km (Nicholls et al., 1983; Howell and Hobson, 2003) although these are generally based on seasonal data and use other semi-conservative tracers. Discharge from the Grand River in April and May of 2014 exceeded the long term (1913–2016) 75th percentile of ~300 m³ s⁻¹ by nearly a factor 2−3 (https://wateroffice.ec.gc.ca/report/historical_e.html?mode=Graph&type=h2oArc&stn=02GB001&dataType=Daily&parameterType=Flow&year=2014&y1Max=1&y1Upper=1&scale=normal), and given that river water temperatures at this time were approximately double those encountered in the lake (~19.5°C compared to 11–12°C), the river plume would have likely been thermally buoyant, spreading thinly over large distances before mixing completely. The FRW was less variable during low flow periods ranging from 0 to ~20%. FRW was generally higher at shallow stations (i.e., the 1 and 3 m depth sites), consistent with the alongshore entrapment of the plume.

Although nutrient concentrations were generally elevated closer to the river mouth and the river water clearly influenced our study sites, we observed strong non-conservative behavior of all forms of phosphorus across the river mouth—lake mixing area. During late May—early June, when concentrations of PP, TDP, and SRP were highest at the river mouth, the measured concentrations in the lake were considerably lower than would be expected based on simple dilution, even at sites in close proximity (< 2 km) with FRW approaching 80% (Figure 3). In fact, losses of PP, TDP, and SRP appeared to be strongest at the sites closest to the river mouth, suggesting rapid sedimentation, consumption, and/or uptake (Figure 3A). Losses of particulate P could reflect the settling of large and coarser organic P particles owing to reductions in flow velocity and turbulence, or perhaps consumption of organic P by dreissenid mussels (Chomicki et al., 2016). Smaller particles such as non-apatite inorganic P (~4–7 μm diameter; Chomicki et al., 2016) would likely remain in suspension for a greater duration and be subjected to more widespread mixing and dispersion. Declines in SRP and to an extent TDP may reflect uptake by phytoplankton as chlorophyll a concentrations increased from ~0.5 μg L⁻¹ at open lake sites (e.g., 502 and 1350) to a maximum of 15.3 μg L⁻¹ at PSN 1371, 2 km east of the river mouth. Similar results have been observed across the Milwaukee River mixing zone in Lake Michigan, where biological processes have been implicated in removal of SRP in close proximity to the point of discharge (Lin and Guo, 2016).

Over the late June to August surveys, losses of particulate P remained consistently strong in close proximity to the river mouth although there were instances of resuspension elevating PP concentrations, particularly at 1 m depth sites (Figure 3B). The responses of TDP and SRP were generally consistent, although the magnitude of difference between observed and expected concentrations was small owing to the relatively small difference between TDP and SRP concentrations at the river mouth and lake sites. At times, TDP and SRP concentrations were larger than expected based on the mixing model calculations. Elevated concentrations of SRP were often observed for LB samples at 10 m sites (Figures 3B–D), but were also observed at 1 m and 3 m sites in July (Figure 3C). The proximate cause of these elevated concentrations could be due to leakage of nutrient rich groundwater (e.g., Robinson, 2015) or perhaps due to the decomposition of Cladophora that had accumulated at the adjacent shoreline as decomposing Cladophora is known to lose phosphorus slowly over a period of several weeks in situ (Paalme et al., 2002).

Oxygen Isotope Ratios of Phosphate Indicate Extensive Phosphorus Cycling

Compared to marine systems (Colman et al., 2005; McLaughlin et al., 2006a,b, 2013; Goldhammer et al., 2011; Jaisi et al., 2011; Joshi et al., 2015; Li et al., 2017) and soils (Zohar et al., 2010; Angert et al., 2012; Tamburini et al., 2012; Granger et al., 2017a), freshwater environments are poorly characterized with respect to δ¹⁸O_P [but see (Elsbury et al., 2009; Granger et al., 2017b; Pistocchi et al., 2017)]. The δ¹⁸O_P we measured in the Grand River during April and early May (+9.4 ± 1.7 to +11.9 ± 0.7‰) are similar to values reported for Lake Erie tributaries in the central and western basins over summer and fall (Elsbury et al., 2009), but lower than those reported for the River Taw, UK (Granger et al., 2017b). δ¹⁸O_P at the river mouth site was near or overlapped δ¹⁸O_{P_EQ} in late April and early May, but declined steadily from +9.4 ± 1.7‰ to a minimum of +5.7‰ through August, with the departure from δ¹⁸O_{P_EQ} increasing to the end of August(Figure 5). Similarly, δ¹⁸O_P at lake sites are considerably lower than those reported by (Elsbury et al., 2009) for the west and central basin of Lake Erie (~+11 to +17‰), and all but one LB sample were below the calculated δ¹⁸O_{P_EQ} based on δ¹⁸O_W and water temperature.

It is important to keep in mind that δ¹⁸O_P measured at a given site at a given time represents a snapshot of a variety of processes occurring along a continuum of P cycling. For example, the P_i pool may be affected to varying degrees by the intensity of (1) P_i that has been equilibrated internally and released by microbes and/or phytoplankton (Blake et al., 1997; Paytan et al., 2002), (2) P_i excreted by zooplankton (Taylor and Lean, 1981), (3) remineralization of DOP excreted by phytoplankton (Bentzen et al., 1992) or liberated via zooplankton grazing and/or viral lysis (Taylor and Lean, 1981), and (4) mixing with an external P_i source with a distinct δ¹⁸O_P (Paytan et al., 2017). The strong negative deviations from theoretical equilibrium values observed here are anticipated in systems where ambient levels of P are limiting and arise due to the large inheritance effects associated with isotopic fractionation during the incorporation of water oxygen into the P_i ion during enzymatic hydrolysis (Liang and Blake, 2006, 2009). Forms of DOP commonly found in aquatic ecosystems include phosphomonoesters, phosphodiesters, and phospholipids produced by organisms in situ (Baldwin, 2013) but may also include phosphomonesters, phosphodiesters and inositol phosphates exported from the surrounding catchment (Monbet et al., 2009). Monoesters are generally considered bioavailable due to the presence of enzymes such alkaline phosphatase (EC 3.1.3.1) (Bentzen et al., 1992), but phosphodiesterases and phytases may make diesters and inotisol phosphates potentially bioavailable (Cheng and Lim, 2006). Although we do not have detailed information on the nature and composition of DOP or enzyme activity in eastern Lake Erie (but see North et al., 2012), if we assume that organic P contained within organisms in the lake is at or near equilibrium (e.g., Paytan et al., 2002, based on known fractionation factors for various DOP compounds, the δ¹⁸O_P of released P_i recycled from the DOP pool is expected to be < ~+14‰, and even lower < +11‰ if substrates such as phytate is not a major component of the bioavailable DOP pool (Table 3). While these are generalizations, they illustrate the relative magnitude expected if regeneration of P_i from DOP is an important P recycling pathway. This is further supported by the significant relationships we observed between δ¹⁸O_P and other proxies of P availability such as %DOP, SRP concentration and seston N:P ratios (Figure 8), and we infer that the recycling of DOP by the plankton community is a major determinant of the isotopic composition of the available P_i pool in eastern Lake Erie.

Oxygen Isotope Ratios Cannot Be Used to Track Riverine Input to Eastern Erie

Concentrations of SRP at the river mouth site were high (up to 6.9 μM) between April and early June. Despite large variation in river discharge volume, δ¹⁸O_P was relatively invariant (Figure 5). This is somewhat surprising given that land use in the Grand River watershed is a mixture of intensive agricultural and urban use (Holeton, 2013). Although we do not have a detailed inventory of δ¹⁸O_P of different sources within the Grand River watershed, based on existing data it seems likely that δ¹⁸O_P might vary according to the relative inputs of P_i from different sources (e.g., Young et al., 2009). In addition,the short residence time of the river (~3 d, Rosamond, 2013) and low temperatures (< 10°C) should have acted to minimize any uptake and cycling of P_i in the river at this time. While it is possible that δ¹⁸O_P within the watershed might not vary across different sources, or perhaps simply reflected thorough mixing of source inputs, we note that our SRP concentrations measured at the river mouth site were up to 2 to 6-fold higher compared to concentrations measured the same day at a nutrient monitoring station ~40 km upstream (Supplemental Information Figure S1), indicating a significant input of P_i below the monitoring station. While further work will be required to identify the proximate source of this P_i, we note that this section of the river is characterized by extensive wetland complexes, some of which at times can have very high TDP concentrations (6–190 μM; (Gilbert and Ryan, 2007), and may explain lack of variation in δ¹⁸O_P during this period.

Despite high concentrations of SRP, a relatively constant δ¹⁸O_P and high discharge during the April to early June we did not observe any evidence of mixing between riverine δ¹⁸O_P and lake δ¹⁸O_P at any of our sites. During the cooler spring months, P uptake by phytoplankton is generally reduced relative to summer months (Lean et al., 1983), yet at sites closest to the river mouth (i.e., those with high FRW), δ¹⁸O_P was no different than lake sites far removed from the river mixing zone. Only LB samples from 10 m depth sites had δ¹⁸O_P that was reasonably close to riverine δ¹⁸O_P, but these samples had low FRW and are likely affected by recycling of P by dreissenid mussels (see below).

Between June and August, δ¹⁸O_P at the river mouth site declined steadily, departing further from δ¹⁸O_{P_EQ} (Figure 7). Such observations are not uncommon in P–rich systems (Paytan and McLaughlin, 2012), however, δ¹⁸O_P moved further away from δ¹⁸O_{P_EQ} as SRP concentrations declined from 6.9 μM in late May to a low of 0.08 μM in July, which is counter to expectation (McLaughlin et al., 2006b; Granger et al., 2017b). On one hand, this could indicate an input of an isotopically light source of P_i such as waste water treatment plant (WWTP) effluent (Young et al., 2009) or P_i liberated via the degradation of glyphosate (Li et al., 2016a). During summer base flow periods, WWTP inputs become a larger component of the total P mass load to the river (Holeton, 2013), however, δ¹⁸O_P of treated WWTP effluent from plants within the Grand River watershed are reported to be higher than values we measured [+11.4 to +22.8‰ (Morrison, 2014)] and P uptake studies downstream of WWTP discharges indicate relatively short uptake lengths (< 5500 m; Barlow-Busch et al., 2006). Similarly, while glyphosate is a common herbicide in southern Ontario, concentrations typically peak during rain events (Struger et al., 2015) and since we observed low δ¹⁸O_P at near base flow conditions, this is unlikely to explain the values measured here. Rather, we speculate that P_i is extensively processed within the lower river by the attendant phytoplankton community. Previous studies of the lower Grand have indicated that an appreciable portion of dissolved P is converted into phytoplankton biomass prior to entering Lake Erie (Kuntz, 2008). Despite the overall eutrophic conditions, the increasing disequilibrium from δ¹⁸O_{P_EQ} at the river mouth during the summer in combination with high chlorophyll a concentrations (25–48 μg L⁻¹) is not inconsistent with a highly productive river estuary where even the use of DOP can be important (Larson et al., 2016; Yuan et al., 2017). Further study of phosphorus sources to the lower river and speciation would be required to confirm such pathways, particularly in order to trace back to specific input sources of P. Nonetheless, even though riverine δ¹⁸O_P appears to be altered prior to discharge into the lake, when our δ¹⁸O_P data are considered together with evidence of strongly non-conservative mixing behavior of SRP, we infer that most of the P_i discharged by the river is quickly assimilated into phytoplankton and rapidly overprinted via DOP recycling.

Isotopic Enrichment at the Lake Bottom Indicates Benthic P Recycling

Samples collected near the lake bed at 10 m depth sites were persistently isotopically heavier than paired surface samples at all sites on all dates, and most samples overlapped the δ¹⁸O_P measured in our excretion experiments (δ¹⁸O_{P_DM}) (Figure 7). In some cases, differences between lakebed and surface SRP concentrations were large (up to 10-fold) indicating enrichment of the near bottom waters while in others, the differences were smaller. Overall, SRP was significantly higher in LB samples compared to surface samples (paired t-test, t = 2.24, p < 0.05) and chlorophyll a concentrations were lower (paired t-test, t = 2.89, p < 0.05) consistent with known impacts of dreissenid grazing of phytoplankton and excretion of soluble P (Ackerman et al., 2001; Conroy et al., 2005; Ozersky et al., 2009). Differences in TDP displayed similar patterns but were not significant (paired t-test, t = 1.69, p < 0.1) and differences in NH₄, NO₃ and PC, PN or PP were not different (data not shown).

The experimentally determined δ¹⁸O_{P_DM} was constrained relative to δ¹⁸O_P collected at lake sites, regardless of where or when dreissenid mussels were collected for incubations. Dreissenid mussels excrete P_i directly (Arnott and Vanni, 1996), DOP (Mosley and Bootsma, 2015), and deposit a larger quantity of nutrients (including P) as partially digested pseudofaeces and feces on the surrounding sediments/lakebed (Nalepa et al., 1991). The smaller range of δ¹⁸O_{P_DM} compared to δ¹⁸O_{P_}from LB samples might be due to differences in the quantities of P_i derived from direct excretion vs mineralization of pseudofecal and fecal material, as our incubation time required to generate sufficient P_i for recovery (~10 μM; 12–36 h) were of sufficient duration to permit accumulation of feces and psuedofaeces, and although we did not measure microbial activity, decomposition of fecal and pseudofecal material generally occurs over a longer time frame than employed here (Roditi et al., 1997). In contrast, our LB samples could contain P_i excreted directly by dreissenids, or derived from the re-mineralization of pseudofaeces and/or feces. P_i excreted directly by dreissenids is likely liberated through the action of acid phosphatases (EC 3.1.3.2) in the acidic gut (Sauey et al., 2016), and based on examples in Table 3, released P_i could range between +7.0 and +14‰. Many of our δ¹⁸O_{P_DM} fall within this range, as do δ¹⁸O_P from LB samples at 10 m sites. The generally consistent overlap between δ¹⁸O_{P_DM} from incubations and δ¹⁸O_P from LB sites is further supported by positive relationships observed between the degree of offset of δ¹⁸O_P from our experimentally determined δ¹⁸O_{P_DM} and dreissenid abundance and shell free biomass (Figure 9) and support the role of dreissenid mussels as active recyclers of P in the benthic environment.

Alternative explanations for isotopic enrichment of δ¹⁸O_P near the lake bed could include possible desorption of P_i from clays or iron minerals, or injection of isotopically heavier P_i from the hypolimnion. Chomicki et al. (2016) suggested that ~50% of the particulate P leaving the Grand River is non-apatite inorganic P (NAIP), of which iron oxides are a dominant component (Mayer and Manning, 1989). While much of this material likely settles out within a short distance of the river mouth, it is nonetheless subject to resuspension and subsequent dispersal by wind and waves (Chomicki et al., 2016). The structurally complex and rough surfaces created by abundant dreissenid mussels may retain settled material in interstitial spaces and enhanced O₂ demand from fecal material and respiratory activity of mussels (Turner, 2010) could promote release of P_i from redox sensitive P bearing sediments. While we did not directly measure the δ¹⁸O_P of redox sensitive NAIP, we consider this to be an unlikely explanation for the persistently heavy δ¹⁸O_P near the lake bed. Jaisi et al. (2010) have shown that P_i adsorbed to mineral surfaces are constantly exchanged with P_i in the surrounding water such that fractionation associated with adsorption/desorption is negligible over a matter of hours. Since NAIP particles near the Grand River are small (4–7 μm; Chomicki et al., 2016) and likely to stay in suspension for several days (DePinto et al., 1981), we would expect P_i retained on NAIP particles that eventually settles within mussel beds to be equilibrated with δ¹⁸O_P in the water column. If P_i is released via redox mediated dissolution, we infer that such P_i would retain δ¹⁸O_P closer to that of lake water, and would be unlikely to explain the higher δ¹⁸O_P values observed near the lake bed.

Entrainment of cooler water from the hypolimnion could theoretically account for near bed enrichment via (1) remineralization of sinking P by bacteria or (2) release of P_i from anoxic sediments. In the open ocean, δ¹⁸O_P shifts toward the theoretical equilibrium with increasing depth as sinking organic matter is re-mineralized and microbes become carbon limited (Colman et al., 2005). Demand for P_i by phytoplankton and activity of hydrolytic enzymes should also be reduced in the cooler and darker hypolimnion compared to the illuminated surface layers (Pick, 1987). Cooler temperatures would favor an isotopically heavier δ¹⁸O_{P_EQ}. Upwelling has been identified as a potentially significant transport mechanism for hypolimnetic SRP in some years (Valipour et al., 2016) which could affect our measures of δ¹⁸O_P. Wind speed and direction data from buoy #C45142, prior to and during our sampling trips did not indicate the presence of conditions favorable for upwelling (Supplemental Figure S2), nor did we observe any deviations in water temperature profiles during surveys associated with hypolimnetic intrusions (data not shown), thus we infer that upwelling cannot explain the heavier δ¹⁸O_P observed near the lake bed. Similarly, Paytan et al. (2017) have linked release of isotopically heavier δ¹⁸O_P in the central basinto sediment release during hypoxic events. However, hypolimnetic oxygen concentrations in the much less productive eastern basin rarely fall below 0.3 mM during the stratified period (Carrick, 2004). Consequently, high rates of sediment P release from internal loading is unlikely to explain the δ¹⁸O_P measured at our 10 m sites.

The lack of a similar degree of isotopic enrichment in the LB samples at 1 and 3 m sites is not easily explained. On one hand, the relationship between δ¹⁸O_P and dreissenid biomass might suggest that dreissenid biomass at these shallow depths was insufficient to generate a strong isotope signal. Rapid uptake of released P_i by Cladophora might also attenuate an isotopic signal, however, Cladophora biomass in 2014 was relatively low, in large part because water clarity in 2014 was poor owing to high inputs of suspended solids from the Grand River (Valipour et al., 2016). It is also possible that the length of time required (~1 h) to collect sufficient sample volume (~200 L) whilst minimizing the disturbance near the lake bed may have compromised our samples to varying degrees if vertical mixing was sufficient to mix surface waters (containing lighter δ¹⁸O_P) with bottom waters (containing heavier δ¹⁸O_P). We evaluated this possibility by calculating the characteristic time scale for diffusive mixing of the water column following Boegman et al. (2008)

where N is the buoyant frequency (calculated as $\sqrt{-g/{\rho }_0(\frac{\partial \rho }{\partial z})}$, H is the water column depth (m), and u_* is the surface shear velocity (calculated as $\sqrt{\frac{{\tau }_w}{{\rho }_0}}$; where τ_w can be estimated from $C_D{\rho }_A{U}_{10}^2{e}^{-\beta z}$; with C_{D =} 1.3 × 10⁻³, ρ₀ = 1,000 kg m⁻³, ρ_A = 1.2 kg m⁻³, and U₁₀ is the average daily wind speed (m s⁻¹) adjusted to a height of 10 m above surface recorded at buoy #C45142). Note since τ_w was evaluated at surface, the exponential term is neglected (Boegman et al., 2008). For many samples, it is clear that the time to completely mix the water column was shorter than our sample collection time (Figure 10)—thus we cannot exclude the possibility that our samples may have been compromised by mixing of lower δ¹⁸O_P from surface waters. Furthermore, the degree of δ¹⁸O_P offset from δ¹⁸O_{P_DM} appeared to be reduced as $T_D^{mix}$increased for samples collected in May and June when density gradients were greatest (Figure 10). Although the low sample size likely precludes statistical significance, these calculations are consistent with our hypothesis.

Conclusions

One of our original objectives was to assess the utility of δ¹⁸O_P to delineate the impact of riverine P_i in the nearshore and determine any potential association with Cladophora blooms. This proved to be challenging as the δ¹⁸O_P of riverine P_i appeared to be rapidly altered by in-river and in-lake processes. Based on the nature and magnitude of deviations from δ¹⁸O_{P_EQ}, δ¹⁸O_P, we infer that internal recycling of DOP by the plankton community is a dominant pathway determining the overall δ¹⁸O_P of the available P_i pool in eastern Lake Erie. In addition, our data provide compelling evidence for the release of P_i at the lake bed, facilitated by dreissenid mussels. Management actions aimed at reducing the extent and severity of Cladophora blooms in Lake Erie will likely require consideration of linking the magnitude of external P inputs to the in-lake growth response in a more complex framework than prior to the advent of dreissenid populations. The ability of dreissenid mussels to function as a P trap in the nearshore as proposed by Hecky et al. (2004) poses distinct challenges. Dreissenid mussels in eastern Lake Erie are highly reliant on phytoplankton (Campbell et al., 2009), and the open waters represents a large reservoir of P that would have been previously unavailable to Cladophora. At the same time, local inputs of P_i from tributaries may be rapidly incorporated into phytoplankton, which while formerly unavailable to Cladophora, must now pass through the littoral zone which is dominated by expansive mussel populations and such P may ultimately be trapped and made available in the nearshore through the feeding activity of mussels. While it is clear that reductions in P inputs are the appropriate management tool, further work will be required to better quantify if reductions in Cladophora growth are achievable and whether these changes will require more stringent P controls than already proposed for Lake Erie.

Author Contributions

DD and VH-B designed the study. DD carried out the study. GK conducted isotope analysis. DD wrote manuscript.

Funding

This study was funded by Environment Canada through the Great Lakes Nutrient Initiative.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors thank the Technical Operations Services at Environment and Climate Change Canada for field and laboratory support. Greg Boulton, Karina Tsang, Danielle Bruyns, and Jacqui Milne provided tireless assistance in the field. Alice Dove provided nutrient data collected at York, ON, and data on dreissenid abundance and biomass Sandra Cooke (Grand River Conservation Authority) provided valuable input on the Grand River ecology. Comments and advice regarding purification protocols from Federica Tamburini, Adina Paytan, and Deb Jaisi are appreciated and were essential for improving sample processing.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2018.00215/full#supplementary-material

References