A National-Scale Assessment of Mercury Bioaccumulation in United States National Parks Using Dragonfly Larvae As Biosentinels through a Citizen-Science Framework.

We conducted a national-scale assessment of mercury (Hg) bioaccumulation in aquatic ecosystems, using dragonfly larvae as biosentinels, by developing a citizen-science network to facilitate biological sampling. Implementing a carefully designed sampling methodology for citizen scientists, we developed an effective framework for a landscape-level inquiry that might otherwise be resource limited. We assessed the variation in dragonfly Hg concentrations across >450 sites spanning 100 United States National Park Service units and examined intrinsic and extrinsic factors associated with the variation in Hg concentrations. Mercury concentrations ranged between 10.4 and 1411 ng/g dry weight across sites and varied among habitat types. Dragonfly total Hg (THg) concentrations were up to 1.8-fold higher in lotic habitats than in lentic habitats and 37% higher in waterbodies with abundant wetlands along their margins than those without wetlands. Mercury concentrations in dragonflies differed among families but were correlated (r2 > 0.80) with each other, enabling adjustment to a consistent family to facilitate spatial comparisons among sampling units. Dragonfly THg concentrations were positively correlated with THg concentrations in both fish and amphibians from the same locations, indicating that dragonfly larvae are effective indicators of Hg bioavailability in aquatic food webs. We used these relationships to develop an integrated impairment index of Hg risk to aquatic ecosytems and found that 12% of site-years exceeded high or severe benchmarks of fish, wildlife, or human health risk. Collectively, this continental-scale study demonstrates the utility of dragonfly larvae for estimating the potential mercury risk to fish and wildlife in aquatic ecosystems and provides a framework for engaging citizen science as a component of landscape Hg monitoring programs.

Introduction

Environmental mercury (Hg) contamination is recognized as a global health threat.[1,2] The severity and scope of environmental and human health risks posed by Hg contamination have motivated landscape-level assessments of variability in mercury bioaccumulation and its drivers[3,4] as well as efforts to minimize fish, wildlife, and human Hg exposure through reductions in environmental Hg releases.[5,6] Biosentinels are an important tool for both landscape scale assessments and effectiveness evaluations of Hg reductions because inorganic Hg loading is often decoupled from methylmercury (MeHg) production, uptake, and biomagnification through food webs.[7−9] However, implementing a biosentinel network at broad scales can be complicated by the availability of the appropriate taxa, and the cost and logistics of executing an appropriate sampling design.

Many animals can serve as effective Hg biosentinels, but their utility is dependent upon both the characteristics of the organism and the inference needs of the data.[10] Fish and aquatic-dependent birds are commonly used as Hg biosentinels because they can directly inform exposure risk to humans or other sensitive wildlife species that rely on aquatic food webs.[11] However, challenges can emerge when consistent taxa are unavailable across the landscape or area of interest,[7] movements obfuscate the spatial and temporal origins of bioaccumulated Hg;[12] or water bodies lack resident fish populations.[13] Vertebrate sampling can also be impeded by logistical, regulatory, or ethical concerns. Therefore, a geographically widespread biosentinel that occurs in a variety of habitat types would be particularly valuable if it provided a reliable measure of MeHg bioaccumulation at fine spatial scales, informed potential wildlife exposure in both fish-containing and fishless environments, and was logistically easy to sample.

Some aquatic insect larvae meet many of the above characteristics, but their widespread use has been limited because they can lack nexus to human exposure, may require technical expertise to identify, and exhibit low and variable percentages of methylmercury (% MeHg).[14] Dragonfly larvae (Odonata Anisoptera) may represent an exception in many respects. Dragonfly larvae are obligate predatory invertebrates and occupy diverse freshwater habitats across six continents,[15] have a narrow trophic range, and have tissues in which most Hg content is as MeHg.[16−18] Dragonfly larvae are common even in waters where other Hg biosentinels, such as fish, do not occur, and their high site fidelity ensures that their tissues reflect the localized food web MeHg availability.[15] Dragonfly larvae total mercury (THg) is well-correlated with both aqueous MeHg and sportfish THg,[17−19] linking dragonflies to both environmental Hg concentrations and potential human exposure pathways. Moreover, dragonfly larvae can serve as vectors of aquatic-derived MeHg to terrestrial predators after dragonfly emergence.[20] Finally, dragonfly larvae are relatively easy to collect and identify, making them ideal for nontraditional monitoring networks, such as citizen science.

