Long-Term Trends in Regional Wet Mercury Deposition and Lacustrine Mercury Concentrations in Four Lakes in Voyageurs National Park

Featured Application: Long-term monitoring of mercury in precipitation and lacustrine ecosystems provides insights into ecosystem responses as a result of mercury and sulfate emissions reductions. Abstract: Although anthropogenic mercury (Hg) releases to the environment have been substantially lowered in the United States and Canada since 1990, concerns remain for contamination in ﬁsh from remote lakes and rivers where atmospheric deposition is the predominant source of mercury. How have aquatic ecosystems responded? We report on one of the longest known multimedia data sets for mercury in atmospheric deposition: aqueous total mercury (THg aq ), methylmercury (MeHg aq ), and sulfate from epilimnetic lake-water samples from four lakes in Voyageurs National Park (VNP) in northern Minnesota; and total mercury (THg) in aquatic biota from the same lakes from 2001– 2018. Wet Hg deposition at two regional Mercury Deposition Network sites (Fernberg and Marcell, Minnesota) decreased by an average of 22 percent from 1998–2018; much of the decreases occurred prior to 2009, with relatively ﬂat trends since 2009. In the four VNP lakes, epilimnetic MeHg aq concentrations declined by an average of 44 percent and THg aq by an average of 27 percent. For the three lakes with long-term biomonitoring, temporal patterns in biotic THg concentrations were similar to patterns in MeHg aq concentrations; however, biotic THg concentrations declined signiﬁcantly in only one lake. Epilimnetic MeHg aq may be responding both to a decline in atmospheric Hg deposition as well as a decline in sulfate deposition, which is an important driver of mercury methylation in the environment. Results from this case study suggest that regional- to continental-scale decreases in both mercury and sulfate emissions have beneﬁtted aquatic resources, even in the face of global increases in mercury emissions.


Introduction
Human activities have considerably increased the amount of mercury in the atmosphere and deposited into aquatic ecosystems [1]. Elevated mercury levels in fish have resulted in widespread fish-consumption advisories across the United States (U.S.), North America, and globally to protect human health [2], and also pose ecotoxicological risk to piscivorous wildlife [3]. High fish-mercury levels can occur anywhere, including remote, relatively pristine ecosystems where the predominant source of mercury is atmospheric deposition. Table 1. Trends in emissions of mercury in the U.S. [9] and in the Canadian provinces of Manitoba and Ontario [10], 1990-2017. [U.S. data were reported in short tons; Canadian data were reported in kilograms [kg]; all data have been converted to Mg yr −1 ; percentage change is from 1990 to 2017.]. Perhaps equally beneficial, from a methylmercury production and bioaccumulation perspective, efforts to control acid precipitation during the 1970s and 1980s led to large decreases in emission and deposition of sulfur oxides across the U.S. since 1970 [11]. Consistent with this continental-scale trend, sulfate deposition within the Voyageurs National Park region has declined substantially. Data from the National Atmospheric Deposition Program's (NADP) site MN16 (Marcell, Minnesota) show annual wet sulfate deposition declining from about 10 kg ha −1 in 1980, to about 6 kg ha −1 in 2000, to about 2.4 kg ha −1 in 2018 ( Figure S1). Drevnick et al. [12] have previously attributed declining fish-mercury levels in lakes at Isle Royale National Park (Lake Superior, USA) to declines in sulfate loading.

Year
Although North American mercury emissions have declined sharply (Table 1 and [13]), and emissions have declined in several other regions as well [13], industrial development and associated increases in mercury emissions-particularly in Asia-have led to a net global increase in anthropogenic mercury emissions of 1.8 percent per year [14]. The relative importance of North American versus global emissions to local deposition in northern Minnesota is unclear. A recent modeling study found that North American emission reductions have a more pronounced effect on wet mercury deposition reductions in the eastern U.S., where there was a greater concentration of sources, than in comparatively remote northern Minnesota [15], but that effort did not focus on ecosystem responses.
A critical question arises: how have lacustrine ecosystems in midcontinental North America responded to regional and North American mercury emissions versus global emissions? This question is particularly important because the Minamata Convention on Mercury [16], a legally-binding international treaty to reduce mercury use and releases to the environment, is being implemented by 128 signatory nations. Thus, assessing ecosystem responses to contemporary mercury reductions may be informative in assessing treaty implementation.
A few other studies have attempted to answer this question. There is mixed evidence for ecosystems responses to changing emissions trajectories, including some literature showing declining fish-mercury levels in response to regional emission reductions [17] and other studies showing long-term (1970s-2000s) declines in fish-mercury concentrations followed by a leveling off or increase in concentrations in recent years [18,19]. Mercury in the feathers of bald eagle (Haliaeetus leucocephalus) nestlings have shown similar temporal trends within the region [20].
This paper addresses the question: how have aquatic ecosystems responded? Here, we report on one of the longest known paired data sets for mercury in atmospheric deposition; mercury, methylmercury, and key other chemical and physical parameters in several lacustrine systems; and total mercury in aquatic biota from the same lacustrine systems. This paper reports an updated analysis of trends in mercury deposition at two northern Minnesota NADP Mercury Deposition Network (MDN) sites, as well as trends in methylmercury and total mercury concentrations in epilimnetic lake water and in biota from four lakes in a minimally disturbed area within Voyageurs National Park, also in northern Minnesota. Through several earlier studies [21][22][23][24], and an ongoing collaboration between the U.S. Geological Survey (USGS) and the National Park Service (NPS), these lakes have been sampled for mercury since 2000 or 2001 (depending on the lake), making these among the longest running data sets that pair total mercury and methylmercury in lake water with biotic mercury. Given the paucity of long-term data sets on aqueous mercury in undisturbed ecosystems, the sensitivity of circumneutral, low-ionic-strength aquatic ecosystems to the effects of mercury and sulfate deposition [7,25], and widespread consumption advisories in this region, the Voyageurs National Park data set provides a useful case study to monitor ecosystem responses to changes in atmospheric inputs. Trends in regional wet deposition, epilimnetic mercury, and mercury in age-1 yellow perch (Perca flavescens) were previously reported through 2012 [24]. This paper examines trends through 2018. In the current paper, the longer time period necessitated a change in trend analysis to account for nonlinearity in the data; also, given a change in organisms collected for biomonitoring, we applied a statistical relationship between mercury in dragonfly larvae and in age-1 yellow perch in order to extend the perch record.

