**Mercury in Hair of Mammoth and Other Prehistorical Mammals as a Proxy of Hg Level in the Environment Associated with Climate Changes**

### **Stella Eyrikh 1,\*, Gennady Boeskorov 2, Tatyana Serykh 1, Marina Shchelchkova <sup>3</sup> and Tatyana Papina <sup>1</sup>**


Received: 5 October 2020; Accepted: 30 November 2020; Published: 3 December 2020

**Abstract:** The paper presents the first results of Hg determination in the hair of prehistorical animals (woolly mammoth, steppe bison, and woolly rhino). Hair of prehistorical mammals can be used as an archive that preserves changes of environmental pollution at the paleoscale. The aim of our study was to assess the levels of Hg exposure of ancient animals and to understand whether Hg concentration in hair could be used as a proxy indicating changes of mercury levels in the environment following global climate changes. We assessed changes of Hg exposure recorded in hairs of seven specimens of mammoth fauna mammals that inhabited the Yakutia region in the period from 45 to 10 ka yr BP. Hg concentrations in hair varied from 0.017 to 0.177 μg/g; the lowest Hg concentration were determined in older specimens (45–33 kyr yr BP). The two highest concentrations belonged sample from the Last Glacial Maximum and the Karginian interstadial (57–24 kyr BP) periods. Our hypothesis is the increase of Hg concentrations in hair reflecting environmental Hg level might be forced by high dust load in cold periods and thawing permafrost in warm climatic periods. Long-term variations of Hg level recovered from Ice Age animals' hair correlate with Hg profiles of concentration and deposition reconstructed from the Antarctica ice core.

**Keywords:** mercury; mammoth fauna mammals; hair; environmental changes; paleoclimate; Pleistocene; Yakutia

#### **1. Introduction**

The content of macro- and microelements in human and animal hairs is a good indicator of their accumulation in the body as a result of environmental exposure, including intake with food and water [1,2]. Hair records the levels of toxic (lead, cadmium, arsenic, etc.) and vital elements (zinc, selenium, iron, etc.), reflecting the elemental status of the whole organism. Hair analysis is used for evaluation of health state, metabolic disorders, mineral maintenance of human and animals, and also the ecological state of the territory where they live [3–5]. The level of toxic metals in the environment indicates a potential risk for the ecosystem and for human and animal health because of bioaccumulation of some metals (particularly mercury) in the body [6,7]. The World Health Organization (WHO) recommends using hair as major biological material for testing the pollution of the human body by heavy metals, since sampling, storage, and analysis of hair samples are easier than they are for other biological materials [8]. The International Atomic Energy Agency (IAEA) uses

hair for the monitoring of global changes in element levels in the environment worldwide [9–11]. Ancient hairs are also "keepers of history" [12–14] which help to assess the level of environmental pollution by the degree of pollutants impact on the body during previous epochs. Hairs of prehistoric animals can be the key to understanding the relationship of environmental changes with climate. The aim of our study was to understand whether Hg concentration in the hair of prehistoric mammals could be used as a proxy indicating that changes of mercury level in the environment reflect climate changes during the Late Pleistocene–Early Holocene. Permafrost and substantial precipitation are the deciding factors in preservation of mammoth soft tissues and hair over tens of thousands of years. Hairs of the woolly mammoth are studied very actively nowadays for decoding and sequencing DNA [15–17]; revealing biologic rhythms [18]; determining the type of nutrition from the balance of stable isotopes of nitrogen, carbon, and phosphorus [19,20]; and their response to short-term (seasonal) environmental changes [21]. It is therefore surprising that studies of trace elements in the hair of mammoth fauna mammals have not been done until now. There are some studies of trace element in museum samples of animal hair [22] and bird feathers [23,24], and seal hairs from a lake sediment core spanning the past 2000 years [25]. The results of these studies are useful for assessing environmental changes and anthropogenic impacts on the environment. However, most of the studies cover the span of the last thousand years, whereas the analysis of mammoth fauna mammals' hair provides a unique opportunity to evaluate environmental changes that were happening tens of thousands of years ago, during different climatic stages. It should be noted that there is a potential problem with reliability of analytical data associated with possible contamination or loss of elements during storage and analysis of samples [26]. Therefore, the development of methodological details of sample preparation and analysis of the prehistoric animals' hair require a special attention.

#### **2. Materials and Methods**

#### *2.1. Study Sites*

All the studied fossil mammals were discovered on the territory of Yakutia (Eastern Siberia, Russia) (Figure 1).

**Figure 1.** Study area map with location of fossil mammals discovered. Study sites were numbered from south to north.

For thousands of years, special climatic and geological conditions prevailed in the territory of Yakutia, the harshest climate in Eurasia: long cold season; low-temperatures anaerobic swamps and floodplains of rivers and lakes; gradual accumulation of precipitation on the floodplains; and permafrost growth [27,28]. These conditions allowed for the preservation of unique ancient specimens of mammals in permafrost deposits.

