Environ. Sci. Technol. 2007, 41, 4568-4573
Time Trends of Methylmercury in Walleye in Northern Wisconsin: A Hierarchical Bayesian Analysis E R I C R . M A D S E N * ,† A N D H A L S . S T E R N ‡ Great Lakes Indian Fish and Wildlife Commission, P.O. Box 9, Odanah, Wisconsin 54861, and Statistics Department, University of California-Irvine, Irvine, California 92697
Methylmercury data from walleye fillets collected by multiple agencies from northern Wisconsin lakes from 1982 to 2005 were examined for regional time trends. Hierarchical Bayesian methods were used to model dependencies and provide probability statements for parameters pertaining to individual lakes and the region as a whole. A missing data mechanism allowed the sex of the fish to be included as a predictor since the sexes grow at different rates. A slight regional decrease in methylmercury of 0.60% annually was found, consistent with declining atmospheric mercury deposition. Methylmercury was estimated to have decreased in 77% of the 420 lakes from which walleye were sampled, although uncertainty regarding time trends was greater for most individual lakes than for the region as a whole. Methylmercury in walleye varied widely from lake to lake, but generally accumulated in the fish at similar rates by length after accounting for differences in sex. Slower-growing male walleye had higher methylmercury concentrations than females for a given length, and skin-on fillets were 16% lower in methylmercury than skin-off fillets.
Introduction Methylmercury (MeHg) contamination of fish has prompted numerous consumption advisories around the world. Atmospheric deposition is thought to be a major source of MeHg contamination, particularly in remote areas (1). There is evidence that atmospheric mercury deposition has declined in recent decades, both worldwide (2) and in the upper Midwest of North America (3). This could lead to declines in MeHg in fish, but information on trends for inland freshwater fish in temperate regions of North America is limited. Methylmercury contamination of fish is especially significant to Ojibwe tribes in northern Wisconsin because of increased harvest since 1985 after federal court decisions reaffirmed Ojibwe off-reservation treaty rights in the ceded territory (Figure S1, Supporting Information), and increased consumption during spring compared to the rest of the year (4). The Great Lakes Indian Fish and Wildlife Commission (GLIFWC) and the Wisconsin Department of Natural Resources (WDNR) have amassed substantial data on MeHg contamination in fish, particularly for walleye (Sander vitreus), which is popular with sportfishers and the species * Corresponding author phone: (608) 265-5639; fax: (608) 2622500; e-mail:
[email protected]. † Great Lakes Indian Fish and Wildlife Commission. ‡ University of California-Irvine. 4568
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 41, NO. 13, 2007
most often harvested by the Wisconsin Ojibwe tribes. Walleye were sampled primarily to provide consumption advice for individual lakes, although walleye were sampled repeatedly from some lakes specifically to evaluate time trends. Time trends in individual lakes are certainly of interest, but a natural question is whether there are trends across the entire region. Data collected from multiple years in a subset of lakes provides valuable information on trends within those lakes and across the region, but lakes with only 1 year of data also contain information on regional trends. Bayesian methods provide a useful and convenient framework to combine information, model dependence, and to identify large-scale characteristics in data such as this. A hierarchical Bayesian approach was used to model time trends and other characteristics of MeHg in walleye from individual lakes as samples from common distributions pertaining to the entire region. This analysis illustrates the ease with which missing data can be handled in a Bayesian framework, which allows the sex of the fish to be used as a predictor for all available samples. This is an important consideration since growth rates affect MeHg concentrations (5) and growth rates differ by sex for walleye and other species (6), but sex is often omitted as a predictor from analyses of MeHg by length or weight. Another advantage of Bayesian statistics is that direct probability statements can be made regarding the parameters since they are viewed as random rather than fixed. This analysis also demonstrates how data from multiple sources can be combined despite different collection methods, and how lake-specific covariates can be incorporated.
Materials and Methods Data Collection. A total of 5994 walleye from 420 lakes were collected and analyzed for mercury in the ceded territory of northern Wisconsin from 1982 to 2005. The data included 3060 walleye fillets collected by GLIFWC since 1989, 2721 fillets collected by WDNR since 1982, and 213 fillets collected by the Lac du Flambeau tribe from 1991 to 1994 and 2001. Most walleye were collected in the early spring shortly after lake ice had thawed. Most fillets collected by GLIFWC were analyzed with the skin off since tribal members usually consume them this way, while most fillets collected by WDNR were analyzed with the skin on. A previous study found MeHg concentrations to be about 10% lower in fillets with skin on compared to those with the skin removed (7). Total mercury was measured as a surrogate for MeHg because at least 95% of the mercury in top predator fish such as walleye is in the form of MeHg (8). The Lake Superior Research Institute at the University of Wisconsin-Superior analyzed 2987 of the walleye collected by GLIFWC based on U.S. Environmental Protection Agency (EPA) Method 245.6. The rest of the walleye collected by GLIFWC were analyzed at the University of Wisconsin-LaCrosse Mercury Laboratory using EPA Method 1631. All samples collected by WDNR were analyzed at the Wisconsin State Laboratory of Hygiene according to protocols found in ref 9. All samples collected by the Lac du Flambeau Ojibwe tribe were analyzed at Era Labs (Duluth, MN) based on EPA Method 245.6. Most of the 420 lakes in the analysis had 20 or fewer walleye sampled from them, and most had data collected in one or two different years (Table S1, Supporting Information). Walleye length was available for all samples. The walleye sampled included 2809 males, 1780 females, and 1405 samples for which sex was either not recorded or recorded as unknown since the sex could not be reliably determined because of sexual immaturity or lack of visible gametes at 10.1021/es0700294 CCC: $37.00
2007 American Chemical Society Published on Web 06/05/2007
the time of collection. Male walleye mature earlier and grow more slowly than females and have higher mortality rates (6, 10, 11). Covariates of pH and lake surface area were incorporated in the model since numerous studies have found them to be related to MeHg in fish (12-14) and they were readily available. For each lake, pH was calculated as an average of historical values from a statewide limnological study conducted from 1966 to 1980 (15) and of mean pH readings collected since 1986, so that both sources were weighted equally if both were available. One source was used if only one was available, so that pH data was available for 310 of the 420 lakes. Bayesian Analysis. Bayesian analysis is described in various texts (16, 17), but a brief overview follows. Traditional frequentist statistics treats data as a random sample drawn from a distribution characterized by fixed but unknown parameters, but Bayesian statistics treats model parameters as random and data as fixed. Prior probability distributions must be provided for the parameters, which are then updated based on the data to give posterior probability distributions. Since the integration necessary to derive exact posterior probability distributions for the parameters is usually prohibitively complex, approximate posterior distributions can be generated much more easily through simulation via Markov Chain Monte Carlo (MCMC) methods. After a “burnin” period for the chain to converge, each iteration of the chain represents a random draw of the parameters from the joint posterior probability distribution. The resulting draws can be summarized to give probability statements regarding the parameters or functions of the parameters that are of interest, essentially eliminating nuisance parameters that are not of interest by averaging over their possible values. For example, the probability that a parameter is greater than some specified value is given by the proportion of the draws of the parameter that are greater than that value. The most likely values of the parameter are characterized by the mean, median, and peak (or peaks) of the posterior distribution, which will all be very close to each other if the distribution is approximately normal. This provides direct probability statements regarding parameter values, in contrast to interpretations of frequentist p-values and confidence intervals that are less straightforward. The model was fit to log-transformed MeHg concentrations using freely available WinBUGS software (18, 19). Walleye length, collection year, pH, and the log of the surface area of the lake were standardized by subtracting their approximate means and dividing by their approximate standard deviations to facilitate convergence (Table S2, Supporting Information). A regression model was used, where Yij ∼ N(µij, σ2i ), where i indexes lakes 1 through 420, j indexes the fish sampled from each lake, Yij is the natural log of the MeHg concentration in fish j from lake i, N(µij, σ2i ) represents a normal distribution with mean µij and variance σ2i , and
µij ) Ri + R ˜ iGij + R*Iij + (βi + β˜ iGij + β*Iij)Xij + γiWij. (1) Observed data in eq 1 includes walleye sex indicator Gij (0 for female and 1 for male), sample form indicator Iij (0 for skin-on fillet and 1 for skin-off fillet), walleye length Xij, and collection year Wij. Parameters in eq 1 include intercept Ri for lake i, slope βi corresponding to walleye length for lake i, R ˜ i and β˜ i for differences in the intercept and length slope for males for lake i, R* and β* for differences in the intercept and length slope for skin-off fillets, and slope γi for the collection year for lake i. A hierarchical structure was used, so that lake-specific parameters were assumed to be drawn from their own parent population distributions, which were characterized by hyperparameters with their own posterior distributions. Pa-
rameters Ri and βi were assumed to be distributed according to regression models, where
Ri ∼ N(µRi ) RR + βRAAi + βRPPi, σ2R)
(2)
βi ∼ N(µβi ) Rβ + ββAAi + ββPPi, σ2β)
(3)
and
where Ai is the natural log of the surface area of lake i, and Pi is pH for lake i. It is prudent to include these covariates since they have been shown to affect MeHg levels in fish and are readily available. Other watershed, geographic, and lacustrine characteristics could easily be included from a modeling standpoint, but obtaining data on them for a reasonable subset of the 420 lakes in this analysis would be costly. It is not necessary to include additional covariates as long as the distributional assumptions in eqs 2 and 3 are appropriate, since factors influencing MeHg levels are not the focus of this analysis. Lake-specific parameters R ˜ i, β˜ i, γi, and σ2i were assumed to be drawn from simple normal and inverse-gamma distributions, so that R ˜ i ∼ N(µR˜ , σR2˜ ), β˜ i ∼ N(µβ˜ , σβ2˜ ), γi ∼ N(µγ, σ2γ), and 1/σ2i ∼ gamma(τR, τβ). In this way, the hyperparameter µγ for the mean of the year slopes γi for each of the 420 lakes serves as a measure of regional trends, and its posterior distribution describes the probability that there has been a regional trend in MeHg levels in walleye in northern Wisconsin over time. Lake-specific parameters R ˜i and β˜ i do not depend on lake surface area and pH since Ri and βi already incorporate these covariates, so R ˜ i and β˜ i express the differences in the intercepts and length slopes when the fish is male rather than female. Variances are parametrized in terms of their inverses in WinBUGS, so the distribution of precision terms 1/σ2i was specified rather than the distribution of σ2i . Parameters R* and β* pertaining to skin-off fillets are not lake-specific, but were estimated as single parameters with vague N(0, 10002) prior distributions. Walleye sex is an important predictor since the differing growth rates of the sexes should lead to differing MeHg concentrations for walleye of a given size (6). However, walleye sex was missing for 23% of the samples, so a model extension was needed to account for the missing data. Missing data procedures are described elsewhere (17, 20), but are briefly summarized here. It must first be decided whether the mechanism that results in missing data is ignorable or nonignorable. In this context, it would be nonignorable if walleye with sex missing or unknown were more likely to be female or male, so that the mechanism that leads to missing data should be modeled explicitly. There is no reason to think that this is the case here, so it was assumed that the missing data mechanism is ignorable and the data is missing at random. However, sex is more likely to be unknown for smaller, sexually immature walleye, so that the data is not observed at random even though it is missing at random. Bayesian analysis can account for such missing values by treating them as parameters whose distributions are influenced by the observed data, since knowing the size of the fish provides information regarding the probability that the walleye is male or female. Inferences regarding parameters of interest are made by averaging over distributions of missing data in the same way as other nuisance parameters. Missing data is easily handled in WinBUGS by treating it as random in the same way as the parameters, but with the initial values for the chain corresponding to the observed data coded as missing. In this way, walleye sex indicator Gij was treated as a random variable that is 1 with probability πij and 0 with probability 1 - πij, where πij was modeled VOL. 41, NO. 13, 2007 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
4569
using logistic regression with walleye length and year as predictors. A similar approach is described in a recent social science example (21). This gives Gij ∼ Bernoulli(πij) for each fish, where
logit(πij) ) Rπi + βπiXij + γπiWij. Lake-specific parameters Rπi, βπi, and γπi were all assumed to be drawn from normal distributions so that Rπi ∼ N(µRπ, σR2 π), βπi ∼ N(µβπ, σβ2π), and γπi ∼ N(µγπ, σγ2π). Missing pH values were assumed to be missing at random, and were modeled using a normal distribution, so that Pi ∼ N(µP, σ2P). Vague N(0, 10002) prior distributions were used for hyperparameters for the means of the population distributions (µR˜ , µγ, etc.), hyperparameters associated with Ri and βi (RR, βRA, etc.), and parameters R* and β* as previously mentioned. Vague gamma(0.001, 0.001) prior distributions were used for τR and τβ, and for hyperparameters for the precisions of the population distributions (1/σR2˜ , 1/σ2γ, etc.). The MCMC chains were thinned to reduce possible autocorrelation so that only every 50th iteration was kept, and 5000 draws were used after discarding the first 1000 as a burn-in period. Full posterior distributions were generated for all parameters and hyperparameters. To provide descriptions of regional trends and characteristics that are easy to understand, the lake-specific parameters in eq 1 were replaced by the hyperparameters for the means of their distributions, giving the function
µ ) RR + βRAA + βRPP + µR˜ G + R*I + (Rβ + ββAA + ββPP + µβ˜ G + β*I)X + µγW. (4) The predictor of interest was then varied over a reasonable range, while the other continuous predictors were held constant at their means. For example, the trend over time for a female walleye using a skin-on fillet was given by µ ) RR + µγW, where W varies from 1982 to 2005, and the other terms were eliminated because of the standardization of the continuous predictors and because the indicator variables were set to 0. One such line was generated for each of the 5000 MCMC iterations, and the median or percentiles of the resulting distributions at each level of the predictor of interest were used to characterize the relationship across the region. Trends and characteristics for individual lakes were derived similarly. An alternative is to replace each of the parameters in eq 4 with the means or medians of their respective posterior distributions. This was not done since it ignores the nature of the joint distribution of the parameters and their uncertainty, but it would have provided virtually identical results for this analysis. The regression coefficients and parameters do not necessarily have intuitive interpretations since log-transformed and standardized data were used. To facilitate interpretation of time trend parameter µγ, MeHg by year from eq 4 was exponentiated to the original scale in µg/g for each of the 5000 simulations. These were then regressed on year using simple least-squares regression, providing approximations of annual changes in MeHg in µg/g that are easier to interpret than log-scale values. Lake-specific time trend parameters, γi, were treated similarly. Exponentiation induces some curvature in the lines, but this was very slight so that the approximations are very close. The lines were parallel when different levels of the indicator variables for sex and sample form were used (or nearly parallel since posterior medians were found for each level of the predictor of interest), but they were not parallel after exponentiation. Thus, results are provided for various levels of the indicator variables since 4570
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 41, NO. 13, 2007
FIGURE 1. Regional MeHg concentrations by length for male and female walleye using skin-off and skin-on fillets. The lines were generated by varying length over its range and finding predicted MeHg for each sex and sample form combination using each of the 5000 MCMC iterations of the parameters, with the other predictors held constant at their means. The lines shown represent the medians of the 5000 MeHg values generated for each length increment. they give slightly different approximations of changes in MeHg over time.
Results and Discussion Sex Differences. Male walleye had higher MeHg concentrations for a given length than females across the region (Figure 1), which was expected since males grow more slowly and tend to be older than females of the same size. Concentrations of MeHg were about 0.06 µg/g higher in males than females at the sample average length of 445 mm. The difference was more pronounced for larger walleye. Male walleye rarely exceed 500 mm in northern Wisconsin, so most walleye larger than this are females. Males were consistently higher in MeHg than females of the same length across all lakes. All 420 lakespecific R ˜ i parameters describing the differences in the intercepts for males had probabilities of at least 0.908 of being greater than 0, and all but two had probabilities of at least 0.949 of being greater than 0. The distributions of the β˜ i parameters describing the differences in the length slopes for males indicated that males tend to increase in MeHg more quickly by length than females. All 420 β˜ i had probabilities of at least 0.701 of being greater than 0, and all but 12 had probabilities of at least 0.875 of being greater than 0. These results are consistent with other findings that MeHg accumulates at similar annual rates for the sexes, with growth rate differences leading to higher MeHg concentrations for the slower-growing sex for a given length or weight (6, 14). Age is a natural alternative predictor of MeHg in fish, but is often unavailable. Length is a readily available surrogate, but should be used on a sex-specific basis to improve regressions of MeHg on length that are frequently done without regard to sex. As a practical matter, accounting for sex is not needed for younger, smaller fish before growth rates begin to diverge. However, many fish that are consumed are sexually mature and larger. Sexual size dimorphism is common in fish consumed by humans, and should be considered when analyzing contaminant data to ensure that all relevant information is included and statistical assumptions are met. Sample Form Differences. Skin-on fillets were about 16% lower in MeHg than skin-off fillets, with a difference of about 0.07 µg/g for males and 0.06 µg/g for females at the reference length of 445 mm based on the posterior medians of R* and β* (Figure 1). This is expected since MeHg tends to accumulate
FIGURE 2. Posterior probability distributions for approximate regional annual change in MeHg based on hyperparameter µγ, shown for male and female walleye using skin-off and skin-on fillets. in muscle tissue rather than skin. The probability that R* is greater than 0 was found to be 1.0, since all simulations of R* were positive. The distribution of β* was centered near 0, indicating a high probability that the difference in MeHg due to the form of the sample does not change with walleye length. Subsequent results are shown for male skin-off fillets and female skin-on fillets at a minimum, since these produce the highest and lowest MeHg concentrations of the four combinations of sex and sample form. Sample form differences are confounded with differences between laboratories and testing procedures, since laboratories used by GLIFWC and the Lac du Flambeau tribe analyzed most fillets with the skin off, and the laboratory used by WDNR analyzed most fillets with the skin on. However, differences between laboratories and testing procedures appear to be minor since the sample form difference of 16% found here is relatively close to the 10% found in the previous study (7). Differences between laboratories and testing procedures could easily be modeled through simple extensions of the indicator variables used here, given the necessary data, to accommodate analyses of data from multiple agencies using different collection and testing methods. Time Trends. A slight decreasing regional trend in MeHg was found. The most likely values of the approximate annual changes were -0.0028 µg/g per year for skin-off fillets from male walleye, and -0.0020 µg/g per year for skin-on fillets from female walleye based on the distribution of hyperparameter µγ (Figure 2). This translates to a regional decrease in MeHg of 0.065 µg/g from 1982 to 2005 for skin-off fillets from male walleye, and a decrease of 0.047 µg/g for skin-on fillets from female walleye (Figure 3). This was a 13.0% decline for each from 1982 to 2005, or an annual decline of 0.60%. There was a 0.9992 probability of a decreasing regional trend in MeHg in walleye since all but 4 of the 5000 simulations of µγ were less than 0. Most individual lakes showed decreasing trends in walleye MeHg over time. Of the 420 lakes, 325 of them (77%) had probabilities of 0.5 or greater of decreasing walleye MeHg over time, having posterior medians of γi less than 0 (Figure 4). However, there was considerably more uncertainty regarding trends on individual lakes than for the regional trend described by µγ (Table S3, Supporting Information). Posterior distributions were generated for time change parameters γi even for lakes with only 1 year of data because of the hierarchical structure used, for which the posterior distributions of γi tended to be wide and centered near the more likely values of µγ. Of the 133 lakes with three or more years of data, where the distributions for γi tended to provide
FIGURE 3. Regional trends in MeHg in walleye fillets in northern Wisconsin. The lines were generated by varying year over its range and finding predicted MeHg for each sex and sample form combination using each of the 5000 MCMC iterations of the parameters, with the other predictors held constant at their means. The solid lines represent the medians of the 5000 MeHg values generated for each year increment, and the dashed lines for male skin-off and female skin-on fillets represent the 5th and 95th percentiles of the 5000 MeHg values for each year increment.
FIGURE 4. Approximate annual change in MeHg for each of the 420 lakes, from the medians of the slopes from linear regressions fit to predicted MeHg by year exponentiated to the original µg/g scale from each of the 5000 MCMC iterations. A reference line is shown at 0. Bars show annual change for skin-off fillets from male walleye and skin-on fillets for female walleye for each lake, subdivided according to whether at least 3 years of data were available. more certainty by being narrower, 95 of them (71%) had probabilities of 0.5 or greater of decreasing walleye MeHg over time. Existing MeHg trend data for inland freshwater fish in North America is limited and inconclusive. Methylmercury declined in three species in Lake Ontario from 1977 to 1988 (22, 23). Methylmercury declined in sport fish in Lake St. Clair from 1970 to 1994, but not in Lakes Huron and Ontario from 1981 to 1994 (24). Methylmercury does not appear to be declining in freshwater fish in the Canadian Subarctic and Arctic (25), although climate change could result in increased MeHg concentrations in this region even with declining atmospheric mercury deposition (26). The decrease in MeHg in walleye fillets in northern Wisconsin since 1982, both regionally and in most individual lakes, is consistent with decreased atmospheric deposition worldwide and in the upper Midwest. This is concurrent with decreased mercury emissions in North America (27) and the Lake Superior basin (28). Atmospheric mercury was estimated to have declined about 3% annually worldwide VOL. 41, NO. 13, 2007 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
4571
FIGURE 5. Posterior medians of MeHg in 445 mm walleye in each of 420 lakes for skin-off fillets from males and skin-on fillets from females, with year held constant at 1995.
from 1990 to 1996 (2). From their peaks in the 1950s through 1970s, mercury measured from sediment cores from four northeastern Minnesota lakes declined an average of 2% annually to the 1980s (3). These declines are also consistent with decreasing aqueous mercury and MeHg observed in Little Rock Lake in Wisconsin (29). In the reference basin of Little Rock Lake (rather than the experimentally acidified basin), a 5% annual decline in aqueous MeHg was observed from 1988 to 1999, and a 5% annual decline was observed in MeHg in yellow perch from 1994 to 2000 (30). No time trends are available for walleye in Little Rock Lake since walleye are not present. The regional annual decline of 0.60% for walleye MeHg in this analysis is small compared to these other measures. Approximate annual percentage declines for MeHg from 1982 to 2005 in skin-off fillets from male walleye were found for each lake based on time trend parameters γi. The median of the annual percentage declines for all 420 lakes was 0.62%, and the greatest annual percentage decline was 4.57%, with only 10 lakes having annual percentage declines greater than 3%. Walleye live longer than perch and typically have a life span of about 7 years in northern Wisconsin with a few reaching 10-12 years old (11), which may partially explain the slower decline of MeHg in walleye. The lakes in this analysis have a mean size of 280 ha and a median size of 132 ha, but Little Rock Lake is relatively small at 18 ha, which could also be related to its more rapid declines. Future investigations could address additional chemical, physical, and environmental reasons for the apparently slow declines of MeHg in walleye from northern Wisconsin. Between-Lake Variability and Covariates. Methylmercury concentrations varied widely from lake to lake, with posterior medians for MeHg for skin-off fillets from 445 mm male walleye ranging from 0.14 to 1.73 µg/g (Figure 5). These posterior medians were greater than 0.47 µg/g for 213 of the 420 lakes. The EPA recommends eating one or fewer meals per month of fish with MeHg concentrations greater than 0.47 µg/g (31), so MeHg levels for many lakes could affect tribal members whose consumption increases during the spring (4). There was a strong inverse relationship regionally between MeHg and pH (Figure S2, Supporting Information) consistent with previous studies (12, 14), with an increase in MeHg of about 0.17 µg/g per unit decrease of pH. The entire distribution of hyperparameter βRP associated with intercepts Ri lay well away from 0. However, the distribution of hyperparameter ββP associated with length slopes βi overlapped 0 so there was a 0.089 probability that ββP was less than 0. This indicates that there is more certainty that pH 4572
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 41, NO. 13, 2007
affects overall MeHg concentrations than that pH affects the rate at which MeHg accumulates in walleye. Methylmercury concentrations in walleye decreased rapidly as lake surface area increased for smaller lakes less than 100 ha or so, but decreased gradually for larger lakes (Figure S3, Supporting Information). The entire distribution of hyperparameter βRA associated with intercepts Ri was less than 0, but the distribution of hyperparameter ββA associated with length slopes βi was centered almost exactly at 0. Most of the variation between lakes occurred through the intercept terms Ri, whose posterior medians had a coefficient of variation (CV) of 42%, rather than the length slope terms βi, whose posterior medians had a CV of 13%. The distributions of the Ri and βi parameters along with their respective hyperparameters indicated that MeHg accumulated in walleye by length at generally similar rates across all lakes after accounting for sex even though overall MeHg concentrations varied widely between lakes. Although the reasons for the variation between lakes were not the focus of this analysis, the results were consistent with other studies that found that between-lake differences in MeHg levels in fish are driven largely by lake-specific watershed and lacustrine characteristics affecting the conversion of atmospherically deposited mercury to MeHg (12, 32). Within-Lake Variability. Within-lake variability characterized by σ2i was high (Figure S4, Supporting Information). Even toward the lower end of within-lake variability, MeHg varied by 0.38 µg/g based on the upper and lower bounds of an interval containing 95% of the samples for 445 mm male skin-off fillets at the fifth percentile of the distribution of the σ2i . The distribution of the σ2i parameters was skewed to the right, so that within-lake variability can potentially be quite high. This implies that when issuing consumption advice for sensitive populations, it would be prudent to consider that some fish may have high concentrations of MeHg even if average levels are low for a particular lake. When within-lake variability is high, it is difficult to detect trends for individual lakes unless sample sizes are large, samples are spread over a wide range of years, or the trends are of large magnitude. Even when trends can be detected for individual lakes, between-lake variability makes it difficult to draw conclusions about regional trends from a limited number of lakes. Regional trends across time and space are more likely to be detected despite high within- and betweenlake variability when all data across the region are incorporated. The Bayesian analysis presented here provides a useful template for accommodating all available data from multiple agencies despite missing data and potential differences in collection and testing methods, which can facilitate analyses across broad spatial scales.
Acknowledgments Thanks to Neil Kmiecik for overseeing GLIFWC’s mercury testing and reviewing the manuscript. Thanks to WDNR and the Lac du Flambeau Water Resources Department for providing mercury data, and to Jen Filbert of WDNR for providing pH data. Thanks to Matt Hudson, Adam DeWeese, and Kam Tsui for their helpful comments, and to Esteban Chiriboga for preparing Figure S1. Testing of fish collected by GLIFWC was supported by various grants from the Agency for Toxic Substances and Disease Registry and the EPA.
Note Added After ASAP Publication Due to a production error, there were minor text changes in the version published ASAP June 5, 2007; the corrected version was published ASAP June 5, 2007.
Supporting Information Available Tables showing the number of fillets analyzed and years sampled for each lake, summaries of predictor variables,
distributions of selected parameters, a map showing the Wisconsin ceded territory, graphs of MeHg by pH and lake surface area, within-lake variability, model checking details, and WinBUGS code. This material is available free of charge via the Internet at http://pubs.acs.org.
Literature Cited (1) Pirrone, N.; Mahaffey, K. R. Where we stand on mercury pollution and its health effects on regional and global scales. In Dynamics of Mercury Pollution on Regional and Global Scales; Pirrone, N., Mahaffey, K. R., Eds.; Springer: New York, 2005; pp 1-21. (2) Slemr, F.; Brunke, E. G.; Ebinghaus, R.; Temme, C.; Munthe, J.; Wa¨ngberg, I.; Schroeder, W.; Steffen, A.; Berg, T. Worldwide trend of atmospheric mercury since 1977. Geophys. Res. Lett. 2003, 30 (23). (3) Engstrom, D. R.; Swain, E. B. Recent declines in atmospheric mercury deposition in the upper midwest. Environ. Sci. Technol. 1997, 31, 960-967. (4) Peterson, D. E.; Kanarek, M. S.; Kuykendall, M. A.; Diedrich, J. M.; Anderson, H. A.; Remington, P. L.; Sheffy, T. B. Fish consumption patterns and blood mercury levels in Wisconsin Chippewa Indians. Arch. Environ. Health 1994, 49, 53-58. (5) Simoneau, M.; Lucotte, M.; Garceau, S.; Laliberte´, D. Fish growth rates modulate mercury concentrations in walleye (Sander vitreus) from eastern Canadian lakes. Environ. Res. 2005, 98, 73-82. (6) Henderson, B. A.; Collins, N; Morgan, G. E.; Vaillancourt, A. Sexual size dimorphism of walleye (Stizostedion vitreum vitreum). Can. J. Fish. Aquat. Sci. 2003, 60, 1345-1352. (7) Dellinger, J.; Kmiecik, N.; Gerstenberger, S.; Ngu, H. Mercury contamination of fish in the Ojibwa diet: I. Walleye fillets and skin-on versus skin-off sampling. Water, Air, Soil Pollut. 1995, 80, 69-76. (8) Bloom, N. S. On the chemical form of mercury in edible fish and marine invertebrate tissue. Can. J. Fish. Aquat. Sci. 1992, 49, 1010-1017. (9) Sullivan, J. R.; Delfino, J. J. The determination of mercury in fish. J. Environ. Sci. Health 1982, A17, 265-275. (10) Henderson, B. A.; Nepszy, S. J. Reproductive tactics of walleye (Stizostedion vitreum vitreum) in Lake Erie. Can. J. Fish. Aquat. Sci. 1994, 51, 986-997. (11) Becker, G. C. Fishes of Wisconsin; University of Wisconsin Press: Madison, WI, 1983. (12) Wiener, J. G.; Martini, R. E.; Sheffy, T. B.; Glass, G. E. Factors influencing mercury concentrations in walleyes in northern Wisconsin lakes. Trans. Am. Fish. Soc. 1990, 119, 862870. (13) Bodaly, R. A.; Rudd, J. W. M.; Fudge, R. J. P.; Kelly, C. A. Mercury concentrations in fish related to size of remote Canadian shield lakes. Can. J. Fish. Aquat. Sci. 1993, 50, 980-987. (14) Lange, T. R.; Royals, H. E.; Connor, L. L. Influence of water chemistry on mercury concentration in largemouth bass from Florida lakes. Trans. Am. Fish. Soc. 1993, 122, 74-84. (15) Lillie, R. A.; Mason, J. W. Limnological Characteristics of Wisconsin Lakes; Technical Bulletin No. 138; Wisconsin Department of Natural Resources: Madison, WI, 1983. (16) Carlin, B. P.; Louis, T. A. Bayes and Empirical Bayes Methods for Data Analysis; Chapman and Hall/CRC Press: Boca Raton, FL, 1996. (17) Gelman, A.; Carlin, J. B.; Stern, H. S.; Rubin, D. B. Bayesian Data Analysis, 2nd Ed.; Chapman and Hall/CRC Press: Boca Raton, FL, 2004.
(18) Gilks, W. R.; Richardson, S.; Spiegelhalter, D. J. Markov Chain Monte Carlo in Practice; Chapman and Hall/CRC Press: Boca Raton, FL, 1996. (19) Spiegelhalter, D. J.; Thomas, A.; Best, N. G.; Lunn, D. WinBUGS Version 1.4 User Manual; Medical Research Council Biostatistics Unit, Institute of Public Health: Cambridge, UK, 2003. (20) Little, R. J. A.; Rubin, D. B. Statistical Analysis with Missing Data, 2nd Edition; Wiley: New York, 1992. (21) Pardoe, I.; Weidner, R. R. Sentencing convicted felons in the United States: a Bayesian analysis using multilevel covariates. J. Stat. Plan. Inference 2006, 136, 1433-1455. (22) Borgman, U.; Whittle, D. M. Contaminant concentration trends in Lake Ontario lake trout (Salvelinus namaycush): 1977 to 1988. J. Great Lakes Res. 1991, 17, 368-381. (23) Borgman, U.; Whittle, D. M. DDE, PCB, and mercury concentration trends in Lake Ontario rainbow smelt (Osmerus mordax) and slimy sculpin (Cottus cognatus): 1977 to 1988. J. Great Lake Res. 1992, 18, 298-308. (24) Scheider, W. A.; Cox, C.; Hayton, A.; Hitchin, G.; Vaillancourt, A. Current status and temporal trends in concentrations of persistent toxic substances in sport fish and juvenile forage fish in the Canadian waters of the Great Lakes. Environ. Monit. Assess. 1998, 53, 57-76. (25) Evans, M. S.; Muir, D.; Lockhart, W. L.; Stern, G.; Ryan, M.; Roach, P. Persistent organic pollutants and metals in the freshwater biota of the Canadian Subarctic and Arctic: An overview. Sci. Total Environ. 2005, 351-352, 94-147. (26) Macdonald, R. W.; Harner, T.; Fyfe, J. Recent climate change in the Arctic and its impact on contaminant pathways and interpretation of temporal trend data. Sci. Total Environ. 2005, 342, 5-86. (27) Pirrone, N.; Allegrini, I.; Keeler, G. J.; Nriagu, J. O.; Rossman, R.; Robbins, J. A. Historical atmospheric mercury emissions and depositions in North America compared to mercury accumulations in sedimentary records. Atmos. Environ. 1998, 32, 929940. (28) Lake Superior Binational Program. Lake Superior Lakewide Management Plan: 1990-2005. Critical Chemical Reduction Milestones; Superior Work GroupsChemical Committee: Toronto and Chicago, 2006. (29) Watras, C. J.; Morrison, K. A.; Hudson, R. J. M.; Frost, T. M.; Kratz, T. K. Decreasing mercury in northern Wisconsin: temporal patterns in bulk precipitation and a precipitationdominated lake. Environ. Sci. Technol. 2000, 34, 4051-4057. (30) Hrabik, T. R.; Watras, C. J. Recent declines in mercury concentration in a freshwater fishery: isolating the effects of de-acidification and decreased atmospheric mercury deposition in Little Rock Lake. Sci. Total Environ. 2002, 297, 229-237. (31) USEPA. Guidance for Assessing Chemical Contaminant Data for Use in Fish Advisories, Volume 2: Risk Assessment and Fish Consumption Limits, 3rd edition; EPA-823/B-00-008; U.S. Environmental Protection Agency: Washington, DC, 2000. (32) Wiener, J. G.; Knights, B. C.; Sandheinrich, M. B.; Jeremiason, J. D.; Brigham, M. E.; Engstrom, D. R.; Woodruff, L. G.; Cannon, W. F.; Balogh, S. J. Mercury in soils, lakes, and fish in Voyageurs National Park (Minnesota): importance of atmospheric deposition and ecosystem factors. Environ. Sci. Technol. 2006, 40, 6261-6268.
Received for review January 5, 2007. Revised manuscript received March 31, 2007. Accepted April 16, 2007. ES0700294
VOL. 41, NO. 13, 2007 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
4573