Study Area
The study area and methods have been previously described [24]. In brief, the four study lakes lie within Voyageurs National Park, a park covering 883 km 2 in northeastern Minnesota ( Figure 1; Table 2). The four lakes are drainage lakes, three of them fed by a small stream discharging from an upstream lake; the fourth (Ryan) fed only by a small headwaters stream. Although the upstream lakes may exert some influence on chemistry of the study lakes, given the small size of the inflowing streams, and relatively long water renewal times (0.6-1.2 year), the lake chemistry likely is governed more by catchment and deposition to the lake surface than by inflows from upstream lakes. Land cover within the park is largely boreal forest, with thin soils, and outcrops of Precambrian bedrock common throughout the park [26]. The lakes have been sampled for methylmercury and total mercury in lake water since 2000 (Shoepack Lake) or 2001 (Brown, Peary, and Ryan lakes) as part of earlier studies [21][22][23][24].
No monitoring of mercury deposition occurs within the park's borders; the two nearest NADP/MDN sites (http://nadp.slh.wisc.edu/mdn/, accessed on 1 February 2021) are used to characterize wet deposition of mercury: sites MN16 (Marcell, Minnesota), approximately 120 km south-southwest of Voyageurs National Park and MN18 (Fernberg, Minnesota), approximately 120 km southeast of the park. Due to the lack of large urban or industrial centers, relatively flat topography (hence lack of orographic effects), and variable wind direction in the region, data from these two monitoring sites are expected Appl. Sci. 2021, 11, 1879 4 of 20 to be representative of mercury and sulfate deposition in the region. We examined trends in wet Hg deposition from 1998-2018; trends in aqueous methymercury (MeHg aq ) and total mercury (THg aq ) from 2000 (Shoepack Lake) or 2001 (Brown, Ryan, and Peary lakes) through 2018; and trends in biotic THg as described below.
No monitoring of mercury deposition occurs within the park's borders; the two near-141 est NADP/MDN sites (http://nadp.slh.wisc.edu/mdn/) are used to characterize wet depo-142 sition of mercury: sites MN16 (Marcell, Minnesota), approximately 120 km south-south-143 west of Voyageurs National Park and MN18 (Fernberg, Minnesota), approximately 120 144 km southeast of the park. Due to the lack of large urban or industrial centers, relatively 145 flat topography (hence lack of orographic effects), and variable wind direction in the re-146 gion, data from these two monitoring sites are expected to be representative of mercury 147 and sulfate deposition in the region. We examined trends in wet Hg deposition from 1998-148 2018; trends in aqueous methymercury (MeHgaq) and total mercury (THgaq) from 2000 149 (Shoepack Lake) or 2001 (Brown, Ryan, and Peary lakes) through 2018; and trends in biotic 150 THg as described below.  Table 2. Selected lake characteristics [22,26]; and selected water-column measurements from the approximate centroid of 155 each lake [24] [27]. Briefly, weekly composited 160 precipitation samples are collected and analyzed using modifications of U.S. Environmen-161 tal Protection Agency (U.S.EPA) Methods 1669 [28] and 1631 [29]. Precipitation samples 162  Table 2. Selected lake characteristics [22,26]; and selected water-column measurements from the approximate centroid of each lake [24]. [W.A., watershed area; % Wetlands, wetlands as a percentage of watershed area; Renewal time = mean hydraulic residence time, in years; ANC, acid-neutralizing capacity; µeq L −1 , microequivalents per liter; TOC, total organic carbon; mg L −1 , milligrams per liter].

