Assessing the utility of sulfur isotope values for understanding mercury concentrations in water and biota from high Arctic lakes

The six lakes were located on Cornwallis Island in the central Canadian Arctic Archipelago near the hamlet of Resolute Bay, NU, Canada; a map is provided in the Supplementary Material (Fig. S5).¹ This region has a polar desert climate, with low mean annual precipitation and temperature (150 mm and −16.4 °C; Antoniades et al. 2011). Detailed water chemistry and lake morphological characteristics are reported in Lescord et al. (2015b). Briefly, these lakes range in size from 39 to 130 ha and in maximum depth from 8 to 27 m, and they are ultra-oligotrophic, with mean chlorophyll-a, total phosphorus (TP), and dissolved organic carbon (DOC) concentrations ranging from 1.28 to 2.23 μg/L, 3.6 to 8.0 μg/L, and 0.5 to 1.9 mg/L, respectively. Dissolved oxygen and temperature profiles showed no thermal or oxygen stratification in the water columns of these systems (Lescord et al. 2015b). Mean pH was similar across the lakes, ranging from 7.9 to 8.1, whereas mean aqueous in surface waters was low but varied by approximately six-fold (0.23–1.35 mg/L across lakes). The food webs have low biodiversity; Arctic char are the only fish species, the benthic invertebrate community consists mainly of chironomids, and the zooplankton community is dominated by copepods (Chételat and Amyot 2009). The lakes from this study (and therefore the char that reside in them) were assumed to be land-locked with no direct marine connection (Lescord et al. 2015b).

Four of the lakes are generally considered un-impacted by human activity (Char, Small, North, and 9-Mile lakes), whereas the catchments of Meretta and Resolute Lakes include the Resolute Bay airport and a former military air base; from 1949 to 1998, untreated waste water from the air base was released directly into the catchment of Meretta Lake (Antoniades et al. 2011), which flows into Resolute Lake. Although recovery from this eutrophication has occurred, these inputs affected both the chemical and microbial characteristics of Meretta Lake (Antoniades et al. 2011). Furthermore, the proximity of both Meretta and Resolute Lakes to the local airport has altered their organic contaminant profiles when compared with more remote lakes in the area (Lescord et al. 2015a; Cabrerizo et al. 2018).