Within a citizen-science framework we used dragonfly larvae as biosentinels for an assessment of Hg bioaccumulation in aquatic ecosystems of National Park Service units and other protected lands (hereafter NPS units). As protected environments generally free from Hg point sources, NPS units serve as optimal locations for evaluating the effectiveness of global Hg reduction efforts on Hg exposure in freshwater ecosystems. The goals of this effort were to (1) implement a national-scale biological Hg monitoring program in protected lands, (2) assess the geographic variation in dragonfly THg concentrations across habitat types, and (3) establish dragonflies as an effective biosentinel to inform potential exposure risk at a continental scale. We developed strict protocols and training approaches that facilitated engagement of citizen-science volunteers for sample collection in NPS units across the contiguous United States (US), Alaska, and Hawaii. We then examined the influence of intrinsic and extrinsic factors on dragonfly THg concentrations, and finally correlated dragonfly Hg concentrations with those from fishes and amphibians sampled from the same locations to determine the efficacy of dragonflies for informing Hg exposure and potential impairment in other taxa.

Methods

The Dragonfly Mercury Project is a national-scale research and monitoring effort that engages citizen scientists paired with trained NPS staff to collect dragonfly larvae for Hg analysis at national parks across the United States (Supporting Information).

Field Sites

During 2009–2018, dragonfly larvae representing six families were collected from 457 unique locations across 100 NPS units (Figure S1). We targeted at least three sites per unit, unless there were insufficient sites in the case of small park units; the average number of sites per unit was 4.4. In total, we compiled data from 877 unique site–year combinations. Site selection was not random, but we strived to incorporate a variety of diverse habitats within each park unit to assess habitat effects of Hg bioaccumulation.

Sampling sites were classified as one of five primary habitats (rivers, streams, lakes, ponds, and wetlands), with these primary habitats further divided into one of 16 secondary habitats based on water permanence, proximity to hydrologically connected wetlands, and waterbody size (Table S1).

Sample Collection

Dragonfly larvae