Wet Deposition
Detailed methods are described by Prestbo et al. [27]. Briefly, weekly composited precipitation samples are collected and analyzed using modifications of U.S. Environmental Protection Agency (U.S.EPA) Methods 1669 [28] and 1631 [29]. Precipitation samples were analyzed for total mercury concentration (Hg precip ) at Frontier Global Sciences in Seattle, Washington. Data included herein met the NADP's acceptance criteria, [27,30] although some further screening of apparent outlier Hg precip values was done as described below.
Inspection of weekly MDN data revealed a limited number of outliers in the weekly data, some of which were suspected as extreme outliers, i.e. unusually high or low, and which may impart an artefact into the calculations of annual wet Hg deposition rate and precipitation-weighted mean concentrations. We used the following approach to identify extreme outliers prior to wet deposition calculations: Weekly data for the study sites were downloaded from the NADP/MDN website. Weekly Hg precip and precipitation volumes were analyzed using a multiple-regression Appl. Sci. 2021, 11, 1879 5 of 20 approach developed previously [31], using the REG procedure in SAS software (v. 9.4. SAS Institute, Carey, NC, USA). Hg precip and precipitation depth data were log-transformed to remove skewness and heteroscedasticity of residuals. A regression model was developed to account for the exogenous effects of weekly precipitation depth; season (Fourier terms); and time (T, in years) (Equation (1)).
log[Hg precip ] = β 0 + β 1 log(Precip) + β 2 sin(2πT) + β 3 cos(2πT) + β 4 T + ε where Hg precip is mercury concentration in weekly precipitation samples in ng L −1 ; Precip is the precipitation depth in mm; sin and cos are the sine and cosine functions; T is time in decimal years. Additional Fourier terms (sine and cosine of 4πT) were significant at many MDN sites across North America [31], but were not significant for the current study sites, thus were not included in the analysis herein. The residual, ε, is the unexplained variation in log[Hg precip ], accounting for exogenous effects of precipitation amount, seasonality, and time. Inspection of the data indicated the presence of outliers. Therefore, two passes of the regression analysis were performed. After the first pass, Hg precip in weekly samples was set to missing if deemed an outlier according to the following rules: Rule 1: absolute value of studentized residual from multiple-regression analysis (Equation (1)) greater than 4, regardless of precipitation amount. This removed extreme outliers for weekly samples with less than 10 mm of precipitation; low-volume events tend to be noisier. From this analysis, four outlier results were removed for MN16 and three for MN18.
Rule 2: for samples with precipitation volume >10 mm, values were set to missing if the absolute value of studentized residual from multiple-regression analysis (Equation (1)) was greater than 3. Using this rule, one Hg precip result was set to missing for MN16 and two for MN18.
Mercury concentrations for the screened samples were deemed extreme outliers for the precipitation amount of the event and were likely an artefact rather than a realistic representation of Hg precip in precipitation. One outlier in particular had an abnormally high Hg precip for a relatively large precipitation event and occurred late in the study period (site MN18; 5 July 2017; reported Hg precip = 89.59 ng L −1 ; Precip = 27.686 mm). The resultant calculated weekly Hg deposition was 2.48 µg m −2 , which is about 25 percent of the annual wet deposition reported by NADP/MDN at this site and would have influenced the trend analysis were it not removed. A table of removed outliers is provided (Supplementary Materials, Table S1).
The above-described screening process retained weekly observations (with precipitation volume retained), but outlier values of Hg precip were set to missing. After removal of outliers, regression analysis was applied (again using Equation (1)), and missing log[Hg precip ] values were predicted by the regression. Predicted Hg precip was calculated by exponentiating regression-predicted log[Hg precip ] values, and multiplying by the mean of the exponentiated residuals to correct for back-transformation bias (Duan smearing estimator, as described in Helsel et al., pp. 256-257 [32]). Weekly Hg precip values, then, consisted of MDN-reported values in nearly all cases, but regression-predicted values where weekly Hg precip was missing (including observations where outliers were screened as described above).
Annual wet Hg deposition was then calculated as the sum of weekly deposition rates (concentration times precipitation depth); and precipitation-weighted mean Hg precip were calculated as annual deposition rate divided by annual precipitation depth. This methodology differs from NADP/MDN network methodology. NADP/MDN effectively estimates missing Hg precip as the precipitation-weighted annual mean concentration, whereas our methodology estimates missing Hg precip using regression-based predictions based on the precipitation volume and season associated with each weekly sample. For the purposes of this paper, regression-based estimates were preferred because, given the dependence of Hg precip concentrations on precipitation depth, precipitation-weighted mean concentrations may underestimate concentrations for small events and overestimate concentrations for large events (which are more important in determining annual loads).
Owing to modest nonlinearity in the time trends, trends in the annual wet Hg deposition rate and in precipitation-weighted annual mean Hg concentration were determined by locally weighted regression, after dealing with extreme outlier Hg concentrations as described above. Locally weighted regression used the LOESS procedure in SAS software, allowing the procedure to select the smoothing parameter by Akaike information criterion algorithm.