All field sampling was conducted between July and August of 2010 and 2011. Details of sampling for water chemistry (e.g., MeHg, total Hg (THg, all forms of Hg in a sample), nutrients, and pH), biota (Arctic char, larval and adult chironomids, and zooplankton) and sediments were reported by Lescord et al. (2015b). Briefly, water samples were collected weekly above the deepest part of each lake, approximately 1 m below the surface, using a Niskin bottle (model 1010, 2.5 L; General Oceanics, Miami, FL, USA). Larval and adult chironomids were collected by kick sweeping shore lines and aspiration from the ice surface, respectively. Periphyton scrapings were taken from littoral rocks in triplicates. Bulk zooplankton samples were collected by towing a net (80 μm mesh) for 3–5 min on each lake; these samples were then size filtered at <250, 250–500, and >500 μm. Large char (>18 cm total length) were collected using gill nets and small char (<18 cm total length) were electrofished along shorelines; these size categories of fish were kept separate for all laboratory and statistical analyses. Fish sampling was done using protocols approved by the Animal Care Committee of the University of New Brunswick (#2010/11-1D-01).

For δ³⁴S analysis, water samples were collected from the surface and deep point of each lake (n = 1 replicate per lake per depth, in each of 2010 and 2011); these data were averaged across depths and years (n = 4 per lake). Samples were kept cool in the field and sulfate was extracted within 6 h of collection at the Polar Continental Shelf (PCSP) base by ion exchange as barium sulfate (Carmody et al. 1998). Surface sediment samples from the Arctic lakes were taken from the top 0.5 cm slice of cores (n = 18, three replicates per lake in 2011) that were collected in the littoral zone of each lake as part of other studies (Lescord et al. 2015a, 2015b).

All laboratory analyses were performed on individual muscle subsamples for large char, whole body homogenates for small char, bulk zooplankton (pooled by sample date, location, and filtration size, n = 2–5 per lake), and benthic invertebrates (all in the family Chironomidae and pooled by sample date and location for adults and larva separately, n = 1–5 per lake). Biotic samples were freeze-dried and hand-ground to a fine powder before all analyses.

Biotic sulfur stable isotope analysis was conducted at the Colorado Plateau Stable Isotope Laboratory (Northern Arizona University, Flagstaff, AZ, USA). For these samples (n = 157), certified reference materials (International Atomic Energy Agency (IAEA) S1, S2, S3, S4, S05, S06) were used and percent differences from expected values (ranging from −34.05‰ for IAEA S06 to 22.67‰ for IAEA S2) averaged 0.82% ± 22.8% (total n = 21). Duplicate samples were also run and the average percent difference between duplicates was 4.21% ± 2.2% (n = 11). Sulfur isotopes (n = 20) and the percentage of sulfur (%S; n = 23; see Table S1¹) in sediments were measured at IsoAnalytical (Crewe, UK). The percent difference from certified reference materials (CRMS; International Atomic Energy Agency IA-R061, IAEA-SO-5, NBS-1577B) averaged 8.49% ± 18.2% (n = 15) and sample duplicates averaged 1.31% ± 3.42% (n = 5). Delta ³⁴S values in aqueous sulfate were measured as barium sulfate at the University of Waterloo Environmental Isotope Laboratory (Waterloo, ON, Canada; average percent difference between duplicates = 0.41% ± 7.9%).

Procedures for δ¹³C, δ¹⁵N, and Hg analyses in water, sediments, and biota were reported by Lescord et al. (2015b), along with associated quality assurance measures. Briefly, stable C and N isotopes were analyzed at the Stable Isotopes in Nature Laboratory (University of New Brunswick, Fredericton, NB, Canada) via Continuous Flow-Isotope Ratio Mass Spectrometry (CF-IRMS). Data on the elemental content of carbon, nitrogen, and sulfur of biota and sediment samples (as percentages of dry weight) were obtained in association with stable isotope analyses for each element. Percent recoveries of CRMs during CN isotope analysis ranged from 97.5% ± 5.8% to 105.1% ± 4.6% (n = 162) and duplicates were <1% different from each other.

Water samples were analyzed for THg and MeHg at the Low-Level Mercury Analytical Laboratory at the Canadian Centre for Inland Waters (CCIW, Environment Canada) following US EPA Methods 245.7 and 1630, respectively. All chironomids and zooplankton were analyzed for [MeHg] using a potassium hydroxide (KOH) digest and detected using Cold Vapour Atomic Fluorescence Spectroscopy (CVAFS) at the Center for Analytical Research on the Environment (CARE) at Acadia University (Wolfville, NS, Canada). Methyl Hg concentrations in blanks were low and consistent (average concentration = 0.01 ± 0.09 ng/g, n = 17) and percent recoveries of CRMs were high (101.1% ± 18.9%, n = 33). Small and large char were measured for [THg] (92%–99% of which was MeHg; Lescord et al. 2015b) using a Milestone Direct Mercury Analyzer (DMA-80) at the University of New Brunswick (Saint John, NB, Canada). Blanks were also low during THg analysis (average concentration = 0.002 ± 0.002 μg/g, n = 66) and 97.5% ± 11.8% of CRM certified values were obtained (n = 108). Total [Hg] in muscle from large char were converted to whole-body equivalent concentrations as per Peterson et al. (2004).

All statistical analyses were conducted using R version 3.4.3 and SigmaPlot version 11, and alpha (α) was set at 0.05 for all tests. Although zooplankton were separated by size factions (<250, 250–500, and >500 μm) in the field, samples were pooled for analysis due to biomass constraints. To account for potential variation in baseline δ¹⁵N values (Jardine et al. 2006), δ¹⁵N values were adjusted (δ¹⁵N_adj) for each biotic sample by subtracting the mean δ¹⁵N value of periphyton from a given lake (Lescord et al. 2015b, n = 3–5 per lake). Delta³⁴S and δ¹³C values were not baseline-adjusted. Because elemental composition has been shown to affect MeHg accumulation in Arctic biota (Chételat et al. 2012), we also considered the percentage of S, C, and N (%S, %C, and %N) of biotic samples in subsequent statistical analyses.

Log₁₀-transformed [Hg] (MeHg in invertebrates, THg in fish) were regressed against corresponding δ³⁴S values within biotic groups as well as in surface waters (MeHg and THg, separately). These relationships were tested using a linear mixed effect model (LME) through the nlme package in R (Pinheiro et al. 2018) including lake as a random factor to account for site-specific variability in the data. To compare the regression slopes of biotic [Hg] versus δ³⁴S between the six lake food webs, analysis of covariance (ANCOVA) models (Log₁₀ Hg [THg (fish, whole-body estimates) and MeHg (invertebrates), dry wt.] = lake + δ³⁴S + δ³⁴S × lake) were run, using biotic Hg as the dependent variable.

To further assess how biotic Hg concentrations related to sulfur isotopes in combination with other elemental measures (δ¹⁵N_adj, δ¹³C, %N, %C, and %S), penalized regression modeling was performed. Specifically, we used the Least Angle Shrinkage and Selection Operator (LASSO) method through the R package lars (Hastie and Efron 2013), similar to Lescord et al. (2018). Four sets of analyses were run: [Hg] in either chironomids (MeHg), zooplankton (MeHg), small char (THg), or large char (THg) were used as dependent variables with six elemental or isotopic predictor variables (δ¹⁵N_adj, δ¹³C, δ³⁴S, %N, %C, and %S). LASSO, a form of least angle regression, performs forward selection of predictors by incrementally reducing a penalty parameter (lambda, λ), allowing more variables to enter the model. Predictors with less influence on a dependent variable retain a coefficient of zero, excluding them from the final model estimates. We chose to use the LASSO method because elemental measures in this study were somewhat correlated (e.g., variance inflation factors of 1.5–27.9 in linear models) and LASSO has been shown to produce more accurate models when sample sizes are low and predictor variables are correlated, compared with other linear methods (Dormann et al. 2013; Hebiri and Lederer 2013; James et al. 2013). Mallow’s Cp values, an information theory measure similar to Akaike Information Criterion (AIC), was used to select the final LASSO models; the model with the lowest Mallow’s Cp value (>2, to prevent overfitting models) was chosen (James et al. 2013). Although the inclusion of a predictor variable in a final LASSO model indicates an effect on the dependent variable (i.e., [Hg] in this study), the R package selectiveInference was also used to estimate the significance of each predictor included in a given model at the corresponding level of λ (Tibshirani et al. 2017). LASSO path plots as well as correlation matrices are included in the Supplementary Material.¹

Average values of δ³⁴S in sediment and aqueous ranged from −13.04‰ to −3.43‰ and from −10.03‰ to −3.10‰, respectively, and absolute differences between these compartments were not consistent across lakes (Fig. 1). The values for Arctic lake sediments overlapped with those of boreal lakes (∼−5‰ to −2‰, Quebec lakes; Croisetière et al. 2009), but not with those of coastal hemiboreal systems (9.6‰–12.5‰, Nova Scotian lakes; Clayden et al. 2017). In addition, δ³⁴S values of aqueous sulfate in the Arctic systems fell within the typical range of terrestrial biogenic and anthropogenic sulfur emissions (Wadleigh et al. 1996) and did not overlap with those from either boreal or temperate systems (∼4‰–10‰ and 10.4‰–10.7‰, respectively), the latter of which closely matched that of in precipitation (Croisetière et al. 2009; Clayden et al. 2017). Mean annual precipitation near these temperate lakes is almost 10-fold higher than at the Arctic lakes (1332 mm versus 150 mm, respectively; Zhang et al. 2010; Antoniades et al. 2011); therefore, it is likely that sulfur in precipitation influences the δ³⁴S values of aqueous in more southerly lakes. Differences in aqueous δ³⁴S values may also be related to differences in bedrock geology among regions (Hitchon and Krouse 1972).

In three of our remote study lakes, mean δ³⁴S values in sediments were more positive than those measured in the water column (Fig. 1). However, in the wastewater-impacted Meretta and Resolute Lakes, aqueous δ³⁴S values of sulfate were 8.4‰ and 7.9‰ more positive than sedimentary values, respectively (Fig. 1). At low temperatures (i.e., 5 °C), δ³⁴S values of sediments where sulfate reduction occurs are expected to be 8‰–14‰ lower than those of water sulfate values due to bacterial fractionation (Canfield 2001); our results may therefore indicate enhanced SRB activity in these two systems, possibly due to the historical wastewater inputs altering the bacterial communities and cycling of nutrients (Antoniades et al. 2011). Indeed, the %S values in the sediments were higher in Meretta and Resolute Lakes (0.19%–0.30%) when compared with three of the remote lakes (0.1%–0.12%), but not 9-Mile Lake, which had sediments with >0.4% (see Table S1¹). Changes to the microbial community can alter sulfur cycling and different strains of SRB fractionate S isotopes at different rates (Wing and Halevy 2014; Bradley et al. 2016). Although the anthropogenic changes to bacterial communities in Meretta have been relatively well studied, any changes to the microbial community in Resolute Lake have not. Furthermore, the historical eutrophication of Meretta Lake may have created more favorable redox conditions in its sediments, enhancing sulfate reduction (Holmer and Storkholm 2001).

An increase in SRB activity may explain why [MeHg] and % MeHg were 2–6 times higher in water samples from Meretta Lake when compared with the remote lakes and with downstream Resolute Lake (see Fig. 1 and Lescord et al. 2015b). Although [MeHg] in fish and chironomids were not higher in Meretta or Resolute Lakes compared with the remote lakes (Lescord et al. 2015b), levels of MeHg in zooplankton from Meretta were four-fold higher than in all other lakes (Table S1 of the Supplementary Material¹). This, in combination with the higher aqueous δ³⁴S values of sulfate and [MeHg], as well as more negative sediment δ³⁴S values, suggests that there could be enhanced Hg(II) methylation in Meretta Lake (Canfield 2001; Alpers et al. 2014). However, recent studies suggest that the rates and degree of fractionation of S isotopes due to SRB vary among bacterial strains/species (Bradley et al. 2016). It is important to note that aqueous sulfate concentrations are very low in these high Arctic lakes (i.e., 0.23–1.35 mg/L across lakes and highest in the unaffected Small Lake); it is, therefore, likely that other processes (e.g., sulfurization; Werne et al. 2008) are also affecting S and Hg cycling, respectively, in these systems.

Although δ¹³C values indicated differences in the habitat use of zooplankton and chironomids in our study lakes (see Lescord et al. 2015b), there was no consistent difference in δ³⁴S values between these groups (Fig. 2A), which precluded the use of isotope mixing models to assess their relative importance in the diets of Arctic char. Delta³⁴S values were also not useful in distinguishing between benthic invertebrates and zooplankton as prey for fish in temperate coastal lakes (Clayden et al. 2017). In our Arctic lakes, δ¹³C values were more negative in pelagic zooplankton than benthic invertebrates in each lake (Fig. 2B) and mixing models indicated that most or all of the carbon (≥76% across lakes) in the diet of char in these systems is derived from benthic rather than pelagic sources (Lescord et al. 2015b), agreeing with the results of Muir et al. (2005) and Gantner et al. (2010) from the same lakes, and of Swanson et al. (2010) for resident fishes from other coastal Arctic lakes. Our current results, therefore, suggest that δ¹³C is a more useful indicator of diet sources for biota from high Arctic lakes than δ³⁴S, which contrasts with results from lakes in boreal Québec, Canada (Croisetière et al. 2009) and those of estuaries in the western United States (Willacker et al. 2017). Overall, far less is known about how δ³⁴S values respond to changes in S cycling and biogeochemical conditions when compared with δ¹³C and carbon cycling. It is possible that physico-chemical factors such as temperature and dissolved oxygen concentrations (and, therefore, thermal and oxic stratification patterns), which vary greatly among temperate, boreal, and Arctic regions, influence S isotope fractionation and therefore its applications to food web delineation. Indeed, these ultra-oligotrophic lakes likely do not undergo seasonal stratification during summer months, which would change the anoxic cycling of S, Hg, and other variables when compared with more southerly systems.

Across all six Arctic lakes, average [MeHg] values in surface waters were positively related to the corresponding δ³⁴S values of aqueous (LME p = 0.022; Fig. 3), which implies that lakes with ³⁴S-enriched had higher production of MeHg due to enhanced SRB activity using more of the ³²S isotope, as discussed previously. MeHg concentrations from deep water samples, however, were not significantly related to their corresponding δ³⁴S values (p = 0.645 and 0.487, r² = 0.027 and 0.049, respectively; data not shown). Because of their oligotrophic nature and lack of riparian cover, these high Arctic lakes do not stratify during the summer months (Lescord et al. 2015b). It is possible that the relationship between δ³⁴S and [MeHg] in surface (but not deep) waters is due to the influence of spring snow- and ice-melt, which are important sources of Hg to high Arctic systems and alter concentrations and δ³⁴S values in receiving waters (Loseto et al. 2004; Mörth et al. 2008; Giesler et al. 2009). Our water samples were collected after the annual spring melt so it is possible the relationship between aqueous [MeHg] and δ³⁴S values of sulfate would change if samples were collected before ice-off. It is also unclear how photochemical reactions in snow or surface waters, which affect Hg retention and loss (Mann et al. 2014), may affect concentrations and δ³⁴S values of spring melt-water.

Total [Hg] in water showed no significant relationship with δ³⁴S values of aqueous (LME p = 0.129; Fig. 3). However, this relationship is largely influenced by Meretta Lake and when removed, the negative relationship between aqueous THg concentrations and δ³⁴S became significant (LME p = 0.040, r² = 0.506; not shown). Because the number of lakes examined in this study was low, we were not able to determine whether Meretta Lake’s influential data were indeed statistical outliers. Studies have shown that the lighter δ³⁴S of DOM complexes are reflective of a terrestrial-based carbon source and to aerobic conditions in riparian soils and limited S isotope fractionation by shoreline plants (Giesler et al. 2009). DOM complexes are important transporters of Hg in aquatic systems; therefore, terrestrial organic matter increases aqueous [THg] in boreal streams and lakes presumably due to enhanced transport from riparian areas (Teisserenc et al. 2011, 2014; Eklöf et al. 2012; Lescord et al. 2018). The negative relationship between δ³⁴S values of sulfate and aqueous THg in the current study may therefore reflect higher Hg input with greater terrestrial DOM. However, sulfur compounds, DOM, and Hg have complex biogeochemical interactions (e.g., Haitzer et al. 2002; Ravichandran 2004) and further study is warranted to understand how these processes alter δ³⁴S values, particularly in ultra-oligotrophic Arctic systems.

Relationships between [MeHg] or [THg] and δ³⁴S values in biota varied among organisms in these high Arctic lakes. Benthic chironomids had positive relationships between their tissue [MeHg] and δ³⁴S (LME p = 0.008; Fig. 3), which is in contrast with the negative relationships between [MeHg] and δ³⁴S for several benthic macroinvertebrate taxa from temperate systems (Clayden et al. 2017). However, zooplankton from these Arctic lakes showed a negative relationship between their [MeHg] and δ³⁴S values (Fig. 3D). This relationship was significant after the removal of Meretta Lake, where Daphnia sp. have been detected, which elevate MeHg in bulk zooplankton (Chételat and Amyot 2009). Daphnia sp. are only found in higher-nutrient Arctic systems and likely thrive in Meretta because of the historical anthropogenic eutrophication (Chételat and Amyot 2009). Similar to benthic chironomids, both large and small char showed positive relationships between [THg] and δ³⁴S (LME p = 0.008 and 0.027, respectively; Figs. 3E and 3F). Other fish studies showed opposing log-THg versus δ³⁴S relationships between species from the same streams (Schmitt et al. 2011), no relationships between these variables across lakes (Ethier et al. 2008; Clayden et al. 2017), or that δ³⁴S was a strong positive predictor of [THg] in fishes in estuarine wetlands (Willacker et al. 2017). In the latter study, the authors attributed the positive relationship to the influence of marine S values as well as the presumed higher rates of Hg methylation and ³⁴S enrichment due to SRB activity in impounded sites (Willacker et al. 2017). The consistent, positive relationships between water [MeHg], invertebrate [MeHg], or fish [THg] and δ³⁴S values may also be related to SRB activity and suggest linkages between Hg and S cycling in these remote high Arctic systems. Studies examining the mechanisms behind these positive relationships are warranted and sampling of additional Arctic lakes would help determine whether the positive trends between biotic [THg] or [MeHg] and δ³⁴S found in this study occur more broadly.

The relationships between log₁₀[THg or MeHg] versus δ³⁴S values of food web organisms varied among lakes (ANCOVA: Lake × δ³⁴S interaction term p = 0.002, F = 4.240), with a positive slope in North Lake, a negative slope in Meretta Lake, and no significant slopes in the other four systems (Fig. 4). These mixed relationships contrast with our previous findings in seven coastal temperate lakes, where relationships were consistently negative (Clayden et al. 2017). It is important to note that the δ³⁴S values in temperate coastal lakes were generally more positive than those presented herein, which would extend the range in values on the x-axis assessed. However, the mixed Hg–δ³⁴S relationships herein also contrast with the strong and consistent relationships between log-THg or MeHg and δ¹⁵N through these food webs (Lescord et al. 2015b) and other aquatic food webs worldwide (Lavoie et al. 2013). Few studies have related [Hg] to δ³⁴S values in biota across a range of trophic levels and, in those that have, both the direction and strength of the relationships varied, particularly among regions. Our results from Arctic systems add to this mixed body of evidence and show that the [Hg]–δ³⁴S relationships in food webs can differ between similar lake ecosystems within a relatively small geographic area (i.e., ∼12 km). It is possible that various physical and chemical factors (e.g., nutrient content and anoxic habitat) may contribute to these differences, and more work is needed to investigate these effects. Interestingly, the two lake food webs that showed the strongest relationships between [Hg] and δ³⁴S had considerably different morphology: Meretta Lake was one of the smallest (0.39 km²) and shallowest lakes (max depth = 9.2 m), whereas Resolute Lake was the largest (1.3 km²) and one of the deepest systems (14.7 m) in our dataset.

Overall, the predictive models showed that various combinations of isotope and elemental measures (including the commonly used δ¹⁵N_adj, and δ¹³C as well as δ³⁴S, %N, %C, and %S) had a significant relationship with [THg] or [MeHg] in biota from these Arctic lakes. Not surprisingly, δ¹⁵N values were a consistent positive predictor of all biotic [Hg], although this relationship was only significant in zooplankton and large char (Table 1). In both large and small char models, the %S in biotic tissue was a positive predictor of [THg] (Table 1, Fig. S2¹), which may be related to the amount of S-containing amino acids, the dominant S-containing compounds in organisms (Brosnan and Brosnan 2006), as they are strong protein-based binding sites for MeHg in tissues (Peng et al. 2016; Bradley et al. 2017). In addition, [THg] or [MeHg] concentrations were positively related to %S through the food webs of five lakes (linear models, p < 0.001–0.003, r² = 0.396–0.894; Supplementary Material¹); the sixth lake, Small Lake, showed no significant relationship between [THg] or [MeHg] and %S (p = 0.130, r² = 0.162; Fig. S4¹). Interestingly, δ¹⁵N values were positively related to %S measures in each of the six food webs (p < 0.001–0.003, r² = 0.375–0.660; Fig. S3¹), implying higher S-containing amino acid content with increasing trophic level. Delta³⁴S values were identified as predictors of [MeHg] or [THg] in invertebrates and small char, although their effect was weak (i.e., coefficients were small) and varied in direction (Table 1). For example, [MeHg] in chironomids and small char were positively related to their tissue δ³⁴S, whereas zooplankton showed a negative relationship between the two variables (Fig. 3). Other elemental measures were excluded from all models (i.e., %N; Table 1) or had mixed relationships between models. Delta¹³C values, for example, were only predictive of fish Hg, with higher [THg] in large char with more positive δ¹³C values; this suggests increased Hg uptake in individuals with more benthic feeding.

Overall, these results imply that [Hg] values in Arctic biota are most strongly related to an organism’s trophic position (i.e., δ¹⁵N) but that measures of S content and δ¹³C also influence [Hg] in some groups of organisms. They further demonstrate the utility of modeling trophic interactions and contaminant transfer using multiple isotopic and elemental tracers, as they simultaneously suggest an increase in %S and [Hg] with trophic level, as well as the influence of S-cycling on Hg bioaccumulation in food webs of high Arctic lakes.