Dragonfly larvae were sampled from each site with dip nets, generally between spring and fall. Sampling was conducted by a combination of researchers, trained NPS staff, and volunteer citizen scientists under the supervision of trained project staff, following strict sample handling and storage protocols (https://www.nps.gov/articles/dragonfly-mercury-project.htm). Between 15 and 20 dragonfly larvae were targeted from each site and event, and efforts were made to ensure that replicates for each dragonfly family were sampled when present (Supporting Information). Each sample was double-bagged in prelabeled, polyethylene zipper-seal bags and held on wet or dry ice until transferred to a freezer within 8 h. To minimize potential contamination, all dragonfly larvae were handled with gloved hands, or clean plastic spoons. We stored all samples at −20 °C or colder until processing and Hg determination could be completed.

Fish

Fish samples were collected from a subset of the sites also sampled for dragonfly larvae. Some sites represented locations not associated with national parks but where dragonfly larvae and fish were both sampled together. From 41 unique locations, we collected several fish species representing four distinct guilds, piscivores (non-salmonid), sunfish, salmonids, and forage fish, that occur across the various habitats sampled (Table S2). Fish were collected using various methods, including electrofishing, hook and line, beach seines, and gill nets. Fish were euthanized under approved animal care and use protocols, placed into clean and labeled polyethylene bags, and stored on ice in the field until returned to the laboratory, where samples were held at −20 °C until processing.

Amphibians

Four species of aquatic-stage salamanders and two species of adult frogs were captured at 50 and 21 sites, respectively, where dragonfly larvae were also sampled (Table S2). Amphibians were captured using dip nets or minnow traps, and then, nonlethal tissue samples (tail clips for salamanders and toe clips for frogs) were collected from each individual[21] and stored frozen at −20 °C until Hg determination.

Laboratory Preparation and Analysis

In the laboratory, all dragonfly, fish, and amphibian samples were dried to a constant mass and either homogenized to a fine powder or prepared for analysis whole. We determined THg concentrations on all samples, following Environmental Protection Agency (EPA) method 7473,[22] and MeHg concentrations on a subset (n = 652) of samples, following EPA method 1630.[23] Methodological details and quality assurance/quality control details can be found in Supporting Information. Total Hg analysis was conducted at the United States Geological Survey (USGS) contaminant ecology research lab (87% of samples; N = 12845), the University of Maine (12% of samples; N = 1842), and Dartmouth College (1% of samples; N = 144). Interlaboratory comparisons were conducted between the two primary laboratories to ensure data comparability, which are summarized in Supporting Information.

Statistical Analyses

All data were natural-log transformed prior to analyses to normalize residuals and meet the assumptions for parametric statistical tests. Unless specified otherwise, we present estimates of the central tendency as geometric or back-transformed least-squares means because of the log-normal distribution of concentration data. All statistical analyses were conducted using JMP V14.0 (SAS Institute, Cary, NC).

Intrinsic and Extrinsic Factors

We examined the importance of key intrinsic (taxonomy and body length) and extrinsic (habitat type and ecoregion[24]) variables on dragonfly THg concentrations using a linear mixed-effects model that included family, habitat type, and level 1 ecoregions as fixed effects and site and categorical sampling year as random effects. We also included total body length (TL) as a covariate, and TL × family, TL × site, and TL × family × site interactions. The initial model results did not support including TL in the model (see Results and Discussion); therefore, we removed TL and its associated interactions from the model. The variables selected for the model are those that have been shown to be associated with biological Hg concentrations at broad scales[7] and represent the underlying processes and mechanisms that drive Hg cycling and bioaccumulation. Although sampling spanned several months, each site was generally only sampled once per year, and sites with similar climates and hydrology were usually sampled within a similar time frame. Therefore, we did not explicitly assess or account for seasonal variation in our models.

Taxonomic Conversions

We next examined relationships in THg concentrations among dragonfly families by correlating paired geometric mean THg concentrations of different families collected concurrently at the same site. For each family pair, we developed linear regression equations for predicting THg concentrations in one family based on those in another. The Aeshnidae family represented the largest proportion (41%) of the data set and was consistently correlated with each of the other families. Therefore, we used the linear regression equations to convert each individual sample from the other families to an Aeshnid-equivalent THg concentration. By normalizing concentrations to a consistent family, the Aeshnid-equivalent concentrations allow for more-representative comparisons across sites and time periods without confounding effects of taxonomy.

Bioindicator Suitability

The suitability of dragonfly larvae as biosentinels of THg exposure in aquatic ecosystems is partly predicated on their ability to provide inference about THg exposure in other aquatic organisms at multiple levels of the food web. To examine this, we paired collections of dragonflies at a subset of sites where either fish or amphibian species were also sampled. At each of these sites, we generated an Aeshnid-equivalent site–year geometric mean THg concentration to ensure representative units were used for dragonfly Hg concentrations.

We similarly calculated site–year geometric mean THg concentrations for each fish guild and both salamanders and frogs. However, limited sample sizes and a diverse range of species precluded adjusting concentrations to a single taxon for each group, as was done with the dragonfly larvae. We chose a guild-based approach because it aggregated related fish species with a similar trophic ecology into discrete categories that are applicable across a wide geographic range and multiple habitat types. After calculating site- and taxa-specific geometric mean THg concentrations for the fish and amphibians, we used linear regression models for each taxonomic group to quantify the relationship in THg concentrations between dragonfly larvae and other aquatic taxa.

Results and Discussion

From 2009 to 2018, more than 4000 citizen scientists contributed to the collection of 14831 dragonfly larvae from 457 individual sites in 100 different NPS units across the US, including Alaska and Hawaii (Figure ). The geometric mean (±geometric standard error) THg concentration [ng/g dry weight (dw)] of all individual dragonfly larvae was 125 ± 2.2 (arithmetic mean ± standard deviation = 181 ± 217). The park units (park unit type abbreviations are defined in Table S3) with the highest geometric mean THg concentrations across all sites included Colorado NM, Bear Meadows NNL, Minute Man NHP, and Maurice WSR, whereas those with the lowest concentrations included Bear Creek Lake Park, Padre Island NS, Fort Worth NNL, and Buffalo NR (Figure ). Individual sites with the highest THg concentrations were in Olympic NP (1411 ± 41), Acadia NP (1319 ± 319), Glen Canyon NRA (1207 ± 86), Maurice WSR (905 ± 38; 854 ± 32; 745 ± 36), Capitol Reef NP (672 ± 92), and Yellowstone NP (626 ± 20), whereas those with the lowest THg concentrations were found in Bear Creek Lake Park (10.4 ± 1.2), Gulf Island NS (16.4 ± 1.6), Padre Island NS (17.6 ± 5.2), Buffalo NR (20.6 ± 2.5), and Forth Worth NNL (21.7 ± 1.8).

Figure 1

Total mercury (THg) concentrations (ng/g dry weight) in dragonfly larvae collected from National Park units across the conterminous US, Alaska, and Hawaii. Bars represent park unit geometric means (±SE), and circles represent geometric mean THg concentrations for individual sampling sites (populations) within each unit. Colors represent aggregations of US Department of Interior regions representative of western (orange), central (blue), and eastern (green) areas, as illustrated in the inset map. Arrows and associated text represent sites with geometric mean THg concentrations beyond the extent of the x-axis. Unit name suffixes are provided in Table S3. Note that the map geography is not to scale. Alterations were made to include Alaska and Hawaii within the map frame. Bear Creek Lake Park is not a NPS administered unit.

In a subset of individuals (N = 652) analyzed for both MeHg and THg, MeHg concentrations ranged from 12.3 to 1870 ng/g dw and were strongly correlated with THg concentrations (R2 = 0.96, p < 0.0001; Figure S2). Additionally, the % MeHg (proportion of THg in the MeHg form) averaged (±SE) 79.9 ± 0.5% of THg concentrations, which is similar to the ranges (83–94%) reported elsewhere.[16,18] Collectively, these findings suggest that, unlike other aquatic macroinvertebrate species where low and variable % MeHg is a major limitation to their use as Hg bioindicators, THg serves as an effective proxy for MeHg concentrations in dragonfly larvae. Therefore, we use THg concentrations for further data analysis.

Variation within and among Parks

Individual dragonfly larvae THg concentrations spanned four orders of magnitude, from 1.0 to 3795 ng/g dw across all sites. Consistent with such a broad range in concentrations, the data set’s geometric coefficient of variation (GCV) was 96%, indicating a high degree of variability across the landscape. However, the average GCV (±SE) within sites (with n ≥ 10) was nearly 3-fold less (35.3 ± 0.95%), and within-site GCVs were not correlated with their paired site-specific geometric mean THg concentrations (p = 0.84, n = 412), indicating that the within-site variation in THg concentrations was relatively similar regardless of the overall Hg condition of a site. These findings suggest that the variation in THg concentrations among sites at a continental scale is substantially higher than within sites. Indeed, site-specific geometric mean THg concentrations ranged 135-fold (10.4–1411 ng/g dw; Figure ), and the GCV among all sites (86%) was more than double the average GCV within sites, where concentrations ranged 5.7-fold, on average (range = 1.2–100-fold). Importantly, even within individual national park units, which is a substantially finer scale than across the US, and where sites are in relatively close proximity to one another, the average (±SE) for among-site GCV (62.7 ± 6.5%) was similar to the GCV among national parks (68.7%). This highlights the importance of site-specific characteristics on Hg bioaccumulation.

Site is often a dominant factor in influencing Hg bioaccumulation,[7,25,26] resulting in substantial geographic variation. Mean fish THg concentrations in 206 species varied 496-fold across >4000 locations throughout western North America,[7] and species-normalized concentrations varied 96-fold in nearly 2000 sites across Canada,[27] similar to the range of variation we found (135-fold) in dragonfly larvae. This variation is not constrained to continental-scale assessments nor likely due to the direct inorganic Hg contamination. For example, fish THg concentrations in 28 high-elevation lakes of eastern Oregon, US varied 18-fold,[28] despite being free from direct watershed disturbance and receiving a similar atmospheric Hg deposition. Furthermore, THg in fish from national parks in the western US varied between 10- and 24-fold among sites.[29,30] Even at more local scales, fish Hg concentrations of individual wetlands in the San Francisco Bay region of California varied by 15-fold despite only being separated by narrow levees.[31] Likewise, we found high variation within individual national parks, where dragonfly larvae THg concentrations varied 33-fold across 13 sites in Olympic NP, 32-fold between two sites in Gulf Islands NS, 13-fold in 85 sites from Mount Rainier NP, 12-fold among 20 sites in Acadia NP, and 11-fold across 17 sites in Yellowstone NP.

Because inorganic Hg contamination is generally decoupled from MeHg production, uptake, and bioaccumulation,[4] the high variability among sites (both within and among national parks) is likely a result of site-specific biogeochemical and ecological factors that control net MeHg production and bioavailability within ecosystems. In fact, the variation in dragonfly THg concentrations (135-fold across sites) was an order of magnitude higher than the variation in THg wet deposition across the US (∼10-fold across sites; http://nadp.slh.wisc.edu/). Thus, we next examined key extrinsic factors known to influence biological MeHg exposure along with intrinsic factors that influence individual-level bioaccumulation.

In our initial global model, the length × family (p < 0.0001) and length × site (p < 0.0001) interactions indicated potential site- or family-dependent relationships between length and THg concentrations. Therefore, we conducted individual regressions for each family and sampling site where there were at least 15 samples from a given year (site–years across five families; Table S4). Length and THg concentrations were significantly correlated (p < 0.05) in only 82 (out of 298; 28%) site–years (Aeshnidae = 34%, Cordulegastridae = 7%, Gomphidae = 17%, and Libellulidae/Corduliidae = 26%). However, of the 82 site–years with a significant correlation between THg concentration and dragonfly length, 54 (66%) were positive correlations, and the remaining 28 (34%) were negative. Moreover, the average R2 value for all significant regressions (positive or negative) was only 0.37, and only 22% of site–years had R2 values >0.50, indicating that the body length was neither an overwhelming nor consistent variable associated with dragonfly THg concentrations. Because of the limited explanatory power and substantial variation in the directional effect of body length on THg concentrations, we did not size-adjust dragonfly THg concentrations, and we removed the length from subsequent models.

Extrinsic Drivers: Habitat and Ecoregion

In the reduced model (i.e., excluding length), dragonfly THg concentrations differed among primary habitats (F4,438.3 = 7.30; p < 0.0001), secondary habitats (F1,1488.4 = 3.54; p < 0.0001), and among ecoregions (F1,2435.4 = 4.99, p < 0.0001). Among the five primary habitat types, dragonfly larvae THg concentrations were highest in rivers and streams (Figure ). Post hoc pairwise contrasts indicated that riverine dragonflies had higher concentrations than those in all other habitat types except streams (streams, F1,432 = 1.17, p = 0.28; lakes, F1,429.8 = 11.3, p = 0.0009; wetlands, F1,432.1 = 9.03, p = 0.003; ponds, F1,435.2 = 15.1, p < 0.0001); stream dragonflies also were higher than those from the other three habitats (wetlands, F1,431.9 = 8.35, p = 0.004; ponds, F1,439.2=18.0, p < 0.0001; lakes, F1,429 = 9.75, p = 0.002), which did not differ from one another.

Figure 2

Least squares mean THg concentrations (ng/g dry weight) in dragonfly larvae among habitat types in national park units across the US. Least squares mean THg concentrations account for the effects of family, ecoregion, site, and year. Patterned bars in secondary habitats represent those bounded by extensive marginal wetlands. Letters represent the significance at α = 0.05 based on pairwise contrasts and test slices.

Among nested secondary habitats (Figure ), THg concentrations differed within rivers (F2,436.1 = 84.1, p < 0.0001), wetlands (F2,427.7 = 7.74, p = 0.0005), and streams (F2.430.8 = 7.69, p = 0.0005) but not lakes (F2,424.6 = 0.29, p = 0.74) nor ponds (F2,857.2 = 1.42, p = 0.24). In rivers and streams, THg concentrations were 71% and 76% higher, respectively, in sites with adjacent riparian wetlands and floodplains (rivers, F1,431.8 = 4.43, p = 0.03; streams, F1,431.4 = 14.9, p < 0.0001). However, perennial and intermittent streams without marginal wetlands did not differ from one another (F1,430.5 = 0.85, p = 0.36), suggesting that the proximity to wetlands may have a stronger influence on MeHg bioaccumulation in streams than water permanence. In contrast, within wetlands themselves, seasonally inundated emergent and forested wetlands had 58% and 112% higher THg concentrations, respectively, than permanently inundated wetlands (emergent, F1,427.9 = 6.81, p = 0.009; forested, F1,425.5 = 6.81, p < 0.0001). Across all of the secondary habitats, dragonflies from perennial rivers with marginal wetlands (the habitat with the highest concentrations) were 3.1-fold higher than those from the lowest habitat (permanent emergent wetlands). Although differences were not consistently significant in all habitat types, those bounded by marginal wetlands averaged THg concentrations that were 35% higher than similar habitats without wetlands.

These habitat differences further highlight the importance of ecological factors as key drivers of spatial variability in biological Hg concentrations. Hg exposure being the highest in lotic and lowest in lentic habitats is consistent with other independent data sets, where fish THg concentrations were 21%–61% higher in lotic than lentic waterbodies across western North America[7] and northeastern US lakes.[32] An important caveat to this is that dragonfly larvae THg concentrations from lakes (particularly large lakes) primarily represent a littoral zone analysis and likely do not effectively integrate pelagic pathways of Hg bioaccumulation. Importantly, water bodies bounded by extensive wetlands and floodplains generally had the highest dragonfly THg concentrations, despite lower concentrations in dragonflies sampled from wetlands themselves (except seasonal forested wetlands) than in either rivers or streams. Wetlands are known sites of MeHg production,[33] and wetland density can be a strong predictor of aqueous MeHg concentrations of streams.[34,35] The apparent discrepancy between these habitats may be a function of dissolved organic carbon (DOC) dynamics. Elevated DOC is common in wetlands, and DOC can both promote MeHg production and transport but also inhibit MeHg bioaccumulation.[36−38] Thus, high DOC in wetlands can result in lower MeHg bioaccumulation factors within wetlands themselves by binding MeHg and making it less bioavailable.[39] At the same time, DOC can facilitate THg and MeHg transport to surrounding rivers and streams where it may become more bioavailable in a lower DOC environment. Indeed, dragonfly THg concentrations were positively correlated with lake DOC concentrations across a series of lakes in the northeastern US, but the DOC inhibition effect was observed in bioaccumulation factors for dragonfly THg across the range of measured DOC.[18]

At a broader spatial scale, the influence of the habitat type and ecological factors on Hg bioaccumulation is illustrated by differences among level 1 ecoregions, which represent broad areas of similar geography, geology, climate, and basic biological organization, which in turn influence the structure and function of specific ecosystems.[40] After statistically accounting for taxonomy, habitat, and site, we found substantial variation in THg concentrations among these ecoregions (Figure S3), with a 3.6-fold range in THg concentrations between the lowest (Great Plains) and the highest (North American Deserts) ecoregions. Arid ecoregions such as North American Deserts and Southern Semiarid Highlands were particularly elevated in THg concentrations, both in the studied dragonflies and in previously reported fishes across western North America,[7] whereas THg concentrations in Taiga were similarly low in both dragonflies and fish. Hydrology can be an important driver of MeHg production and bioaccumulation,[41] and arid ecoregions represent areas where many waterbodies are only seasonally inundated and anthropogenic water impoundments (i.e., reservoirs) are widespread. Drying and rewetting of littoral sediments inherent to reservoirs can exacerbate Hg methylation and increase MeHg bioaccumulation in food webs both within and downstream of impoundments,[42−44] which may influence dragonfly THg concentrations in these arid areas. Similarly, we found that dragonflies from seasonal wetlands had substantially higher THg concentrations than those from permanently inundated wetlands, findings consistent with those of others from managed wetland systems where cyclic wetting and drying patterns allowed for the reoxidation of terminal electron acceptors that are important for MeHg production.[9,45,46]

Intrinsic Drivers: Taxonomy and Body Size

Taxonomy can influence Hg concentrations across phylogenetic categories because it incorporates differences in foraging ecology, habitat use, and bioenergetics. Consistent with this, we found that THg concentrations (ng/g dw) differed among larval dragonfly families (F4,14585 = 251.0, p < 0.0001; Figure S4). Aeshnidae (141 ± 11.0) had the highest THg concentrations, followed by both Macromiidae (131 ± 10.9) and Cordulegastridae (131 ± 10.7), which did not differ from one another. Libellulidae and Corduliidae (116 ± 9.0), which were categorized together because of the uncertainties in their differentiation at this resolution, were significantly lower than the preceding three families but higher than Gomphidae (95.0 ± 7.5), which had the lowest THg concentrations, and were 1.5 times lower on average than Aeshnids (Figure S4).

The mechanisms for taxonomic differences require more study but are likely tied to some combination of foraging ecology, habitat use, and physiology.[47] Dragonfly larvae are commonly described in terms of one of four categories: claspers (Aeshnidae), sprawlers (Macromiidae, Corduliidae, Libellulidae), hiders (Cordulegastridae), and burrowers (Gomphidae), which represent their preferred microhabitats and activity levels.[15] Although some species in each family fall into different categories, these generalized microhabitat distinctions likely influence both diet and energetics. We limited our identification to the family level, which may add unexplained variation if species–specific differences in THg bioaccumulation are pronounced within families. This requires further investigation, but there was low variation among Gomphidae species across 17 lakes in the Laurentian Great Lakes region.[17]

Resilience to the confounding effects of intrinsic factors that influence Hg concentrations is among the most important aspects of effective Hg biosentinels. Of the two common intrinsic factors that we examined (body size[48] and taxonomy[7]), which integrate other important factors such as foraging ecology and physiology, only taxonomy consistently influenced THg concentrations. However, paired geometric mean THg concentrations of dragonfly larvae sampled from the same location at the same time were strongly correlated among all family pairs (Figure ; equations S1–S5), except Macromiidae and Cordulegastridae, a pair we could not test because they were concurrently sampled at only five sites. These relationships between families at a national scale facilitate a simple approach to reliably convert THg concentrations from one family to those of another family, making it possible for equivalent comparisons among locations or over time where the same families were not sampled. Similar approaches have been used for both wildlife and fish[27,49] to ensure robust comparisons among locations or to estimate the likely wildlife Hg exposure on the basis of their prey.[50,51] To facilitate the use of dragonfly THg data across landscapes where the same families may not be sampled, we provide regressions for each family pair in Supporting Information.

Figure 3

Relationships of total mercury (THg) concentrations in dragonfly larvae among families. Each data point represents a paired site–year geometric mean THg concentration for each family. The x-axis represents geometric mean concentrations for the family listed in each pane; y-axis represents paired THg concentrations in the Aeshnidae family. Libellulidae and Corduliidae are combined as a single family group because their distinction can be uncertain without the identification of the genus.

Relationships with Fish and Wildlife Hg Exposure

Aeshnid-equivalent THg concentrations were positively correlated with THg concentrations in all four fish guilds (Figure ): piscivores (F1,21 = 20.9; p = 0.0002, N = 22 sites), sunfish (F1,12 = 10.4; p = 0.008, N = 13 sites), salmonids (F1,9 = 15.5; p = 0.004, N = 10 sites), and forage fish (F1,12 = 65.1; p < 0.0001, N = 13 sites), as well as with those in both salamanders (F1,50 = 45.1; p < 0.0001, N = 51 sites) and frogs (F1,21 = 12.6; p = 0.002, N = 22 sites) (Figure ).

Figure 4

Relationship between Aeshnid-equivalent dragonfly larvae THg concentrations (ng/g dw) and THg concentrations (μg/g wet weight (ww)) in four fish guilds (left panel), salamander tail clips (ng/g dw; top right panel), and frog toe clips (ng/g dw; bottom right panel. Each data point represents the paired geometric mean THg concentrations for sites where both dragonfly larvae and either fish or amphibians were sampled. Fish THg concentrations are reported in wet weight to facilitate comparisons to common health benchmarks, whereas dragonfly larvae and amphibian concentrations are reported on a dry weight basis to reduce variance associated with external moisture.

The ability to predict THg concentrations in other components of aquatic food webs is a particularly important aspect of dragonflies as Hg biosentinels. These relationships, spanning classes of organisms as well as guilds within them, at a national scale, suggests transferability to broad scales and across ecosystems. These findings provide the underpinnings of a tool for estimating the potential risk due to Hg exposure across multiple taxa. Haro et al.[17] similarly showed strong relationships between site-specific THg concentrations in Gomphidae and those in both forage fish and predatory fish. However, they estimated that Gomphidae THg concentrations of 40 ng/g of dw were associated with muscle THg concentrations in predatory fishes equivalent to the US EPA MeHg criterion value for the protection of human health (0.30 μg/g of wet weight (ww)). Our findings for piscivorous fishes were substantially more conservative, suggesting that the EPA criterion was reached when Aeshnid-equivalent THg concentrations exceeded 162 ng/g of dw, which corresponds to approximately 112 ng/g in Gomphidae. The reason for this discrepancy is unclear but could be associated with differences in the geographic scale of these studies (i.e., models derived from similar habitats and a narrower geographic region (Great Lakes region) as opposed to those from a more diverse assortment of habitats at a continental scale). Alternatively, it could simply be the result of different statistical approaches used to generate site-specific mean THg concentrations (geometric means vs arithmetic means). Regardless, the similar strength and slopes of the relationships between the two studies highlight the utility of using dragonflies as an index of fish Hg exposure across a variety of scales.

Integrated Impairment Index

Dietary- and tissue-based impairment benchmarks provide estimates of the potential toxicological risk of Hg to fish and wildlife[52−55] as well as to humans through fish consumption.[56] For each fish guild, we used the linear regressions described above to model the Aeshnid-equivalent THg concentrations that would be indicative of levels exceeding published health benchmarks that span a range of potential hazard (Table S5). We also modeled the Aeshnid-equivalent THg concentrations that would correspond with each fish guild (except forage fish) exceeding US EPA methylmercury criterion (0.30 μg/g of ww) for the protection of human health.[56]

Using these benchmarks, in concert with the corresponding dragonfly THg concentrations for different fish guilds, we classified Aeshnid-equivalent THg concentrations into five indices of progressively increasing severity that incorporated both the magnitude of potential impairment as well as the number of guilds exceeding these individual benchmarks (Figure , Table S5). As opposed to simply evaluating individual species, this approach considers the potential risk across multiple guilds, allowing for a broader assessment of the potential ecological and human health exposure risks in freshwater ecosystems. The established impairment benchmarks are still lacking for other taxa, such as amphibians, but their correlations with dragonfly larvae suggest that the development of benchmarks could lead to further refinements of integrated impairment assessments.

Figure 5

Integrated impairment indices for potential ecosystem risk to mercury, in Aeshnid-equivalent units (ng/g dry weight total mercury (THg)). Integrated impairment index categories (specified by different colored boxes) represent concentrations corresponding to exceedances of a range of individual published toxicity benchmarks[53−55] across several fish guilds (piscivores, sunfish, salmonids, and forage fish). Within each of these guilds, the potential impairment was assessed based on (1) fish diet (i.e., risk to other fish as prey), (2) fish tissue (labeled as fish health; i.e., risk to fish themselves based on their tissue concentrations), and (3) avian diet (i.e., risk to fish-eating birds as prey), with three possible levels of severity (low, moderate, and high as white, grey, and black arrows, respectively) for each of these benchmark types. Potential human risk is incorporated through the US EPA methylmercury criterion for the protection of human health, which was converted to a whole-body equivalent.[7] The integrated impairment indices are defined by the number and severity of benchmarks exceeded as well as the number of guilds exceeding them. The symbols and lines represent the modeled Aeshnid-equivalent THg concentrations associated with the THg concentrations of each impairment benchmark for each fish guild. Circles represent fish dietary benchmarks, triangles represent fish health benchmarks, diamonds represent avian dietary benchmarks, and squares represent US EPA MeHg criterion. White, gray, and black symbols reflect low, moderate, and high-severity benchmarks, respectively.

This integrated approach allows resource managers to simultaneously evaluate the potential exposure and hazard across a range of ecological receptors to more thoroughly assess the threat Hg may pose to ecosystem function. Overall, 10% of the 877 site–years had geometric mean Aeshnid-equivalent THg concentrations that were below any of the deleterious effect benchmarks, and 22% were associated with the low-hazard index. In contrast, more than half of the site–years (56%) were categorized as moderate hazards. Finally, 11 and 1% of sites had geometric mean Aeshnid-equivalent THg concentrations that were high enough for them to be classified in the high- or severe-hazard indices, respectively. Geographically, the distribution of risk categories was heterogeneous across the continent (Figure ). Moreover, within individual park units there were often sites that spanned the entire range of benchmark categories (Figure ).

Figure 6

Integrated impairment indices for all 457 sites sampled between 2009 and 2018. Integrated impairment indices are derived from Aeshnid-equivalent geometric mean THg concentrations for each site–year and their corresponding association with Hg exposure in fish and other wildlife. The inset for the cumulative frequency distribution illustrates the proportion of sites and years that fall into each of the five categories. Note that the map geography was altered to include Alaska and Hawaii within the map frame.

The national-scale citizen-science network leveraged herein represents the most expansive assessment of the MeHg distribution and the potential impairment risk using aquatic macroinvertebrates as biosentinels. We show that dragonfly larvae provide an effective model for understanding important factors of Hg bioaccumulation and can be used to estimate potential ecosystem impairment. We also demonstrate that widespread Hg contamination in protected NPS lands poses risks to trust resources within them. Utilizing a citizen-science framework also provides a cost-effective approach for biological sampling while also engaging the public in experiential learning opportunities. Future applications of dragonfly larvae Hg monitoring could be particularly insightful for better understanding the landscape, biogeochemical, and climatological drivers of Hg availability to aquatic food webs; predicting the potential risk to vertebrates; and evaluating the effectiveness of Hg reductions, both regionally and globally.