Lake Water
Lake water sampling and analytical methods have been described previously [24,33]. Unfiltered, epilimnetic lake water was sampled two to three times (typically three) per year, between May and September, during most years of the study period. Field crews generally sampled the upper ca. 5 cm of water out of a canoe heading upwind, in the approximate center of the lake, using trace-metal clean sampling protocols described in [21]. Water was collected in a pre-cleaned Teflon ® FEP (fluorinated ethylene propylene) bottle, preserved by addition of HCl (to a final normality of ca. 0.02 N), and shipped to the USGS Mercury Research Laboratory (MRL) in Middleton, Wisconsin, for analysis. Detailed descriptions of all the analytical procedures used by the USGS MRL are available at the following website: https://wi.water.usgs.gov/mercury-lab/research (accessed on 1 February 2021) and the descriptions are summarized below. Total Hg determinations in lake water (THg aq ) were determined by U.S.EPA Method 1631 [29]. The USGS MRL typically achieves a daily detection limit (DDL) for total mercury analytical runs of about 0.06 pM and the precision, measured as the relative percent difference (RPD) between analytical duplicates averages 10 percent. Methylmercury (MeHg aq ) samples were analyzed by U.S.EPA Method 1630 ( [34,35] with the added advancement (starting in 2006) of including isotope dilution by adding a small known amount (about 30 picograms) of isotopically labeled methylmercury (Me 200 Hg) to each sample, which allows for more accurate measures of sample recovery rates.
Water samples for other constituents were also collected from the same lake locations with a 2 m long, 3.2 cm inner diameter PVC tube that integrates a 2 L sample from the upper 2 m of the water column. Samples collected once per summer, typically in July, were analyzed for several major ions including sulfate by White Water Associates, Inc. (Amasa, Michigan) (2006-2013) and CT Laboratories (Baraboo, Wisconsin) (2014-2018). Sample processing, handling, and quality assurance and quality control procedures are described in Elias et al. [36].
For the 2001-2012 time period [24], trends in MeHg aq conformed reasonably to a linear-regression model for three of the four lakes. However, inspection of data for the longer time (through 2018) revealed nonlinearity at some sites. Therefore, we used locally weighted regression using the LOESS procedure in SAS, again allowing the procedure to select the smoothing parameter by the default Akaike Information Criterion algorithm within the LOESS procedure.
To calculate the percent change in concentrations from 2001-2018, we used LOESSpredicted concentrations for arbitrary dates approximately at the beginning and end of the period of data collection (1 July 2001 and 1 July 2018). No slope parameter (and hence, no p-value for significance) is calculated in LOESS. Relevance of the trends can be assessed by examining the magnitude of change viewed along with the variability of the data; and by examining whether the 95-percent confidence interval for the LOESS smooth at the end of the time period includes, or does not include, the LOESS-predicted value for the start of the time period. Lake water chemistry and lake level data from 2000-2007 are available from the USGS's National Water Information System web retrieval (https://doi.org/10.5066/F7P55KJN, accessed on 1 February 2021), for the following USGS site identification numbers: Brown Lake (483059092474501); Peary Lake (483129092462001); Ryan Lake (483109092422601); and Shoepack Lake (482951092531601). Starting in 2006, sample collection and data archiving was led by the NPS, and data are available for retrieval from the National Water Quality Monitoring Council's (NWQMC) water-quality data portal (https://www.waterqualitydata.us/portal/, accessed on 1 February 2021), using the organization ID of 11NPSWRD_WQX, and the following site identification numbers: Brown Lake (11NPSWRD_WQX-VOYA_12); Peary Lake (11NPSWRD_WQX-VOYA_14); Ryan Lake (11NPSWRD_WQX-VOYA_17); and Shoepack Lake (11NPSWRD_WQX-VOYA_05). Data for mercury in yellow perch, and dragonfly larvae from 2008-2012 are also available at the NWQMC's water-quality portal for the same sites. Dragonfly larvae (Odonata, Anisoptera) data from 2014-2018 are available within a USGS data release [37].

Lake Levels
In 2006, the NPS established reference points on the shore of each lake. After establishing these, lake levels were determined relative to an arbitrary datum at each lake's reference point. Water-level anomaly was calculated as the difference between water level on a given sampling date and the initial water level relative to local datum.

Mercury in Lake Biota
Details of sampling for age-1 yellow perch and dragonfly larvae have been described elsewhere [24,38,39]. In brief, both yellow perch and dragonfly larvae were sampled annually from each lake during spring. Because total mercury concentrations (THg) can vary among families [38], we normalized THg in dragonfly larvae to those of a single family (Aeshnidae) following Eagles-Smith et al. [38]. This ensures a consistent unit of dragonfly larvae THg for each site year. Fish sampling occurred from 2000 to 2012, whereas dragonfly larvae were sampled from 2009 to 2018. To facilitate temporal comparisons across the entire study period we first examined the relation in THg concentrations between paired samples of yellow perch and Aeshnid-equivalent dragonfly larvae from 14 lakes in national parks in the western Great Lakes region because previous findings have shown them to be correlated [40], using linear regression of the geometric mean THg concentrations of each taxa where they were collected together. This analysis indicated that dragonfly larvae THg concentrations well correlated with those in yellow perch (See Results); therefore, we used the linear regression equation to estimate yellow perch THg concentrations for years where yellow perch were not sampled.
As with aqueous MeHg concentrations, examination of the temporal trends in yellow perch THg revealed nonlinearity in the three lakes with a complete temporal data set (Brown Lake, Ryan Lake, and Parry Lake). Therefore, we similarly used the LOESS procedure as described above to estimate the change in biotic Hg concentrations over time, though the lower data density necessitated a higher degree of smoothing than with the higher resolution MeHg aq sampling.

Wet Deposition Trends
Rates of both annual wet Hg deposition and precipitation-weighted mean Hg concentrations at both MN16 and MN18 decreased over the time period 1998-2018 ( Figure 2; Table 3); much of the decreases occurred prior to 2009, with relatively flat trends since 2009. Of interest from an ecosystem perspective is the change in wet-depositional loading. For reasons previously noted, removal of outliers, then use of regression-predicted weekly Hg precip values to calculate weekly and annual wet Hg deposition rates, was preferable to using NADP/MDN-reported annual deposition rates. Expected values, from locally weighted regression analysis, of annual wet deposition rates decreased by 20 and 24 percent, respectively for MN16 and MN18, over the 1998-2018 period. These percentage declines are smaller than those reported for 1998-2012 [24]; the change in magnitude may be due to a different analytical approach and effect of removing influential outliers, as well as a flattening of the trends starting around 2009. The flatter trends since~2009 may reflect a leveling-off of emissions-reductions in the U.S. and Canada over the last decade; i.e., large reductions in mercury emissions occur early in this time period, with more modest reductions in recent years (Table 1). Table 3); much of the decreases occurred prior to 2009, with relatively flat trends since 321 2009. Of interest from an ecosystem perspective is the change in wet-depositional loading. 322 For reasons previously noted, removal of outliers, then use of regression-predicted weekly 323 Hgprecip values to calculate weekly and annual wet Hg deposition rates, was preferable to 324 using NADP/MDN-reported annual deposition rates. Expected values, from locally 325 weighted regression analysis, of annual wet deposition rates decreased by 20 and 24 per-326 cent, respectively for MN16 and MN18, over the 1998-2018 period. These percentage de-327 clines are smaller than those reported for 1998-2012 [24]; the change in magnitude may be 328 due to a different analytical approach and effect of removing influential outliers, as well 329 as a flattening of the trends starting around 2009. The flatter trends since ~2009 may reflect 330 a leveling-off of emissions-reductions in the U.S. and Canada over the last decade; i.e., 331 large reductions in mercury emissions occur early in this time period, with more modest 332 reductions in recent years (Table 1) Because precipitation amount is used in the calculation of wet Hg deposition rate, a trend in precipitation could drive a trend in Hg deposition. Linear regression of precipitation depth versus time shows no significant trend for MN16 (p = 0.80) and a weak positive trend at MN18 (p = 0.10) (see Figure S2: precipitation volume trend plots). The sites display considerable interannual variability in total precipitation depth (ranges: 554-908 and 508-832 mm yr −1 for MN16 and MN18, respectively), driving interannual variability in wet Hg deposition. The weak, positive trend in precipitation amount at MN18 likely drove the relatively larger decline in precipitation-weighted mean concentrations at that site, compared to MN16, as larger precipitation events tend to have lower Hg precip .
Given the lack of a significant trend in precipitation volume, we conclude that the observed overall trend of declining wet Hg deposition rate is likely driven by reductions of mercury emissions in North America and not trends in precipitation. In addition, the observed declines in mercury deposition are synchronous with known declines in North American mercury emissions since 1990 [13,41], although global emissions have been comparatively constant [13]. A recent modeling study indicated that in northern Minnesota, emission reductions in North America are roughly equally important in comparison to emission reductions in the rest of the world in determining wet Hg deposition trends [15].