We studied seven specimens of Ice Age mammals; their detailed description (location and time of the discovery, estimation of the geological age by radiocarbon dating, sex, and physiological age of the animal) is presented in Table 1. The calibration of radiocarbon data was carried out by using the IntCal20 calibration curve, the version of program OxCal 4.4 (https://c14.arch.ox.ac.uk/oxcal/OxCal.html).


**Table 1.** Description of the studied Ice Age mammals.

\* Climatic stages are explained in Section 2.2. LG, Late Glacial; LGM, Last Glacial Maximum; KARG, Karginian interstadial; HAS, Hasselo stadial; DEN, Denekamp interstadial.

#### *2.2. Features of Discovered Species and Individuals*

The woolly mammoth (*Mammuthus primigenius*) is an extinct species that lived during the Pleistocene, until its extinction in the early Holocene epoch. Woolly mammoths lived in open grassland biomes, the mammoth steppe. High-productivity grasses, herbs, and shrubs dominated there. Stomach contents clearly show that the diet of the woolly mammoth was mainly grasses and sedges [36], although tree bark and twigs also constituted a small part of their winter diet [37,38]. Yukagir mammoth (Figure 2a) is an old male who lived during the Last Glacial Maximum (LGM), the maximum Sartanian glaciations in Siberia. He died by falling into a hole. Malolyakhovsky mammoth is an old female discovered on the Maly Lyakhovsky Island (New Siberian Islands in the Laptev Sea). Radiocarbon dating of bones and hairs demonstrated good agreement, attributing the lifetime of the mammoth to the Karginian interstadial (KARG). It was a relatively warm climatic phase in Late Pleistocene in Siberia (57–24 kyr BP), which overlaps with the Denekamp interstadial warming in Europe. Berelyokh mammoth was found in one of the biggest fossil sites, the Berelyokh Mammoth Cemetery located in the basin of the Indigirka River. The entire back leg, 175 cm long, was discovered (Figure 2c), with the longest hairs reaching 120 cm. Probably it belongs to the adult female which lived in the Late Glacial (LG). Oymyakon mammoth is a baby female; she died by falling into a permafrost crack. Only the upper part of her body is preserved well. This is the oldest of the samples

from the beginning of the Karginian interstadial in Siberia (Hasselo stadial in European classification). Stocky limbs and thick wool of the woolly rhino (*Coelodonta antiquitatis*) are well suited to the cold and arid steppe–tundra environment prevalent during the Pleistocene glaciations. "Churapcha rhino" (Figure 2d) is an adult female who lived during the Last Glacial Maximum (LGM); she died by falling into a coastal swamp shortly before it froze.

**Figure 2.** Photographs of the studied fossils of mammoth fauna mammals: (**a**) Yukagir mammoth head, (**b**) Yukagir bison, (**c**) posterior leg of Berelyokh mammoth, and (**d**) Churapcha rhino skeleton.

The steppe bison (*Bison priscus*), the mammoth, and the woolly rhino were the last largest herbivores that survived in Eurasia by the turn of the Pleistocene and Holocene. Steppe bison ate sedges, cereals, and plants from the forbs group. Fossil bison mummies are extremely rare. Only four well-preserved corpses are known, with two of them having been discovered in Yakutia: Yukagir and Malykhchin bison. The Yukagir bison (Figure 2b) is a young male, a complete frozen mummy that was the youngest of the studied fossil belonging to early Holocene, confirming that the bison survived an abrupt climate change at the Pleistocene–Holocene Boundary. The Malykhchin bison is a young female found on the right bank of the Indigirka River, where she died by getting stuck in coastal mud in the summer. She ate forest herbs, branches, and foliage. The presence of mosses in her stomach indicates the existence of wetland biotopes during the life of bison in the warm period of the Karginian interstadial.

Isotopic biogeochemistry helps to reveal the ecological structure of the mammoth steppe fauna. Isotopic differences reflect different dietary choices by herbivores. Woolly rhinoceros and bison grazed fresh grass, and mammoths consumed dry grass. Despite some differences in nitrogen and carbon isotopes, woolly mammoths and woolly rhinoceros are considered globally similar in diet (grass) and physiology (monogastric) [19]. Thus, it can be assumed that the differences in the concentration of mercury between the hairs of different animals reflect the changes in concentration of Hg in the environment.

Hair samples of woolly mammoth, woolly rhino, and steppe bison were obtained from the Geological Museum of Diamond and Precious Metals Geology Institute (Yakutsk). Hair records the cumulative exposure to mercury in the short- to medium-term, depending on the length of the hair sample. Whereas human hair growth rate is about 1 cm per month, and the concentration of metal in hair can show the level of mercury exposure that has occurred over many years, the animal's hair is replaced every 1–1.5 year (it holds for most of animals, both ancient and modern). Therefore, the full length of hair represents a continuous record of the elements intake over this period. Mammoth hair grows approximately 31 cm/year; the longest hair ever found covers 39 month of a mammoth's life [18,21]. Thus, the hair of mammoths and other mammoth fauna mammals reflects the environmental situation in the last years of their lives.

#### *2.3. Sample Preparation*

The determination of mercury in the hair of prehistoric animals and museum exhibits is associated with a number of difficulties: obtaining a representative sample, choosing an appropriate method given the small amount of sample, and reliable analytical determination of Hg concentration in it. Loss of volatile Hg and sample contamination are possible during long-term storage and transportation of the sample. Samples of studied fossil animals' hair were stored in museum in glass cases or wooden boxes in conditions excluding their mercury contamination during storage and therefore, they are suitable for analysis.

The amount of prehistoric animals' hair is very limited, and a single procedure of sample preparation must be developed not only for Hg, but also for a wide range of other trace elements. Here we aimed to determine the total concentrations of mercury in the hair of prehistoric animals both endogenous and exogenous in origin reflecting the intake from food and water, as well as from the air.

Methylmercury easily incorporated into hairs as it grows and its concentration in the hair is proportional to the blood concentration. The high affinity of hair for metals is mainly due to the presence of cysteine or sulfhydryl (SH) groups [39]. Elemental mercury may also bind to the hydrophobic core of the melanin polymer in the hair structure [40]. The IAEA recommends hair washing procedure using acetone and deionized water [10]; it is not suitable for Hg because the fat and keratin structures of the hair are destroyed by acetone, which leads to the loss of endogenous mercury. Washing the hair with HCl solution can leach methyl mercury from hair samples [41]. We used a chemically inert detergent ("SYNERGETIC Baby", fragrance and color free) which removes only surface grease and dust from hair samples without disturbing their structure. All reagents were tested for Hg content and purified if necessary. Nitric and hydrochloric acids were purified using a Savillex DST-1000 distillation system (Savillex, Eden Prairie, MN, USA). Ultrapure water (MQ-water) was obtained using a Simplicity UV water purification system (Millipore SAS, Molsheim, France). All stages of sample preparation were carried out in a "clean room" equipped with outdoor air handlers that use progressively finer filters including high-efficiency particulate air (HEPA) filter and charcoal mercury filter, which remove particulate and elemental mercury from the incoming air.

The washing procedure for hair comprises the following steps:


Microwave system MARS-5 (Thermo Fisher Scientific, Waltham, MA, USA) was used for digestion of hair samples by the program previously optimized for Hg analysis in biological objects (Table 2) [42]. We tested acid and acid–peroxide digestion and demonstrated applicability of acid–peroxide digestion (2 mL HNO3 + 1 mL H2O2) for Hg analysis. The latter method was used for the analysis of samples.


**Table 2.** Optimized parameters of microwave digestion of biotic samples.

After microwave digestion the samples were cooled to 25 ◦C, the pressure was brought to <50 psi, and sample volume was adjusted to 10–12 mL with MQ water. The procedure of sample preparation (washing, cutting, and digestion) was developed and tested by using hair samples of modern yak living in the Barnaul Zoo and Certified Reference Material of human hair (CRM, Hair NSC DC 73347, China). The developed sample-preparation procedure is suitable for both for Hg and multi-element analysis of prehistorical animals' hair and blood.

#### *2.4. Hg Analysis*

Content of mercury in hair and blood samples was determined by Mercur Duo Plus Analyzer (Analytik Jena, Jena, Germany), combining atomic fluorescence with the cold vapor method and amalgamation on gold collector. Analytical characteristics of the method are presented in Table 3. The accuracy was confirmed by using Certified Reference Material of human hair (CRM, Hair DC 73347, China). Optimization of the instrumental parameters [43], using ultrapure reagents and clean conditions, allowed us to achieve a method detection limit of up to 0.4 ng/L for liquid samples and 0.003 μg/g for hair samples (0.03 g, dry weight). Split sampling and analyzing the same samples at different times and by different operators were used for assessing precision, recovery, and reproducibility. Good spike recovery values were demonstrated for samples of yak hair and CRM. The confidence interval for low-concentration samples did not exceed 17%.

**Table 3.** Analytical characteristics of method.


<sup>1</sup> Certified value is for human hair, Certified Reference Material (CRM) NSC DC 73347, China.

#### **3. Results and Discussion**

#### *3.1. Hg Concentration in Fossil Animals' Hair*

Mercury concentrations determined in the prehistoric animals' hair varied from 0.017 to 0.177 μg/g; average concentrations and ranges of Hg content in different types of prehistoric and modern animals are presented in Table 4, together with reference values and intervals. Mercury coming directly from water, air, and food tends to accumulate in both plants and animals, being toxic to most life forms. WHO guidance established 2 μg/g for total Hg in human hair as the reference level for risk evaluation [44]. The US Environmental Protection Agency (EPA) sets the reference dose for human hair and wildlife toxicity at 1 μg/g [45]. All of the prehistoric animals' hairs have Hg concentrations significantly below these levels. Moreover, they do not exceed the background level of mercury in hair of non-seafood consumers (0.5 μg/g). As far as we are aware, there is no background assessment of Hg level in herbivore prehistorical animals. We can compare our results with Hg concentrations in the hair of modern herbivore animals and their reference intervals for Hg (discussed in Section 3.2).


**Table 4.** Hg concentrations in hair of modern and prehistoric animals and reference values.

\* μg/L; ND—not detected.

#### 3.1.1. Woolly Mammoth Hair and Hemolyzed Blood

Hairs of four fossil mammoths were studied; the Hg concentrations are shown in the Figure 3a. Higher concentrations were in hair samples of Berelyekh and Yukagir mammoths from the Last Glacial and Preboreal warming, respectively, whereas the lowest Hg concentrations were in older specimens, Oymyakonsky and Malolyakhovsky mammoths (45–33 ka yr BP).

**Figure 3.** Mercury concentration in hair of different fossil specimens: (**a**) mammoth; (**b**) bison and rhino.

There was rare opportunity to sample enough freshly thawed hair material from different parts of the body of the same animal (Malolyakhovsky mammoth). Samples of hair taken from the neck and back leg demonstrated slightly different concentrations: 0.017 and 0.024 μg/g, respectively. The difference between Hg concentrations sampled from different body areas for modern animals varied insignificantly for yak from 0.004 to 0.007 μg/g (this study) and for beef cattle from 0.062 to 0.070 μg/g [46]. The comparison of guard and down yak hair demonstrated a negligible difference in Hg concentrations: 0.0054 and 0.0058 μg/g, respectively. Cattle hair also yielded similar Hg concentrations: 0.010 and 0.008 μg/g for guard and down hairs. It was shown that the elemental composition and concentrations of most elements of beef cattle hair on different body areas does not differ notably [46]. Of course, in the ideal case, hair sampled from the same body areas should be used for the comparison of the animals. However, unfortunately, there are very few opportunities to get hair samples from the same part of body for different prehistoric animals. Therefore, comparison of prehistoric animals' hair sampled from different body surface areas is evidently acceptable due to minor differences in the concentrations as confirmed by our results (no more than 30%).

There was rare chance to determine Hg concentration in hemolyzed blood of Malolyakhovsky mammoth. Mercury measurement in the whole blood provides information about recent exposure (~1–2 months) to both organic and inorganic mercury through ingestion of food and drinking water and inhalation of elemental mercury vapor in ambient air. The level of mercury in blood indicates recent exposure, but it does not reflect historical exposure or variations in exposure. Here we determined Hg concentration in mammoth hemolyzed blood 0.69 ± 0.07 μg/L. Hg levels in mammoth blood are below background levels of mercury in the blood of people who do not consume fish

(<2 μg/L) [45]. It is comparable with relatively low levels of total mercury in the blood of modern animals: The median value of Hg in blood of Galician cows from NW Spain was <0.438 μg/L, with the range of <0.438 to 15.4 μg/L [49]. This is similar to Hg concentrations in blood of surveyed dogs from Alaska (0.16–12.38 ng/g) [50].

Hair to blood Hg ratio in Malolyakhovsky mammoth was 36, which is closer to the hair to blood ratio in fish-fed dogs (59 ± 7.6), harbor seals (22–40), and polar bears (100) [51], than to the WHO value for humans (250) used for risk assessments to predict blood Hg from hair concentrations [8]. The differences in ratios may be due to the differences in relative surface area and hair density for different animal species and humans. Unfortunately, we have a single blood sample for mammoth, so we can only assume that this ratio level is characteristic for mammoths and other Ice Age animals.

#### 3.1.2. Steppe Bison Hair

Hg concentrations in the hair of two steppe bison differ by about five times: 0.035 ± 0.001 and 0.177 ± 0.026 μg/g (Figure 3b). The concentration of Hg in Yukagir bison hair is at a comparable level to the Berelekh mammoth hair, whereas Hg concentration in hair of Mylakhchin bison is the highest among all studied Ice Age animals. We hypothesize that there was an increase in mercury concentrations in the environment during this period. Mylakhchin bison lived in Karginian interstadial, in conditions of climate warming. The latter might be responsible for enhanced release of Hg due to thawing permafrost. Reconstruction based on palynological data revealed that, during the Karginian interstadial, there were stages with warmer and milder-than-today climate conditions, and the amplitude of climate fluctuations was different for different regions of Siberia [52]. Modern bison are very similar to the prehistoric ones in terms of nutrition, wool structure, etc. It has been observed that bison can find food under deep snow layers (>50 cm) [53]. Unfortunately, we could not find mercury concentrations in the hair of modern bison, although very low levels of hepatic Hg in the liver of captive and free-ranging European Bison from two different sites (0.003 μg/g) indicate a low mercury load [54]. The levels of such vital trace elements, such as iron, titanium, and vanadium, in the hair of a modern European bison are much lower than in hair of both prehistoric fossils [55].

#### 3.1.3. Woolly Rhino Hair

There was only one sample of woolly rhino hair, and it had a high Hg concentration of 0.092 ± 0.003 μg/g (second highest of all studied samples). This is supposedly related to Hg variability based on climatic stage and is discussed in detail in Section 3.4.

#### *3.2. Comparison with Modern Animals' Hair and Reference Intervals*

Unfortunately, we cannot compare Hg levels in modern animals to historical levels in the same animal species, because there are no modern animals identical to the mammoth mammals. Elephants are closest to mammoths genetically, but they have a different body and habitat. We compared the Hg concentration in hair of mammoth and yak, as they have similar characteristics (nutrition, long hair, etc.). A yak (*Bos mutus*) from the Barnaul Zoo has a low Hg level (0.006 ± 0.001 μg/g) (Table 4). There are no data about Hg concentration in the hair of yak from other regions of the world, but concentrations of other elements were found to be comparable between Altai and Asian yak, indicating that the differences in their exposure to metals are insignificant for the vast territories of their habitats [56,57]. Hg concentrations in all mammoth mammals' hair samples were significantly higher than Hg concentration in unpolluted hair samples of modern animals such as cattle (0.0066 ± 0.0002 μg/g) (Table 4). Methods for determining reference ranges in hair by using results from a large human population are described in detail elsewhere [11]. The reference intervals and 90% confidence intervals for the lower and upper limits were calculated for hair trace-element content in cattle (*Bos taurus*) per the recommendations of the American Society for Veterinary Clinical Pathology Quality Assurance and Laboratory Standard Guidelines [48,58]. Concentrations of Hg in mammoth fauna mammals' hair mostly lie within the optimal reference range for cattle (Table 4), excluding two highest concentrations,

which apparently reflect high environmental exposure to mercury in these mammals during the last periods of their lives. For plant-eating animals, vegetation is one of the main factors characterizing the living conditions [20], although Hg accumulation by animals depends both on their diet and habitat.

#### *3.3. Hg Levels in Arctic Animals and Humans (Historical and Modern)*

In historical samples of hair of human mummies of the Aleutian Islands (Alaska) dating 1450 AD, mean total mercury concentration (5.8 ± 0.9 μg/g) is comparable to the levels observed in hair of modern residents of the northern polar territories (Alaska, Canada, Faroe Islands) [59]. That confirms the main contribution of the traditional nutrition based on fish and meat of marine mammals to the accumulation of mercury for residents of these territories. The Egyptian, Chilean, and Peruvian mummies had mercury exposures below the US EPA reference level of 1 μg/g and were considerably lower than that of northern pre-industrial populations [60]. Hg concentrations in hair samples of historical (10.42 ± 1.31 μg/g) and modern (10.42 ± 2.45 μg/g) arctic foxes were similar and strongly correlated with ecotype and available food source [22]. Unlike humans and foxes, Hg concentrations in the hair of Greenland polar bears showed a significant increase from 0.52 to 4.9 μg/g (from 1300 to 2000 years) [61]. Comparison between Hg levels in the hair of the ancient dogs of the Seward Peninsula (0.657 ± 0.273 ng/g [62]) and the modern Alaska fish-fed dogs (0.54 ± 0.11 μg/g, [51]) did not show significant difference. Thus, when environmental exposure (atmosphere and water) to mercury is low, the increased levels of mercury in the bodies of ancient humans and animals are primarily associated with their diet. Mammals of the mammoth fauna have low Hg levels, since they are herbivores that get mercury from plants and accumulate it in their bodies (and hair). Biomagnification along the food chain (as seen in aquatic ecosystems and fish-eating animals and humans) is not observed.

#### *3.4. Comparison with Other Paleoarchive Data*

Environmental archives such as lake and marine sediments, peat bogs, glacial ice, and tree rings are widely used to reconstruct Hg accumulation at the local, regional, and global scale. All archives have their advantages and disadvantages, but none of them is a definite record of past mercury levels, because of the complexity of the mercury cycle's being influenced by various processes in each archive [63,64]. Most archives record the past several hundred to several thousand years (ice cores, peat bogs, lake and marine sediment cores, and tree rings), whereas long-term paleorecords recording up to a hundred thousand years are scarce (ice cores of Antarctica and Greenland, sediments cores and speleothems). Ice Age animals in this study lived 45 up to 10.5 kyr BP. Animals were exposed to mercury from the diet and the environment. The natural sources of Hg emissions were volcanoes, air–sea and soil–vegetation–air exchange, biomass burning (wildfires), and the revolatilization of deposited Hg from the soils (including release associated with permafrost thawing due to climate change).

The highest Hg concentration recorded in this study dates to 33.930 cal kyr BP. The other peak of Hg concentration at 23.292 cal kyr BP falls into the LGM period. An increase of mercury concentrations in hair coincides with variations of Hg concentrations and depositions recorded in the Antarctica Dome C ice core (Figure 4). Due to constant snow accumulation, Hg concentrations and fluxes change synchronically. Total Hg and Hg2<sup>+</sup> concentrations are also characterized by similar trends, except during the initial period from 15 to 2 kyr BP. The highest peak occurred during the Karginian interstadial of the Late Pleistocene, the period of maximum insolation in 200 ka years [65]. The presence of mosses in Malykhchinsky bison's food masses indirectly indicated a significant wetland area in this period because of significant climate warming. Hg increase in environment caused by rapid release of mercury during thawing periods was recorded in other paleoarchives, such as sediments in Limnopolar Lake (South Shetland Islands), where extraordinary high Hg enrichment was observed [65]. Research based on about 600 samples from soil permafrost cores (Alaska) discovered that the active layer is the largest Hg pool on the planet. The Northern Hemisphere permafrost region contains 1656 ± 962 Gg Hg, of which 793 ± 461 Gg is frozen in permafrost. The active layer and permafrost contain nearly twice as much Hg as all other soils, the ocean, and the atmosphere combined [66]. This allows us to assume that, in the past warm climatic periods, thawing permafrost caused significant mercury to be released into the environment, from the active layer of permafrost.

**Figure 4.** Concentrations (**a**) and fluxes (**b**) of total mercury (HgT) and inorganic mercury (Hg2<sup>+</sup>) in the EPICA Dome C ice core [67] and total Hg concentrations in mammoth fauna mammals' hair. Concentrations and fluxes of Hg below the Method Detection Limit (MDL) are presented in the graph as 1/2 of MDL.

The exogenous Hg in hair of ancient animals is mainly due to its sorption on the hair surface, from the atmosphere. The endogenous Hg in hair also can be due to Hg0 influx. Elemental mercury (up to 80% of inhaled Hg0 vapors) is absorbed in the lungs, quickly diffused into the blood, and distributed to all organs of the body; it also accumulates in growing hair.

The second peak corresponds to the Last Glacial Maximum, where the Antarctic record also shows a drastic increase in Hg concentrations during the LGM. It was found for mercury in the Antarctic [67] that the oxidation of gaseous mercury by sea-salt-derived halogens occurred in the cold atmosphere. The oxidized mercury compounds were then transferred to the abundant mineral dust particles and deposited. A significant correlation between dust concentrations and changes of temperature during glacial periods was confirmed by comparing dust and stable isotope, up to 90% of the dust variability can be explained by the temperature variations. The deposition of dust in Antarctica during glacial periods is about 20 times higher than during interglacials [68]. The cooling marked in three independently dated North Atlantic marine sediment cores is synchronous with the sharp increase in dust flux recorded in the Greenland ice cores, an increase in dust transport from Asia to Greenland observed during few Greenland stadials [69,70]. Deposition of mercury with dust on the surface of the land and its accumulation by the plants and snow might be the reason of increased dietary Hg exposure of herbivores. A high content of loamy particles in fecal samples of mammoths indicates an occasional or deliberate lithophagy [71]. Moreover, when thick ice completely covers the water, animals eat snow. Thus, changes of mercury concentrations in the hair of prehistoric animals are in good agreement with global changes of mercury concentrations recorded in other paleoarchives of Northern Hemisphere. It should be noted that the question remains open: Is there is a real difference in the deposition, distribution, and conservation of mercury in the Northern and Southern Hemispheres, or are Hg changes global? All archives preserve Hg differently and present changes in global Hg cycle at various spatial and temporal scales [64].

#### **4. Conclusions**

Mercury content in the hair of mammoths and other prehistoric animals allows us to estimate changing mercury levels between 40,000 to 10,000 years ago. Since the amount of ancient hair is very limited, we suggest comparing the hair of different animal species which are similar in diet and habits, as well as hair samples taken from different parts of the animals' body. All prehistoric animals have a low Hg level in their hair. This level is below concentrations associated with toxicity in wildlife and do

not exceed background levels of mercury in hair of non-seafood consumers (0.5 μg/g). Most of the Hg concentrations in the hair of prehistoric animals were within the reference range for modern cattle.

There are many advantages to using ancient hair as an indicator of environmental pollution, and now we present a new application of hair as an indicator of climatic changes. We hypothesize that Hg concentrations in hair reflect the variation in Hg level in the environment changing with climate changes, and can be used as a proxy for climate change assessment. The increase of Hg concentration in hair during the coldest climatic stages (such as LGM) coincides with the increase in Hg deposition on the Earth's surface, associated with the highest atmospheric dust loads. Moreover, mercury can be released to the atmosphere because of permafrost thawing during interstadial warming; the highest Hg concentration coincides with the Karginian interstadial of the Late Pleistocene, the period of maximum insolation and warming. Climate changes in warm and cold climatic stages were oscillatory; relatively warm periods alternated with cooler periods during each glacial and interglacial. For example, Karginian interstadial consisted of five periods (three warming and two cooling), in which features of the distribution and boundaries of permafrost are still under study.

Mammoth fauna mammals' hair, together with other natural archives, will be useful in assessing the response of Hg cycle to climate change. More paleo data are necessary to confirm our first finding, and to clarify whether these changes will differ for the Northern and Southern Hemispheres' archives, so we are planning further studies of mammoth fauna mammals on a wide spatiotemporal scale.

**Author Contributions:** This work was carried out in the framework of cooperation between Institute for Water and Environmental Problems from Barnaul, Diamond and Precious Metals Geology Institute, and Ammosov's North-Eastern Federal University from Yakutsk. Conceptualization and investigation, S.E. and G.B.; mammoth samples and morphological information, G.B.; methodology, S.E.; investigation and analysis, T.S., S.E., and M.S.; supervision T.P.; resources and funding acquisition, G.B. and T.P.; writing—original draft preparation, S.E.; writing—review and editing, S.E., G.B., and T.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was performed under a State Assignment of the IWEP SB RAS (project no. AAAA-A17- 117041210242-1). Studies by G.G. Boeskorov were performed under a State Assignment of the Institute of Diamond and Precious Metal Geology, Siberian Branch, Russian Academy of Sciences (project no. 0381-2019-0002); studies by M.V. Shchelchkova were performed under a State Assignment of the Ministry of Education and Science of the Russian Federation (project no. 37.7935.2017/6.7).

**Acknowledgments:** We are grateful to the Director of the Barnaul Zoo Sergey Pisarev for his help in obtaining yak hair samples for method testing. We expressed our deep gratitude to Eugenia Bakunova for kind help in English editing of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Mercury in Marine Mussels from the St. Lawrence Estuary and Gulf (Canada): A Mussel Watch Survey Revisited after 40 Years**

### **Daniel Cossa 1,\* and Anne-Marie Tabard <sup>2</sup>**


Received: 29 September 2020; Accepted: 22 October 2020; Published: 27 October 2020

**Abstract:** Various species of marine mussels have been used, in the last 50 years, as sentinel organisms for monitoring metal contamination along marine coasts. There are two main reasons for this: these mollusks concentrate metals in their soft tissue and they are geographically widespread. In practice, trace metal concentrations in mussel soft tissue reveal (after some correction for biotic effects) the contamination level of their surrounding environment. We present the results of a mercury (Hg) survey in *Mytilus* spp. collected in the summers of 2016, 2018, and 2019 at 51 stations distributed along the coasts of the Estuary and Gulf of St. Lawrence. Mercury concentrations ranged from 0.063 to 0.507 <sup>μ</sup>g g−<sup>1</sup> (dry weight, dw), with a grand mean of 0.173 <sup>±</sup> 0.076 <sup>μ</sup>g g−<sup>1</sup> dw (±1 standard deviation), and a median of 0.156 μg g−<sup>1</sup> dw for the 504 individuals analyzed. Mercury contents per individual mussel were significantly (*p* < 0.01) related to shell length and dry tissue weight, with the smaller individuals having the highest Hg concentrations. To take into account these biotic effects, we normalized Hg concentrations of the mussel soft tissue for constant shell length (L) and soft tissue weight (TW) based on the log-log relationships between Hg content and L or TW. The normalized Hg contents of mussels varied from 10.9 to 66.6 ng per virtual individual of 35 mm length and 0.17 g dry weight. A similar normalization procedure applied to 1977–1979 data, yielded a very similar range: 12 to 64 ng. This observation suggests that the Hg bioavailable to marine mussels in the study area did not change over a span of 40 years. Regional Hg distribution patterns indicate a gradual decrease of Hg content in mussels downstream from freshwater discharges to the St. Lawrence Estuary and the Baie des Chaleurs, suggesting that rivers constitute a significant Hg source in these estuarine systems. Atmospheric Hg deposition and concentration in marine waters of the Atlantic Ocean are known to have decreased in the last decades. However, in coastal environments, the response to these changes does not seem to be rapid, probably because of the long residence time of Hg in soils before being exported to coastal areas.

**Keywords:** mercury; mussel; mussel watch; *Mytilus*; St. Lawrence

#### **1. Introduction**

Marine mussels (*Mytilus* spp.) have been successfully used over the past fifty years as sentinel organisms for monitoring metal contamination along marine coasts [1–9]. The reason for this is that this mollusk genus is geographically widespread in sub-boreal and temperate environments and it concentrates metals in its soft tissue in proportion with the concentration in surrounding waters [10–12]. Thus, trace metal concentrations in the soft tissue of the blue mussel reveal the contamination level of the waters of its environment. Applied to the monitoring of temporal and geographical trends of chemical contamination, this approach has been named "Mussel Watch", and has been adopted as one of several coastal environmental quality-monitoring strategies by United Nations programs [3,13]. However, biological factors related to mollusk growth rate also control metal uptake and excretion and must be taken into account in order to optimize the use of mussels as sentinel organisms [4,10,14–16]. Several Mussel Watch programs have been carried out to monitor trace metal contamination along the eastern coasts of Canada and US [8,17–19], and especially mercury (Hg) contamination in the Estuary and Gulf of St. Lawrence [1,20,21]. These programs have not been maintained over more than a few years and, consequently, temporal trends are not as well documented. Nevertheless, they constitute baselines against which future assessments can be compared [22].

Here, we present the results of a Hg survey carried out in the summers of 2016, 2018, and 2019 in the St. Lawrence Estuary and Gulf, forty years after the first Mussel Watch was completed in this area. We used a normalization model to minimize the effects of biological factors on the Hg content of the mussel soft tissues. Our observations suggest that the amount of Hg bioavailable to marine mussels in the study area is similar to what it was 40 years ago. Sub-regional Hg distribution patterns indicate a gradual decrease of Hg content in mussels downstream from the main freshwater tributary (St. Lawrence River) to the St. Lawrence Estuary.

#### **2. Materials and Methods**

#### *2.1. Sampling and Pre-Treatment*

Adult specimens of *Mytilus edulis* ranging from 23 to 84 mm in length were sampled at 51 sites, at mid-tide level along the shores of the Estuary and Gulf of St. Lawrence (Figure 1). Geographical coordinates are given in Appendix A Table A1. Ten to twelve specimens (except for Station 36 with only 4 individuals) were collected at each site between 2–10 August 2016 along the South Shore of the Lower St. Lawrence Estuary (LSLE) and around the Gaspé Peninsula, between 27 August and 8 September 2018 along the North shore of the LSLE and Moyenne Côte Nord, and between 22–30 August 2019 along the Baie des Chaleurs, Northumberland Strait and Cape Breton (Figure 1). A few individuals with a shell deformity or green soft tissue were excluded to avoid the species *Coccomyxa*-infested *Mytilus trossulus*, for which background information about their trace element bioaccumulation properties is scarce [23–25]. In the field, mussel shell lengths were measured with a Vernier caliper, weighted, and their soft tissues freed from the shell. Soft tissues were kept at +4 ◦C in 2 mL of an ethanolic solution (90%, *v*/*v*) during the sampling periods, then kept at −18 ◦C until freeze-dried.

#### *2.2. Chemical Analyses*

The freeze-dried soft tissue of each mussel was individually analyzed using an automated atomic absorption spectrometer designed for Hg determinations (AMA-254, Altec, Czech Republic). The analytical procedure is described as EPA Method 7473 [26]. The analytical accuracy was checked every 10 determinations, using Certified Reference Materials (CRM): BCR-278R from the Institute for Reference Materials and Measurements, or IAEA-336 from the International Atomic Energy Agency, or DORM-4 from the National Research Council of Canada. BCR-278R is a mussel tissue, whereas IAEA-336 and DORM-4 are lichen and fish muscle, respectively. Analytical results were always within the target range of the certified values. Reproducibility, calculated as the variation coefficient (i.e., confidence interval/mean) of 6 replicate analyses of CRMs, varied between 2% and 5%. The detection limit, defined as 3 times the standard deviation of the 6 blank replicates, was 0.002 μg g−<sup>1</sup> (dry weight, dw). The amount of Hg dissolved in the ethanolic solution was measured in 20 samples taken at random; it never exceeded 2% of the total Hg burden in the soft tissues.

#### *2.3. Normalization Procedure and Statistics*

Mercury mussel content depends upon mussel shell length and soft tissue mass. The former because Hg accumulated throughout the lifespan of the mussel, the latter mainly because of its seasonal weight gains and losses [4,10]. To minimize the influence of mussel size on their Hg content, the normalization procedure described by Cossa and Rondeau [1] was applied. Normalization consists of correcting Hg raw data based on a multilinear regression function relating logHg content and shell length (logL) and soft tissue dry weight (logTW) (see Section 4.3). This normalization minimizes the influence of size on Hg content of the animals, thereby enabling the interpretation of mussel Hg content in terms of Hg bioavailability in their surrounding environment. According to Cossa and Rondeau [1] a two-fold Hg content distinction, environmentally-sound in terms of Hg bioavailability, can thereby be achieved.

The statistical analyses were performed with Xlstat software from Addinsoft (https://www. xlstat.com/). As their variance tends to increase with their mean, the Hg concentration values were log10-transformed before any statistical treatment.

#### **3. Hydrological and Ecological Settings of the Area of Study**

The studied area comprises two main regions, the LSLE and Gulf of St. Lawrence (Figure 1). The LSLE stretches from the mouth of the Saguenay Fjord (Station 26) to Pointe des Monts (Station 35); the south and north shores of the LSLE strongly differ in their hydrographic characteristics. The south shore is characterized by an estuarine circulation as the brackish waters, flowing out the St. Lawrence Estuary, run along the south shore to form the Gaspé Current, whereas the north shore of the LSLE is affected by upwellings that bring deep waters of the Gulf of St. Lawrence to the surface. The Gulf is an epicontinental sea open to North Atlantic water inputs by two straits: the Cabot Strait between Cape Breton and Newfoundland, and the Strait of Belle Isle between Labrador and Newfoundland. The Northeastern shore of the Gulf, called "Moyenne Côte Nord" (Figure 1), is under the hydrographical influences of upwellings (especially near the Anticosti Gyre), freshwater inputs from large rivers draining the Canadian Shield, and Labrador surface waters inflowing the Gulf westward through the Strait of Belle Isle. The estuarine waters of the LSLE that spread along the Gaspé Peninsula are diluted by more salty waters from the Anticosti Gyre (Figure 1). The shores of the Baie des Chaleurs, the Northumberland Strait, and the west coast of Cape Breton are characterized by relatively shallow waters and warmer water temperatures than in the northern part of the Gulf. In addition, the surface waters of the Baie des Chaleurs are impacted by river inputs, mainly from the Restigouche River, located at the western end of the Bay. A geographical partition of the Estuary and Gulf of St. Lawrence has been proposed based on ecological criteria [27]. These authors concluded that the most biologically significant hydrological features of the system are the LSLE, the Gaspé Current that hugs the coast of the Gaspé Peninsula, and the Northwestern Gulf, i.e., the Anticosti Gyre whose waters intercept the Moyenne Côte Nord west of Anticosti Island. These findings are supported by fluorescence distributions in the surface waters of the St. Lawrence Estuary and Gulf [28]. Based on the above-listed hydrological and ecological criteria, we decided to distinguish the following six regions: (1) North shore of the LSLE, (2) Moyenne Côte Nord, (3) South shore of the LSLE, (4) Gaspé Peninsula, and (5) Baie des Chaleurs, and (6) Northumberland-Cape Breton continuum.

**Figure 1.** Location of the sampling sites and hydro-ecological regions. LSLE: Lower St. Lawrence Estuary. P.-E.I. Prince Edward Island. Brown arrows indicate main surface water circulation.

#### **4. Results and Discussion**

#### *4.1. Geographical Patterns of the Growth*

Total shell length (L), width (l), height (h), and tissue mass (M) measurements on mussel allow us to calculate allometric growth indices, that may vary depending on the *Mytilus* species and growth conditions [29,30]. In the St. Lawrence Estuary and Gulf, L, l, and h are linearly correlated, whereas L and M are better related by power functions [17]. These findings are confirmed by the present data (Table 1, Figure 2). In this study (2016–2019), allometric growth indices (h/L, l/L ratios, and the parameters of the L vs M power functions) were similar to the earlier sampling (1977–1979) (Table 1). This similarity in the mussel morphometry suggests that the possible environmental changes in the St Lawrence Estuary and Gulf [31] during the last decades, such as temperature, did not generate shell shape differences despite the high plasticity of the mussel shell [32]. Nevertheless, the parameters of the L vs M equations, derived for each of the six regions, display some variations (Figure 2, Table 1). Such slight regional morphological shell differences are consistent with the small regional differences in the absolute growth rate estimated in 1977–1979 [17]. We can infer from these allometric index comparisons that mussel growth conditions have not differed significantly for the last forty years, with slight regional differences in the L *vs* M index persisting (Table 1). This finding allows us to compare the Hg load of the mussel soft tissues between the two time periods without any significant bias affecting Hg bioaccumulation due to the change of mussel growth rates.

**Figure 2.** Length-mass relationship of *Mytilus edulis* shell from six hydro-ecological regions of the St. Lawrence Estuary and Gulf (Figure 1). Mussels sampled in 1977–1979 (n = 1222) and in 2016–2019 (n = 504) periods.

**Table 1.** Allometric relationships in the mussel shell shape and their changes between 1977–1979 and 2016–2019 periods. L: length (cm), l: width (cm); h: height (cm); M: mass (g). All the relationships were statistically significant with a *p* < 0.001. AGR: Absolute growth rate (cm y<sup>−</sup>1) from ref. [17]. 1977–1979 statistics were calculated on 1222 individual mussels vs. 504 in 2016–2019.


#### *4.2. Estuary and Gulf of St. Lawrence Mussel Watch in the North Atlantic Context*

The Hg concentrations in mussel soft tissues ranged from 0.06 to 0.51 μg g−<sup>1</sup> dw (n = 504), with a grand mean of 0.17 <sup>±</sup> 0.08 <sup>μ</sup>g g−<sup>1</sup> dw (<sup>±</sup> 1 standard deviation), and a median of 0.16 <sup>μ</sup>g g−<sup>1</sup> dw. In 1977–1979, the Hg concentrations were very similar with a mean of 0.16 <sup>±</sup> 0.05 <sup>μ</sup>g g−<sup>1</sup> dw, and a median of 0.15 μg g−<sup>1</sup> dw (n = 442) [33]. Very few Mussel Watch programs have been carried out uninterrupted for longer than a few years with the notable exception of the US and French monitoring programs that have been running yearly since the 1980s [34,35]. The median concentration in the St. Lawrence system is slightly higher than the most recent available medians: 0.11 μg g−<sup>1</sup> dw (2005–2012, n = 298) for the US Mussel Watch and 0.12 μg g−<sup>1</sup> dw (2000–2004, n = 303) for the French Mussel Watch programs [34,35]. It is similar to the annual median concentrations of the French coasts

Mussel Watch for the period 1980–94 (0.15 μg g−<sup>1</sup> dw, n = 937) [36,37]. In summary, the current median Hg concentrations of the St. Lawrence Mussel Watch is similar to the long-term representative surveys carried out along the Atlantic shores in the last 40 years.

#### *4.3. Temporal and Regional Hg Trends*

To reach an optimal discriminating capacity for detecting temporal and regional trends, we used a normalization procedure minimizing bias due to biological factors (see above). The relationships between the Hg content of the 504 analyzed mussels and the dry mass of their soft tissue (TW) or shell length (L) were highly significant (*p* < 0.001). The Equation used for the normalization model built with the entire sampling set is:

$$\log \text{Hg} = 1.39 \pm 0.15 \log \text{L} + 0.49 \pm 0.04 \log \text{TW} - 0.36 \text{ (R}^2 = 0.87; \text{n} = 504) \tag{1}$$

compared to:

$$\text{logHg} = 1.43 \pm 0.49 \,\text{logL} + 0.34 \pm 0.16 \,\text{logTW} - 0.48 \,\text{(R}^2 = 0.91; \text{n} = 143) \tag{2}$$

obtained for the St Lawrence Mussel Watch performed in 1977–1979 [1]. Regression coefficients, obtained forty years apart, are not statistically different (*p* < 0.01).

Normalized Hg content (individual of 35 mm length and 0.17 g dry weight) varied from 10.9 to 66.6 ng per individual, with a grand mean of 42.0 ± 2.1 ng (±1 standard deviation), and median of 44.2 ng, for the 504 individuals. Expressed as Hg concentration in the soft tissue, the range becomes 0.06–0.39 μg g−<sup>1</sup> (dw). A similar normalization procedure applied to the 1977–1979 data, gave a very similar range: 12 to 64 ng per individual mussel [1] for the same stations. This finding strongly suggests that the Hg bioavailable to marine mussels in the study area has not changed over 40 years. At first glance, these observations are surprising since a decrease in Hg inputs to the St. Lawrence Estuary and Gulf waters could be expected (see below).

First, riverine Hg inputs to the St. Lawrence Estuary should have decreased due to the implementation of the International regulation on the Great Lakes basin [38–40]. Few time-series observations of Hg inputs from riverine sources, but a high-frequency sampling experiment of water was performed over 18 months in 1995–96 in the St. Lawrence River. Mercury export to the LSLE was estimated at ~1.2 Mg y<sup>−</sup>1, for a mean dissolved Hg concentration of 0.60 <sup>±</sup> 0.46 ng L−<sup>1</sup> [41]. Since that period, no systematic study has been published on temporal variations of Hg concentrations in the St. Lawrence River. According to a recent government report, the Hg flux to the St. Lawrence River has not changed between 1995–1996 and 2004–2008 [42]. A comparison of mussel mean Hg contents at stations from the head of the LSLE (Stations 53–56), where freshwater influence is maximum, fails to show a statistically significant difference (*p* < 0.01) in Hg tissue levels between 1977–1979 (44.6 ± 6.9 ng, n = 5) and 2016–2019 (37.4 ± 6.6 ng, n = 5). It could be interpreted to imply that Hg bioavailability in freshwaters of the LSLE has not changed in the last 40 years.

Secondly, atmospheric Hg concentrations and wet deposition in Eastern North America are thought to have declined during the 1990–2010 period [43–45]. Regulations for reducing Hg emissions were implemented in New England and Eastern Canada [46]. In an assessment of Hg sources and fate in a marine environment very close to the Gulf of St Lawrence (The Gulf of Maine), Sunderland et al. [47] reported that: "Temporal patterns in sentinel species (mussels and birds) have in some cases declined in response to localized point source mercury reductions but overall Hg trends do not show consistent declines". Likewise, studies of coastal Massachusetts, New Hampshire, and Maine, reveal that no significant temporal trends in mussel Hg concentrations were found between 1990 and 2010 at 12 of the 15 monitored stations [19]. Our results are also consistent with those of Hg trends in herring gull eggs from Atlantic Canada collected between 1972 and 2008 [48]. These authors reported that, after adjusting Hg trends for dietary shifts, environmental Hg in coastal ecosystems had remained relatively constant in Eastern Canada over the previous 36 years. According to Sunderland et al. [43], reductions in atmospheric Hg deposition from North American sources could have been offset by increased deposition from global Hg sources.

Figure 3 illustrates the station-to-station variability in mussel tissue concentration. Strikingly, the continuum on the South shore of the LSLE- Gaspé Peninsula (Stations 56 to 14) exhibits a gradual Hg decrease from the brackish estuarine water originating from the St Lawrence Estuary to the Gaspesian coast as they are diluted by the Gulf waters originating from the Anticosti Gyre. This pattern strongly suggests that St. Lawrence River waters are a significant source of bioavailable Hg. Seaward, along the Gaspé Current, mussels exhibit very low Hg contents and variability (Stations 48, 49, 16 to 24, Figures 1 and 3). This would imply that Gaspé Current waters are impoverished in bioavailable Hg compared to waters of the LSLE. Increasing Hg reduction and evasion in the atmosphere in the productive waters (Anticosti Gyre and Gaspé Current, see Section 3) may favor this Hg depletion. It should be noted that the Gaspesian coastline is also an area where riverine inputs are small compared to other LSLE and Gulf shores. Another Hg dilution structure is visible along the Baie des Chaleurs (Stations 6 to 3, Figures 1 and 3) where riverine inputs are also important (see Section 3). Interesting to note is that, during the 1077–79 Mussel Watch, high Hg contents in mussel soft tissue were also observed at stations along the southern shores of the LSLE and of the Baie des Chaleurs where brackish waters are present (see Figure 4 in Ref. [1]).

**Figure 3.** Mercury content (ng) of individual mussels after normalization to 35 mm-long and 0.17 g virtual specimens. Mean contents are indicated by a horizontal bar and the median by a cross. Each station combines 10 to 12 individual mussels. Stations are eastward distributed on the x-axis.

The North Shore of the LSLE and the Moyenne Côte Nord coastlines display inter-station large variations in mussel tissue Hg content, even within the same oceanographic entity (e.g., Moyenne Côte Nord, Figures 1 and 3). These fluctuating distributions are difficult to interpret since no steady geographical trend can be observed in Hg content. Nevertheless, some hypotheses can be put forward to explain these. Mercury content in the coastal waters of the North Shore of the St. Lawrence Estuary and Gulf are most likely influenced by upwellings and organic-rich river waters that drain through the Canadian Shield soils [49,50]. The organic matter of these circumneutral pH waters contains strong Hg-binding functional groups (e.g., thiol) that favor Hg transport in solution [41,51–53]. On the other hand, upwellings bring saline low Hg waters of Atlantic origin to the surface, especially in the Minganie region (Stations 44–47, Figures 1 and 3) [54]. These Atlantic waters enter the Gulf through

Cabot Strait and are known to contain low picomolar Hg levels [55]. Conversely, waters entering the Northeastern Gulf through the Strait of Belle Isle originate from the Labrador Sea Current, known for its relative high Hg concentrations due to the inputs of the organic-rich rivers outflowing the Arctic Canadian Archipelago [56]. These waters, brought by the Labrador Sea Current, are transported westward along the J. Cartier Passage, a remote region without any known trace-metal sources but with relatively high Hg content in sampled mussels (Figures 1 and 3).

It is interesting to note that low Hg content mussels are located along the Gaspé Peninsula (Stations 48, 49, 16–24) and Minganie (Stations 44–47) coasts, where commercial or experimental mussel farms are installed (mapaq.gouv.qc.ca/fr/Peche/aquaculture) [57,58]. This latter area is known for its low-temperature upwellings in summer [54].

#### *4.4. Current and Former Hot Spots*

Industrial areas distributed along the Estuary and Gulf coastlines are likely to shelter polluted sites. These would include a former chlor-alkali plant, pulp and paper mills, and smelters [33]. The highest Hg concentrations of the present survey were found at Station 32 near the city of Baie Comeau. The city hosts an active industrial complex and the bay sediments are enriched in trace metals from current and legacy sources [59,60]. Our study fails to highlight other potential hotspots that had been identified by previous Mussel Watch Programs. A chlor-alkali plant located near the mouth of the Restigouche River in the Baie des Chaleurs has been shown to generate local Zn and Hg contamination [61], but its impact on mussel tissue Hg contents was not revealed by our study. This might be explained by may be the result of the 20-km distance between Station 6 and the plant and/or the fact that it was shut down in 2008 [62]. Likewise, the previously reported Hg contamination [20] at the mouth of the Saguenay Fjord (Station 26), is no longer being observed, as the main Hg source, the Arvida chlor-alkali plant was shut down 40 years ago [63].

#### **5. Summary and Conclusions**

A Hg Mussel Watch Program was carried out at 52 stations distributed along the intertidal zone of the Estuary and Gulf of St. Lawrence, in the summers of 2016, 2018, and 2019. This survey took place forty years after a previous Mussel Watch survey was performed in this area (1977–1979), at almost all stations currently sampled. The same normalization procedure applied 40 years ago was also applied to the current data set. It allows the minimization of the effects of biological factors on the Hg content of the mussel soft tissues and optimizes the interpretation of Hg distribution in terms of the availability of Hg in coastal waters. Current Hg concentrations in the soft tissue of mussels from the St. Lawrence System are similar to those of the perennial Mussel Watch Programs implemented along the US and French Atlantic coasts [34,35]. The results indicate that the normalized Hg contents of mussels in the Estuary and the Gulf of St Lawrence were similar at sampling times 40 years apart. In addition, sub-regional Hg distribution patterns indicate a gradual decrease of Hg content in mussels downstream from freshwater inputs, which suggests that the spatial distribution of Hg concentrations in the soft tissues of marine mussels is, in part, governed by riverine Hg sources. Atmospheric Hg deposition and concentration in marine waters of the Atlantic Ocean are known to have decreased in the last decades [44,64]. However, in coastal environments, the response to these changes does not seem to be rapid. This probably results from the long residence time of Hg in soil and land cover before being exported to coastal areas with freshwater discharges. Indeed, different environmental reservoirs have different time scales as regards Hg mobility [65,66]. Furthermore, coastal sediment resuspension/sedimentation cycling may also contribute to a longer retention time of Hg in the intertidal marine environments, making mussel habitats conducive to the retention of legacy Hg. The changes in atmospheric Hg inputs are thus damped in terms of Hg availability for sessile benthic animals. A renewed Mussel Watch survey in a few decades may elucidate the long-time trend of Hg load in the St Lawrence Coastal System.

**Author Contributions:** D.C. designed the experiment, performed sampling and analyses, and wrote the article. A.-M.T. participated in sampling and sample preparation. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** Thanks are due to A. Mucci for providing some of the sampling material and S. Guédron for his support in mercury analyses performed within the analytical chemistry platform of ISTerre (OSUG-France). Special thanks to E. Bourget, A. Mucci, J.-M. Sévigny, and M. Vautour for their comments on the manuscript.

**Conflicts of Interest:** The authors declare that there is no conflict of interest regarding the publication of this paper.

#### **Appendix A**

**Table A1.** Station coordinates. See also map in Figure 1. Mean normalized Hg content of (±1 SD) in ng per individual.



**Table A1.** *Cont*.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

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

**Mark E. Brigham 1,\*, David D. VanderMeulen 2, Collin A. Eagles-Smith 3, David P. Krabbenhoft 4, Ryan P. Maki <sup>5</sup> and John F. DeWild <sup>4</sup>**


**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 fish 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 (THgaq), methylmercury (MeHgaq), 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 flat trends since 2009. In the four VNP lakes, epilimnetic MeHgaq concentrations declined by an average of 44 percent and THgaq 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 MeHgaq concentrations; however, biotic THg concentrations declined significantly in only one lake. Epilimnetic MeHgaq 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 benefitted aquatic resources, even in the face of global increases in mercury emissions.

**Keywords:** mercury; methylmercury; lakes; wet deposition

#### **1. 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.

**Citation:** Brigham, M.E.; VanderMeulen, D.D.; Eagles-Smith, C.A.; Krabbenhoft, D.P.; Maki, R.P.; DeWild, J.F. Long-Term Trends in Regional Wet Mercury Deposition and Lacustrine Mercury Concentrations in Four Lakes in Voyageurs National Park. *Appl. Sci.* **2021**, *11*, 1879. https://doi.org/ 10.3390/app11041879

Academic Editor: Stéphane Guédron

Received: 31 October 2020 Accepted: 7 February 2021 Published: 20 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Although deposited primarily as inorganic mercury, within aquatic ecosystems some fraction of the total mercury load is converted to methylmercury [3]. Mercury in fish tissues is predominantly in the methylmercury form [3]; thus, factors that control methylmercury production in aquatic ecosystems are of interest. Two such factors, relevant to the current paper, are the availability of inorganic mercury and sulfate. Increased mercury loading to aquatic ecosystems results in a proportional increase in methylmercury production and accumulation by aquatic biota [4]. Furthermore, because sulfate reducing bacteria are important as methylators of mercury [5,6], elevated sulfate deposition also can result in increased methylmercury production in aquatic ecosystems [5,7]; it is likely, therefore, that methylmercury concentrations in aquatic ecosystems (including fish) are not only greater than background, preindustrial levels, but greater than one might expect with solely a proportional response to mercury inputs.

With an increased understanding of the global extent of atmospherically driven mercury contamination in aquatic ecosystems, the U.S. and Canada have undertaken significant efforts to reduce anthropogenic mercury emissions. Lake-coring evidence indicates that mercury inputs to remote, mid-continental lakes peaked in the 1960s–1970s and declined subsequently [8]. More recently, emissions reporting has shown substantial declines in U.S. mercury emissions of approximately 87 percent from 1990–2017 [9], with the largest decreases occurring in the mid 1990's following regulations limiting emissions from hospital and municipal waste incinerators and chlor-alkali facilities; similarly large declines in mercury emissions occurred in Canada over the same time frame [10] (Table 1).

**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<sup>−</sup>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−<sup>1</sup> in 1980, to about 6 kg ha−<sup>1</sup> in 2000, to about 2.4 kg ha−<sup>1</sup> 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.

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–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.

#### **2. Materials and Methods**

#### *2.1. 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 km2 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–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

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 (MeHgaq) and total mercury (THgaq) from 2000 (Shoepack Lake) or 2001 (Brown, Ryan, and Peary lakes) through 2018; and trends in biotic THg as described below.

**Figure 1.** Location of four study lakes within Voyageurs National Park, and nearby National Atmospheric Deposition Program precipitation-monitoring sites MN16 and MN18.

**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<sup>−</sup>1, microequivalents per liter; TOC, total organic carbon; mg L<sup>−</sup>1, milligrams per liter].


#### *2.2. 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 (Hgprecip) 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 Hgprecip 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 Hgprecip and precipitation volumes were analyzed using a multiple-regression

approach developed previously [31], using the REG procedure in SAS software (v. 9.4. SAS Institute, Carey, NC, USA). Hgprecip 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[\text{Hg}\_{\text{precip}}] = \beta\_0 + \beta\_1 \log(\text{Precip}) + \beta\_2 \sin(2\pi \text{T}) + \beta\_3 \cos(2\pi \text{T}) + \beta\_4 \text{ T} + \varepsilon \tag{1}$$

where Hgprecip is mercury concentration in weekly precipitation samples in ng L<sup>−</sup>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[Hgprecip], 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, Hgprecip 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 Hgprecip 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 Hgprecip in precipitation. One outlier in particular had an abnormally high Hgprecip for a relatively large precipitation event and occurred late in the study period (site MN18; 5 July 2017; reported Hgprecip = 89.59 ng L−1; Precip = 27.686 mm). The resultant calculated weekly Hg deposition was 2.48 μg m<sup>−</sup>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 Hgprecip were set to missing. After removal of outliers, regression analysis was applied (again using Equation (1)), and missing log[Hgprecip] values were predicted by the regression. Predicted Hgprecip was calculated by exponentiating regression-predicted log[Hgprecip] 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 Hgprecip values, then, consisted of MDN-reported values in nearly all cases, but regression-predicted values where weekly Hgprecip 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 Hgprecip were calculated as annual deposition rate divided by annual precipitation depth. This methodology differs from NADP/MDN network methodology. NADP/MDN effectively estimates missing Hgprecip as the precipitation-weighted annual mean concentration, whereas our methodology estimates missing Hgprecip 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 Hgprecip 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.

#### *2.3. 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 (THgaq) 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 (MeHgaq) 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 (Me200Hg) 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 MeHgaq 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].

#### *2.4. 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.

#### *2.5. 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 MeHgaq sampling.

#### **3. Results**

#### *3.1. 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 Hgprecip 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).

**Figure 2.** Trend plots of annual wet mercury (Hg) deposition rates for MDN sites MN16 (**A**) and MN18 (**B**); and precipitationweighted annual mean Hg concentrations for MDN sites MN16 (**C**) and MN18 (**D**) (Trend line is locally weighted regression; shaded region is 95-percent confidence interval. LOESS smoothing parameters are as follows: (**A**) 0.833; (**B**) = 1.0; (**C**) 1.0; and (**D**) 1.0.

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−<sup>1</sup> 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 Hgprecip.

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].

**Table 3.** Trend analysis results for wet mercury (Hg) deposition rate and precipitation-weighted mean concentrations for two National Atmospheric Deposition/Mercury Deposition Network sites in northern Minnesota. [Hg deposition μg m<sup>−</sup>2, annual Hg deposition rate in micrograms per square meter; Precip.-weighted [Hg], ng L<sup>−</sup>1, annual precipitation-weighted mean Hg concentration in nanograms per liter; EV1998 and EV2018 are expected values for 1998 and 2018, respectively, from locally weighted regression of the variable of interest versus time in years. % change, percentage change in expected values from 1998–2018. Lower and upper 95 percent confidence limits shown in brackets.].


#### *3.2. 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 MeHgaq 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 MeHgaq (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 MeHgaq. The magnitude of the declines in MeHgaq 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).

**Table 4.** Epilimnetic lake water trends in methylmercury (MeHgaq), total mercury (THgaq), and sulfate, from locally weighted regression analysis, for the period 2001–2018. (EV2001 and EV2018 are expected values for concentrations of MeHgaq and THgaq, in ng L−1; and sulfate, in mg L−1, from locally weighted regression analysis, for 1 July 2001, and 1 July 2018, respectively. % change, percentage change in expected values from 2001–2018. Lower and upper 95-percent confidence limits shown in brackets).


**Figure 3.** Time series plots for methylmercury (MeHgaq) in epilimnetic lake water. Smooth lines are locally weighted regression lines; smoothing parameters, selected as described in the methods section: (**A**) Brown Lake, smoothing parameter = 0.713; (**B**) Peary Lake, smoothing parameter = 0.424; (**C**) Ryan Lake, smoothing parameter = 0.489; and (**D**) Shoepack Lake, smoothing parameter = 0.633. Gray shading indicates the 95-percent confidence interval. A high outlier for Brown Lake (reported concentration 0.51 ng L−<sup>1</sup> for 30 August 2012) is omitted from the plot and regression.

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 THgaq in Ryan Lake appears to be a significant decline with clear separation of confidence intervals in the beginning versus end of the study period.

**Figure 4.** Time series plots for total mercury (THgaq) in epilimnetic lake water. Smooth lines are locally weighted regression lines; smoothing parameters, selected as described in the methods section: (**A**) Brown Lake, smoothing parameter = 0.713; (**B**) Peary Lake, smoothing parameter = 0.641; (**C**) Ryan Lake, smoothing parameter = 0.424; (**D**) Shoepack Lake, smoothing parameter = 0.3. Gray shading indicates the 95-percent confidence interval.

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].

**Figure 5.** Time series plots for sulfate in epilimnetic lake water. Smooth lines are locally weighted regression lines; smoothing parameters, selected as described in the methods section: (**A**) Brown Lake, smoothing parameter = 0.711; (**B**) Peary Lake, smoothing parameter = 0.75; (**C**) Ryan Lake, smoothing parameter = 0.711; and (**D**) Shoepack Lake, smoothing parameter = 0.70. Gray shading indicates the 95-percent confidence interval.

#### *3.3. Lake Level Fluctuations*

MeHgaq correlates modestly (R<sup>2</sup> = 0.36) with lake-level anomaly for Brown Lake (Figure 6). However, none of the other lakes in this study showed significant correlations between MeHgaq and lake-level anomaly (not shown).

#### *3.4. Mercury Trends in Lake Biota*

Total mercury concentrations in age-1 yellow perch were well correlated with Aeshnidequivalent dragonfly larvae THg concentrations (R<sup>2</sup> = 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.

**Figure 7.** Total mercury (THg) concentrations in age-1 yellow perch in relation to dragonfly larvae (Aeshnid equivalent) THg concentrations. Shaded area shows the 95-percent confidence interval. [Regression equation: ln(THgYPE) = <sup>−</sup>1.756 + 1.582 <sup>×</sup> ln(THgAE), r<sup>2</sup> = 0.66, *<sup>p</sup>* < 0.0001, N = 40, where ln is the natural logarithm; THgYPE, total mercury in age-1 yellow perch equivalent; THgAE, total mercury in Aeshnid equivalent.].

For the three lakes with long-term biomonitoring, temporal patterns in biotic THg concentrations were similar to patterns in MeHgaq 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 MeHgaq, 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 MeHgaq, 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 MeHgaq concentrations during the same time period.

**Figure 8.** Total mercury (THg) in age-1 yellow perch versus time in (**A**) Brown Lake, (**B**) Peary Lake, and (**C**) Ryan Lake. Solid symbols represent geometric mean concentrations from yellow perch, open symbols indicate that yellow perch concentrations are derived from dragonfly larvae using the linear regression shown in Figure 7. Shaded area shows the 95-percent confidence interval. [Smooth lines are locally weighted regression. Smoothing parameters: (**A**) Brown Lake = 0.75; (**B**) Peary Lake = 0.83; and (**C**) Ryan Lake = 0.76. Gray shading indicates the 95-percent confidence interval.

#### **4. 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 MeHgaq concentrations declined by an average

of 44 percent and THgaq 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 MeHgaq and THgaq 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 MeHgaq 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 MeHgaq in a 2001–2002 survey of small lakes in Voyageurs National Park [21], could be responsible for the increase in MeHgaq in Brown Lake during the 2001–2012 time period. The apparent trend reversal, i.e., the decline in MeHgaq in Brown Lake that occurred starting around 2010, coupled with the correlation between lake level and MeHgaq, supports the hypothesis.

Higher observed MeHgaq 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 MeHgaq production within a lake's watershed, irrespective of loading from upstream lakes. However, the remaining lakes in this study showed no correlation between MeHgaq 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 MeHgaq or THgaq for Brown and Peary Lakes but not Ryan Lake. Whereas MeHgaq 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 MeHgaq or THgaq 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 MeHgaq 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 MeHgaq, compared to epilimnetic THgaq is likely driven by both decline in wet Hg deposition (and thus in-lake THgaq), 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 MeHgaq declines, in comparison to declines in THgaq, are consistent with the notion that MeHgaq 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.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2076-341 7/11/4/1879/s1, Table S1: Outlier samples identified by regression analysis of log[Hgprecip] versus logPrecip, seasonal terms, and time. Figure S1: Annual wet sulfate deposition for the National Atmospheric Deposition Network site at Marcell, Minnesota. Figure S2: Precipitation volume trend plots for National Atmospheric Deposition Network/Mercury Deposition Network sites Marcell, Minnesota (MN16) and Fernberg, Minnesota (MN18).

**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.S. Geological Survey (USGS). Sampling of yellow perch was funded by the NPS Great Lakes Inventory and Monitoring Network under Task Agreement J2105080012 of the Great Lakes-Northern Forest Cooperative Ecosystem Studies Unit and by the Great Lakes Restoration Initiative, Environmental Protection Agency Project Number 222, under Task Agreement J2105100001 of the Great Lakes-Northern Forest Cooperative Ecosystem Studies Unit—both under Cooperative Agreement H6000082000 between the NPS and the University of Minnesota. The USGS/NPS Water Quality Partnership supported lake-water data collection in 2001–03. The USGS Toxic Substances Hydrology program supported the analysis of samples and staff time at the USGS Mercury Research Laboratory to process the samples, archive the data, and assist in the data interpretation.

**Institutional Review Board Statement:** Previous fish research incorporated into this study was conducted in accordance with Animal Care and Use Procedure (ACUP 306.02) published by Cornell University and adopted by researchers from the University of Wisconsin-La Crosse and the U.S. National Park Service.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data used for this study are available, as detailed in the Methods section.

**Acknowledgments:** Since 2008, lake sampling was conducted by National Park Service staff, particularly Jaime LeDuc, with Richard Damstra and Alex Egan assisting with data management. In addition, since 2008, mercury analyses of lake-water samples have been provided by the U.S. Geological Survey's Mercury Research Laboratory (MRL); Jacob Ogorek (USGS/MRL) is acknowledged for laboratory support. Assistance with the sampling and analyses of fish and dragonfly larvae was provided by University of Wisconsin at La Crosse students including Sean Bailey, Matthew Brantner, Adam Hyer, John Kalas, Kevin Miller, Kristen Thoreson, Jeffrey Ziegeweid; and faculty Roger Haro, Mark Sandheinrich, and James Wiener (retired). We also thank the Dragonfly Mercury Project team (Colleen Flanagan Pritz, Sarah Nelson, James Willacker, Amanda Klemmer, Meghan Hess, Colleen Emery) and USGS Contaminant Ecology Lab (Branden Johnson, John Pierce) for laboratory support and data management. Jennifer Murphy (USGS) offered helpful advice on statistics. We thank Austin Baldwin (U.S. Geological Survey) and three anonymous reviewers for helpful reviews of earlier versions of this manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