Trends in Epilimnetic Methylmercury, Total Mercury, and Sulfate Concentrations
Both methylmercury and total mercury concentrations declined in epilimnetic lake water over the 2001-18 period, although the declines at some lakes were small, relative to variability. For Brown Lake, a high methylmercury outlier previously identified [24] was omitted from trend analysis. The overall decline in MeHg aq for Brown Lake is modest (32%; Table 4), but noteworthy for two reasons. First, the previously reported trend for Brown Lake [24] was positive, but weak; and second, concentrations have declined sharply since peaking around 2010, about the end of the time frame for the previous trend analysis. Peary and Ryan Lakes had the largest declines in MeHg aq (Table 4), with much of the change occurring in the first few years of record. The comparatively large declines in these two lakes is consistent with earlier findings [24]. Shoepack Lake also exhibited a modest decline in MeHg aq . The magnitude of the declines in MeHg aq for Brown and Shoepack Lakes is modest, in relation to the variability in concentrations at these two lakes; in addition, the 95-percent confidence intervals for the start and end of the period of study overlap for these two lakes suggesting that the decline is not significant (Figure 3).  Epilimnetic sulfate concentrations decreased in each lake over the study period (Ta-385 ble 4; Figure 5). The trends were nonlinear, revealing midtime series peak concentrations 386 around 2007 for Brown, Ryan, and Shoepack Lakes, and somewhat later (~2010) for Peary 387 Lake. The modest trends reported here (mean decrease of 45%) follow much larger de-388 creases for these lakes from the 1980s to 2000, as reported by Kallemeyn et al. [26], based 389 on a 1980s lake survey reported by Payne [42].  Aqueous total mercury declined modestly in Brown, Peary, and Shoepack Lakes, with an overlap of the 95-percent confidence intervals for the start and end period of the study (again, indicating perhaps a lack of statistical significance) (Figure 4). The 47% decline in THg aq in Ryan Lake appears to be a significant decline with clear separation of confidence intervals in the beginning versus end of the study period. Epilimnetic sulfate concentrations decreased in each lake over the study period (Table 4; Figure 5). The trends were nonlinear, revealing midtime series peak concentrations around 2007 for Brown, Ryan, and Shoepack Lakes, and somewhat later (~2010) for Peary Lake. The modest trends reported here (mean decrease of 45%) follow much larger decreases for these lakes from the 1980s to 2000, as reported by Kallemeyn et al. [26], based on a 1980s lake survey reported by Payne [42].

Lake Level Fluctuations
MeHg aq correlates modestly (R 2 = 0.36) with lake-level anomaly for Brown Lake ( Figure 6). However, none of the other lakes in this study showed significant correlations between MeHg aq and lake-level anomaly (not shown). MeHgaq correlates modestly (R 2 = 0.36) with lake-level anomaly for Brown Lake (Fig-414  ure 6). However, none of the other lakes in this study showed significant correlations be-415 tween MeHgaq and lake-level anomaly (not shown).

421
Total mercury concentrations in age-1 yellow perch were well correlated with Aesh-422 nid-equivalent dragonfly larvae THg concentrations (R 2 = 0.66, p < 0.0001, N = 40; Figure 423 7), facilitating converting dragonfly THg concentrations to those of yellow perch for years 424 when fish were not sampled. Yellow perch data for all four lakes were collected through 425 2012 as summarized previously [24]. After yellow perch collections ceased, only Brown, 426

Mercury Trends in Lake Biota
Total mercury concentrations in age-1 yellow perch were well correlated with Aeshnidequivalent dragonfly larvae THg concentrations (R 2 = 0.66, p < 0.0001, N = 40; Figure 7), facilitating converting dragonfly THg concentrations to those of yellow perch for years when fish were not sampled. Yellow perch data for all four lakes were collected through 2012 as summarized previously [24]. After yellow perch collections ceased, only Brown, Peary, and Ryan Lakes were sampled for dragonfly larvae, so only those three lakes are considered here. Peary, and Ryan Lakes were sampled for dragonfly larvae, so only those three lakes are 427 considered here. 428 For the three lakes with long-term biomonitoring, temporal patterns in biotic THg 429 concentrations were similar to patterns in MeHgaq concentrations (Figures 3 and 8); how-430 ever, biotic THg concentrations declined in only Peary Lake. Expected values for yellow 431 perch THg for Brown Lake increased by 4.6% between 2000 and 2018, but the 95-percent 432 confidence intervals overlapped between those years indicating that the difference is not 433 significant; similar to MeHgaq, there was a substantial 54% increase in THg concentrations 434 in yellow perch between 2000 and 2010, followed by a 46% decrease between 2010 and 435 2018. As with MeHgaq, Peary Lake had the greatest decline in biotic THg (45%), which was 436 primarily driven by the 31% decrease between 2000 and 2010. Ryan Lake showed an initial 437 decline of biotic THg until about 2010, followed by an increase; yellow perch THg concen-438 trations in Ryan Lake increased by 5% between 2000 and 2018, similar in magnitude to 439 Brown Lake, and not apparently significant. However, in contrast with Brown Lake, at 440 Ryan Lake there was a substantial decline (38%) between 2000 and 2010, followed by a 441 69% increase between 2010 and 2018-a considerably sharper increase than MeHgaq con-442 centrations during the same time period. For the three lakes with long-term biomonitoring, temporal patterns in biotic THg concentrations were similar to patterns in MeHg aq concentrations (Figures 3 and 8); how-ever, biotic THg concentrations declined in only Peary Lake. Expected values for yellow perch THg for Brown Lake increased by 4.6% between 2000 and 2018, but the 95-percent confidence intervals overlapped between those years indicating that the difference is not significant; similar to MeHg aq , there was a substantial 54% increase in THg concentrations in yellow perch between 2000 and 2010, followed by a 46% decrease between 2010 and 2018. As with MeHg aq , Peary Lake had the greatest decline in biotic THg (45%), which was primarily driven by the 31% decrease between 2000 and 2010. Ryan Lake showed an initial decline of biotic THg until about 2010, followed by an increase; yellow perch THg concentrations in Ryan Lake increased by 5% between 2000 and 2018, similar in magnitude to Brown Lake, and not apparently significant. However, in contrast with Brown Lake, at Ryan Lake there was a substantial decline (38%) between 2000 and 2010, followed by a 69% increase between 2010 and 2018-a considerably sharper increase than MeHg aq concentrations during the same time period.

Discussion
Wet Hg deposition at two regional MDN sites (Fernberg and Marcell, Minnesota) declined by an average of 22 percent from 1998-2018, with much of the decline occurring prior to 2010. In the four lakes, epilimnetic MeHg aq concentrations declined by an average of 44 percent and THg aq by an average of 27 percent. Although the magnitude of trend in some lakes was small, it is noteworthy that for all the lakes both MeHg aq and THg aq show declines for the 2001-2018 time period, including the latter part of that period when wet Hg deposition rates leveled off, suggesting a lag related to watershed inputs. Epilimnetic MeHg aq may be responding both to a decline in atmospheric Hg deposition as well as a decline in sulfate deposition, which is an important driver of mercury methylation in the environment. The long-term reduction in epilimnetic sulfate concentrations in the lakes also reflects declines in sulfate deposition, as has been observed elsewhere [43]. This observation is a good example of the importance of collecting data on other known key factors (for example, sulfate) that control mercury cycling in the environment when the goal is to accurately attribute the drivers of change.
Environmental mercury data sets that include long-term monitoring of multiple media (atmosphere, surface water, and biota) in a relatively small area are exceedingly rare. As such, the opportunities to assess baselines and trends in mercury levels in aquatic ecosystems, especially in the lead-up to expected globally driven emissions changes from the Minamata Treaty [16], are likewise uncommon.
Previously, it was hypothesized that inflowing water from a lake upstream from Brown Lake (Oslo Lake), which yielded relatively high concentrations of MeHg aq in a 2001-2002 survey of small lakes in Voyageurs National Park [21], could be responsible for the increase in MeHg aq in Brown Lake during the 2001-2012 time period. The apparent trend reversal, i.e., the decline in MeHg aq in Brown Lake that occurred starting around 2010, coupled with the correlation between lake level and MeHg aq , supports the hypothesis.
Higher observed MeHg aq concentrations coincident in time with higher lake levels is consistent with the generally held conceptual understanding from the mercury literature that wetter conditions and cyclical inundation and draining of low-lying areas (such as wetlands) leads to increased MeHg aq production within a lake's watershed, irrespective of loading from upstream lakes. However, the remaining lakes in this study showed no correlation between MeHg aq and lake level. This null finding is in contrast to the coherence of water level and MeHg in water and fish observed for lakes in northern Wisconsin [44,45]. However, it should be noted that the ecological setting in northern Wisconsin is quite different from the Voyageurs National Park region, especially in regard to hydrology. The Northern Highlands of northern Wisconsin are characterized by high permeability due to deep outwash sands and gravel that yield poorly integrated surface drainages. Our study area, in contrast, is more of a classical boreal system with shallow soils overlying bedrock and highly integrated flow systems. As such, the lack of concurrence between findings of Watras et al. [44,45] and our study is not surprising.
Trends in fish-tissue THg concentrations moderately tracked MeHg aq or THg aq for Brown and Peary Lakes but not Ryan Lake. Whereas MeHg aq concentrations often correspond to biological mercury uptake and accumulation in many water bodies, there can be substantial variability in the efficiency of transfer into and through food webs due to the context dependence of site-specific bio-geochemical and ecological drivers. Also, whereas MeHg aq or THg aq is an instantaneous measure of conditions, biological tissues integrate exposure over much longer time periods, including the Odonates that are generally several years old. This disparity can complicate interpretations of the effectiveness of decreasing mercury emissions and deposition. In addition, for boreal-like settings, the connectivity to terrestrial soils and their legacy accumulation of decades of mercury deposition is well understood; however, the internal time lags of how long this large pool of mercury will continue to yield meaningful amounts of mercury to downstream aquatic ecosystems remains unknown. This finding does not imply that declining mercury emissions and deposition (and subsequent MeHg aq production) offer limited benefits for mercury risk reduction. Instead, it emphasizes the need to interpret long-term environmental mercury data sets in the context of a complex set of pathways and processes that control mercury cycling in the environment, and the need for multimedia (deposition, water, and biota) and multiconstituent data (beyond just mercury and methylmercury measurements) for effective trend analysis for mercury.
Results from this case study suggest that regional-to continental-scale decreases in both mercury and sulfate emissions have benefitted aquatic resources, even in the face of global increases in mercury emissions. The reductions in atmospheric pollutant loading may be of considerable benefit to human and ecosystem health, considering that mercurybased fish-consumption advisories are in place for all lakes of Voyageurs National Park and many lakes in the region, and northern pike (Esox lucius) mercury levels in park lakes have exceeded thresholds for detrimental effects to fish reproduction [46,47].
A number of MDN sites across North America had substantial declines in wet mercury deposition from the late 1990s through early 2000s, followed by a leveling-off and in some cases increase in Hg deposition starting around 2010 [31]. The two northern Minnesota sites appear to fit this broader geographic pattern. Trend analysis by locally weighted regression showed a relatively sharp decline in the period from 1998 to about 2010, followed by a leveling-off of deposition rates. Although the MDN has a data review and quality assurance program in place, the data-screening procedure employed herein identified a small number of extreme outlier mercury concentrations that appeared unreasonable. Because extreme outliers for an individual sample can bias annual wet-deposition calculations, our screening procedure (or similar ones) warrants further consideration.
The larger declines in epilimnetic MeHg aq , compared to epilimnetic THg aq is likely driven by both decline in wet Hg deposition (and thus in-lake THg aq ), as well as declines in sulfate deposition. As noted previously [24], in response to emission controls related to the Clean Air Act, sulfate deposition has declined dramatically in northern Minnesota, as well as more broadly across North America [11,48]. Other research has shown that adding sulfate to wetlands greatly increases methylmercury production [7], whereas decreased sulfate loading results in decreased methylmercury production [49].
Owing to long-term atmospheric deposition of anthropogenic mercury and sulfate, lakes in Voyageurs National Park, and regionally, surely have elevated methylmercury levels in both water and biota, relative to pre-industrial conditions. It is encouraging, however, that declines in anthropogenic mercury and sulfur emissions have translated to declines in wet mercury and sulfate deposition, which in turn appear to have resulted in declines in methylmercury contamination in lake ecosystems within the park. The relatively large MeHg aq declines, in comparison to declines in THg aq , are consistent with the notion that MeHg aq levels are influenced by both anthropogenic mercury as well as anthropogenic sulfate deposition. It is worth noting that the response of lake ecosystems to decreased mercury inputs is expected to include both a rapid component owing to direct deposition to the lake surface and a slow component driven as previously deposited mercury slowly re-equilibrates from wetland and upland soils [50].
Lastly, as emphasized previously, there are relatively few published long-term, multimedia data sets that include atmospheric mercury-deposition monitoring coupled with methylmercury and total mercury in lake water and mercury in lake biota. This is particularly important for undisturbed settings where methylmercury production and bioaccumulation are largely governed by natural processing of atmospheric pollutant loads. Watras et al. [44] reported trends for aqueous total mercury and methylmercury and biotic mercury for two lakes in northern Wisconsin (Little Rock Lake, 1988-2017 and Trout Bog, 1999-2017) that are relatively close (ca. 275 km southeast of VNP), yet the two studies yielded trend analyses that are notably different. This variability in temporal trends is consistent with the overarching influence that within-lake and watershed bio-geochemistry can have on mercury transport and methylmercury production, potentially decoupling them from trends in mercury loading. This highlights the importance of considering the context of each ecosystem and supports the notion that recovery from many decades of sustained mercury emissions is unlikely to be a linear process. Data sets like the one used in this study, while rare, will serve as critically important baselines for executing effectiveness evaluations associated with the post-Minamata-Treaty implementation. Although more extensive networks have been proposed to monitor ecosystem responses to controls on anthropogenic mercury emissions [51], in the absence of such programs leading up to the global change expected from the Minamata Treaty, researchers might better coordinate small-scale, long-term research efforts so that broader-scale assessments can be made.
Author Contributions: M.E.B. analyzed precipitation and lake-water chemistry data; C.A.E.-S. analyzed the biotic data. J.F.D. contributed to water chemistry methods and quality assurance of the water chemistry data. All authors contributed to preparation of the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: Preparation of this manuscript was funded by the National Park Service (NPS). National Atmospheric Deposition Network/Mercury Deposition Network sites were supported by the U.S. Forest Service, U.S. Environmental Protection Agency, and the Minnesota Pollution Control Agency (MPCA). Funding earlier research that contributed data to the current paper is detailed in [24]. The University of Wisconsin at LaCrosse's collection of some of the biotic data was funded by the Great Lakes Restoration Initiative. Partial financial support for lake water and biotic chemistry data was provided by the MPCA and the U.