# Evidence that Pacific tuna mercury levels are driven by marine methylmercury production and anthropogenic inputs

## Significance

Humans are exposed to toxic methylmercury mainly by consuming marine fish. New environmental policies under the Minamata Convention rely on a yet-poorly-known understanding of how mercury emissions translate into fish methylmercury levels. Here, we provide the first detailed map of mercury concentrations from skipjack tuna across the Pacific. Our study shows that the natural functioning of the global ocean has an important influence on tuna mercury concentrations, specifically in relation to the depth at which methylmercury concentrations peak in the water column. However, mercury inputs originating from anthropogenic sources are also detectable, leading to enhanced tuna mercury levels in the northwestern Pacific Ocean that cannot be explained solely by oceanic processes.

## Abstract

Pacific Ocean tuna is among the most-consumed seafood products but contains relatively high levels of the neurotoxin methylmercury. Limited observations suggest tuna mercury levels vary in space and time, yet the drivers are not well understood. Here, we map mercury concentrations in skipjack tuna across the Pacific Ocean and build generalized additive models to quantify the anthropogenic, ecological, and biogeochemical drivers. Skipjack mercury levels display a fivefold spatial gradient, with maximum concentrations in the northwest near Asia, intermediate values in the east, and the lowest levels in the west, southwest, and central Pacific. Large spatial differences can be explained by the depth of the seawater methylmercury peak near low-oxygen zones, leading to enhanced tuna mercury concentrations in regions where oxygen depletion is shallow. Despite this natural biogeochemical control, the mercury hotspot in tuna caught near Asia is explained by elevated atmospheric mercury concentrations and/or mercury river inputs to the coastal shelf. While we cannot ignore the legacy mercury contribution from other regions to the Pacific Ocean (e.g., North America and Europe), our results suggest that recent anthropogenic mercury release, which is currently largest in Asia, contributes directly to present-day human mercury exposure.

Mercury (Hg) is a global pollutant that affects both ecosystems and human health. Mercury is emitted to the atmosphere by natural processes (e.g., volcanism) and anthropogenic activities (e.g., fossil fuel combustion and artisanal and small-scale gold mining) (1) mainly as gaseous elemental Hg0. Atmospheric Hg can be deposited or taken up at the ocean surface (2). Rivers deliver significant amounts of Hg to estuarine and coastal oceans, but only a small fraction of this Hg is transported to open ocean regions (3). In the ocean, a substantial fraction of inorganic Hg is naturally transformed to methylmercury (MeHg) in the subsurface waters (4, 5). The neurotoxic MeHg naturally bioaccumulates and biomagnifies in aquatic food webs, leading to elevated concentrations in higher-trophic-level organisms such as tunas. Humans are exposed to MeHg mainly by the consumption of marine fishes, and tunas in particular, with MeHg representing the major chemical form (>91%) of total Hg in tuna muscle (68). The United Nations Environmental Program Minamata Convention (131 parties in June 2021) aims to protect human health and the environment from anthropogenic Hg. Governments and policy makers must therefore balance economic, environmental, and health interests. Understanding the factors that drive MeHg concentrations in tunas, and how climate change interferes in the process is therefore critical to policy makers and, ultimately, to limiting human exposure to MeHg.

In line with the spatial distribution of Hg emissions, observations show higher atmospheric Hg0 levels in the northern hemisphere compared to the southern hemisphere (9). Sediment and peat records of atmospheric Hg deposition also suggest a larger anthropogenic Hg enrichment in the northern hemisphere (10). Anthropogenic Hg emissions are thought to have tripled the concentrations of total Hg in the thermocline waters (100 to 1,000 m) of the global ocean relative to deeper older waters (11), yet still unclear is how these anthropogenic Hg inputs are converted into methylmercury in oceans. Subsurface MeHg concentrations are higher in the eastern Pacific Ocean, with coastal and upwelling areas exhibiting the highest proportion of MeHg relative to total Hg (4, 12).

Methylmercury uptake and biomagnification associated with fish age (and length), feeding ecology, and habitat use are known factors influencing Hg bioaccumulation in tunas (8, 1315). A recent global-scale study documenting the distribution of tuna Hg concentrations at 11 different locations found up to an eightfold difference among sites but did not identify the drivers of this variability (16). Another study conducted at high spatial resolution in the southwestern Pacific Ocean suggested that ocean physics, local ocean biogeochemistry, and tuna habitat play important roles in determining tuna Hg concentrations (8, 17). Meanwhile, temporal trends in tuna Hg show contrasting results between oceanic basins. Decreasing Hg concentrations were observed in Atlantic bluefin tuna (Thunnus thynnus) from the northern Atlantic Ocean between 1990 and 2000 (18). In the Pacific Ocean, stable Hg concentrations between 2001 and 2018 were observed in tropical tunas (i.e., yellowfin, Thunnus albacares; bigeye, Thunnus obesus; and skipjack, Katsuwonus pelamis) from the southwestern Pacific (15), while increasing Hg levels in yellowfin and bigeye tunas since 1970 were documented in the central north Pacific (CNPO) (19). These contrasting trends suggest that tuna Hg may be sensitive to spatiotemporal changes in anthropogenic Hg emissions and ocean uptake and/or regional variations in biogeochemical processes or foraging ecology that affect tuna exposure to MeHg. Studies of large spatial coverage are therefore needed to investigate the control and relative importance of atmospheric Hg sources, oceanic Hg methylation, and foraging ecology on tuna Hg concentrations.

Here, we report the highest spatially resolved dataset of Hg concentrations of skipjack tuna, with records from six distinct areas of the Pacific Ocean (n = 661, from ∼45°N to 25°S, Fig. 1A). As the third most consumed marine fish species worldwide (20), it is a species of high commercial and human health importance. Although classified as highly migratory, this epipelagic and fast-growing tuna species shows relatively restricted movements and high site fidelity rates in the Pacific Ocean (21, 22) and thus should be able to reflect spatial MeHg patterns in surface water food webs. We applied generalized additive models (GAM) to 1) explore geographical patterns, including hemispheric differences, of tuna Hg concentrations and 2) disentangle the driving factors of tuna Hg spatial variability. In particular, we used atmospheric Hg0 model estimates, ecological tracers in tuna, biogeochemical estimates at the base of food webs, and oceanographic variables to explore the relative contribution of anthropogenic Hg emissions, ocean biogeochemistry, and tuna foraging ecology to observed tuna Hg concentrations.

Fig. 1.

Skipjack tuna sample provenance and mercury concentrations. (A) Map of skipjack sampling locations, with the size of the circles proportional to the number of individuals collected. (B) Power-law relationship between log(observed Hg) and skipjack fork length. The relationship was fitted on samples from the SWPO and WCPO only, but Hg residuals were calculated for all samples (symbolized by the black arrows). Oceanic areas: NWPO, CNPO, NEPO, EPO, SWPO, and WCPO.

## Results and Discussion

### Spatial Variability of Tuna Mercury Concentrations.

As Hg is known to bioaccumulate with fish size (23), we fitted a length-based model to the observed Hg concentrations to calculate length-standardized Hg concentrations (Fig. 1B, at mean length = 60 cm fork length; see Materials and Methods), hereafter defined as standardized Hg concentrations, in contrast to observed Hg concentrations. This normalizing approach allowed us to explore the spatial distribution of Hg concentration anomalies not explained by Hg bioaccumulation and fish size differences among individuals but likely due instead to other processes (e.g., anthropogenic emissions, ocean biogeochemistry, or tuna foraging ecology) (8, 15).

Observed Hg concentrations in skipjack ranged from 0.01 to 3.15 µg ⋅ g−1 dry weight (0.56 ± 0.51 µg ⋅ g−1 dw, mean ± SD; Fig. 2A and SI Appendix, Fig. S1), which is comparable to values reported near Hawaii and in the western Indian Ocean for the same species (2426). The relatively low Hg levels in skipjack tuna are consistent with the ecology and biology of this species. Compared to other tropical tuna species that exhibit higher Hg concentrations, skipjack tuna are characterized by a lower trophic position, a shallower vertical habitat that allows access to epipelagic prey with lower Hg concentrations, and a shorter lifespan (14, 27).

Fig. 2.

Spatial variability of skipjack mercury concentrations. Smoothed spatial contour maps of (A) observed and (B) standardized Hg concentrations (micrograms ⋅ grams−1, dw) in skipjack white muscle samples from the Pacific Ocean. The black dots represent the location of skipjack samples. Ocean areas correspond to the sample origin: NWPO, CNPO, NEPO, EPO, SWPO, and WCPO. The transparent dots represent the location of seawater samples with available and published MeHg data (see Materials and Methods).

Observed and standardized Hg concentrations in skipjack tuna were both highly variable in the Pacific Ocean. Standardized Hg concentrations were significantly higher in the northern than in the southern hemisphere (Fig. 2B), which is consistent with latitudinal Hg gradients observed in other Pacific Ocean reservoirs (i.e., atmosphere, seawater, and marine sediments) (Fig. 3 and SI Appendix, Fig. S2). Significant spatial differences of standardized Hg concentrations were also detected at smaller spatial scales within the Pacific, among the defined oceanic areas (Kruskal–Wallis test, P < 0.001; Fig. 2B and SI Appendix, Fig. S1). The highest standardized Hg concentrations were found in the northwestern Pacific (NWPO), where they were 1.5 to 2 times higher than in the CNPO and eastern regions (i.e., northeastern [NEPO] and eastern [EPO] Pacific) and 4 to 5 times higher than in the southwestern (SWPO) and the western central (WCPO) Pacific Ocean (Fig. 2B). The only two available spatial studies of Hg in tropical yellowfin and bigeye tunas also suggested higher Hg concentrations in the EPO compared to the central and western regions (8, 28).

Fig. 3.

Hemispheric mercury gradients in different Pacific Ocean reservoirs. The boxplots illustrate the hemispheric gradient of (A) atmospheric Hg0 model estimates (nanograms ⋅ meters−3) extracted at tuna sampling locations (see Materials and Methods), (B) observed total Hg concentrations (picomolar) in seawater (12, 34, 38, 53, 54, 56), (C) observed MeHg concentrations at peak (picomolar) in seawater (12, 34, 5357), (D) standardized total Hg concentrations (micrograms ⋅ grams−1, dw) in skipjack tuna (this study), and (E) observed total Hg concentrations (nanograms ⋅ grams−1, dw) in marine sediments (2, 7378). ** indicates significant differences between the northern and southern hemispheres (Kruskal–Wallis test; P < 0.01).

The main drivers of this spatial variability were identified with GAM. The optimal model explained 72% (deviance explained [DE]) of the spatial variability of observed tuna Hg concentrations and included fish length, proxies for the depth of seawater MeHg peak production (i.e., the depth of net heterotrophy and [O2] in subsurface waters; see Materials and Methods) and atmospheric Hg0 model estimates (Fig. 4). Baseline-corrected tuna δ15N values were used as a proxy for the tuna trophic position to investigate trophic processes and Hg biomagnification through the pelagic food webs. They were, however, not selected in the best model, indicating that spatial trophic effects (i.e., geographical changes in foraging ecology) have a limited influence on the spatial variability of tuna Hg concentrations, as already observed in the WCPO (8). Furthermore, the absence of the effect of particulate organic matter δ15N estimates in our optimal GAM also indicates that spatial changes of the nitrogen cycle at the base of the marine food web may be of weak importance to explain spatial changes of tuna Hg concentrations at the large Pacific Ocean scale. Conversely, given the selected drivers in the optimal GAM, the spatial variability of skipjack Hg concentrations is more likely to result from 1) variable seawater MeHg concentrations at the base of the food web among oceanic areas and/or 2) variable amounts of surface ocean inorganic Hg loading from the atmosphere, as discussed in the two following sections.

Fig. 4.

Optimal GAM predicting observed mercury concentrations in skipjack tuna. The main drivers of the spatial variability of log(observed Hg) were (A) fish length, (B) the depth of net heterotrophy, (C) [O2] in subsurface waters, and (D) Hg0 model estimates. The DE percentage of the model is reported on the top of the figure. Colored dots are the partial residuals of the GAM by oceanic areas. Black lines show the expected values, and gray bands show CIs for the expected value, which are twice the SE. Oceanic areas: NWPO, CNPO, NEPO, EPO, SWPO, and WCPO.

### Natural Marine Biogeochemical Drivers of the Spatial Variability of Skipjack Mercury Concentrations.

A key aspect when investigating Hg variability in marine biota is its relationship with seawater MeHg concentrations, which may vary as a function of biogeochemical processes. In the Pacific Ocean, observed MeHg data in seawater are scarce and are available at limited and low spatial resolution in comparison with our skipjack Hg dataset (Fig. 2). To overcome this limitation, we used a combination of biogeochemical and physical proxies of seawater peak MeHg concentrations to assess the influence of ocean biogeochemistry on skipjack Hg burden. In the optimal model, the spatial variability of skipjack Hg concentrations was largely explained by the depth of net heterotrophy (i.e., the depth at which autotrophy is overcome by heterotrophy), which is a proxy for the depth of peak seawater MeHg concentrations, and [O2] in subsurface waters (Fig. 4 B and C). Response curves indicated that these two variables combined might be of particular importance to explain the west-to-east gradient of tuna Hg content. Higher tuna Hg concentrations found in the east (i.e., NEPO and EPO) were linked to low subsurface [O2] and shallow net heterotrophy, while lower standardized Hg concentrations in the west (i.e., SWPO and WCPO) were associated with higher subsurface [O2] and deeper net heterotrophy (SI Appendix, Figs. S3 and S4). This pattern suggests that oxygen-poor conditions and shallow bacterial loops both associated with the intense remineralization of organic matter are responsible for enhanced concentrations and/or bioavailability of MeHg at the base of marine food webs.

Marine MeHg concentrations generally peak in low-oxygen thermocline waters, and in the Pacific Ocean, slightly higher MeHg concentrations were measured in the suboxic oxygen minimum zone of the Peruvian upwelling system (12). Although most known Hg methylators are strict anaerobic microorganisms (e.g., sulfate- and iron-reducing bacteria) (29), there is abundant evidence that Hg can be methylated in oxic marine waters (3032), potentially due to anaerobic microniches in otherwise oxygenated environments that support methylation (33). Our results indicate that high seawater MeHg concentrations are not only induced by ambient low-oxygen conditions but also by more intense and shallower remineralization that could be responsible for generating anaerobic microniches with active methylation. Indeed, when exploring observed depths of seawater peak MeHg concentrations in relation to our combination of biogeochemical and physical model estimates with GAMs, variables selected in the optimal model were [O2] in subsurface waters and the depth of net heterotrophy (SI Appendix, Fig. S5). Also consistent with this hypothesis is the fact that a shallower peak of MeHg concentrations was measured in the eastern part of the equatorial Pacific Ocean compared to the western part (34). Knowing that skipjack tuna remain in the upper 100 m and occupy the surface trophic niche throughout the Pacific Ocean, a shallower peak of MeHg concentrations in the EPO would contribute to higher MeHg bioavailability in surface pelagic food webs in this region, which translates here into skipjack spatial Hg distribution. Overall, these correlations between tuna Hg concentrations and subsurface [O2] alongside the depth of net heterotrophy strongly suggest that the biogeochemical conditions of the Pacific Ocean naturally induce large spatial variability of Hg concentrations in marine fish.

### Anthropogenic Releases Enhancing Tuna Mercury Concentrations in the NWPO.

In addition to fish length and biogeochemical and physical variables, Hg0 model estimates were selected in the optimal GAM (Fig. 4D). Higher Hg concentrations in tuna, found in the NWPO and, to a lesser extent, in the CNPO and NEPO, were associated with higher atmospheric Hg0 model estimates (Fig. 4D), suggesting a link between atmospheric Hg emissions and Hg in tuna in the northern Pacific. In particular, the highest standardized Hg concentrations in the NWPO occurred under the Asian atmospheric outflow, characterized by high levels of Asian anthropogenic Hg emissions (35). Knowing the very shallow heterotrophic conditions in the NWPO, one could argue that the biogeochemical processes alone (i.e., higher seawater MeHg concentrations combined with a shallower MeHg peak), as explained previously, are responsible for higher skipjack Hg concentrations in this oceanic region. To account for these potential confounding effects, we defined a modified version of our optimal GAM deliberately excluding Hg0 model estimates (i.e., using only fish length, depth of net heterotrophy, and subsurface [O2]) and used it to predict skipjack Hg concentrations across the Pacific. When mapping the difference between observed and predicted Hg concentrations (Fig. 5), tuna Hg concentrations appeared to be correctly predicted in the entire Pacific except in the NWPO, which highlights that biogeochemical processes alone, combined with fish length, do not allow predicting the spatial variability of observed skipjack Hg concentrations in this region. In the NWPO, GAM output values indeed underestimate the skipjack Hg concentrations compared to what is observed. This highlights the presence of an additional contribution most likely governed by anthropogenic Hg inputs along the Asian coasts, in particular in the East China and Yellow Seas. Enhanced Hg concentrations in tuna from the NWPO are in accordance with enriched Hg levels in both seawater and sediments from the same area (3638).

Fig. 5.

Relative contribution of current local anthropogenic emissions to tuna mercury concentrations. Smoothed spatial contour maps of the difference between observed and GAM predicted mercury concentration (micrograms ⋅ grams−1, dw) in skipjack white muscle samples from the Pacific Ocean. The GAM included fish length, depth of heterotrophy, and subsurface [O2], explaining most of tuna mercury variability, except an Asian outflow mercury contribution in the NWPO. The black dots represent the location of skipjack samples. Ocean areas correspond to the sample origin: NWPO, CNPO, NEPO, EPO, SWPO, and WCPO.

It is worthwhile to mention that we did not consider the potential influence of riverine inputs in our modeling approach. These inputs, although generally considered negligible in the open ocean (39), deliver significant amounts of anthropogenic Hg near coastal areas, especially in the Eastern China and Yellow Seas (3, 37). Nevertheless, we consider that Hg0 model estimates are good surrogates to investigate the anthropogenic contribution of Hg in skipjack tuna, since high atmospheric Hg levels near industrialized landmasses should also translate into enhanced river inputs. Also important when discussing the burden of anthropogenic Hg in tuna is the contribution of legacy Hg to current marine predator Hg levels. Legacy Hg is the stock of historically emitted Hg that is presently stored in terrestrial soils and in the intermediate and deep ocean. The legacy of past North American and European anthropogenic emissions and deposition are estimated to account for about 30% of global anthropogenic Hg in the surface ocean (40) and to be also responsible for high modern Hg0 concentrations in the atmosphere (41). This suggests that despite the clear regional contribution of currently enhanced Hg inputs from the Asian Pacific region, skipjack tuna from the Pacific in general also display a fingerprint of legacy Hg originating from other historic contributing regions that continue to supply the atmosphere and surface ocean with legacy Hg. Yet the surface ocean is a small reservoir (2,600 Mg) of similar size to the atmosphere (4,400 Mg) (1) and has been pointed out to be particularly sensitive to recent anthropogenic emissions (42), suggesting that recent Asian emissions may be partly responsible for the observed Hg hotspot in tuna from the NWPO. Overall, our findings suggest that pelagic biota do respond to recent anthropogenic Hg release but with a significant effect only detectable at the regional scale in the NWPO.

### Implications for Mercury Monitoring in the Global Ocean and Marine Biota.

This study provides an assessment of broad-scale and high-resolution spatial patterns of Hg concentrations in skipjack tuna from the Pacific Ocean. Our spatial interpolation of length-standardized Hg concentrations revealed two primary areas where tuna are relatively enriched in Hg in the Pacific Ocean: the northwest near Asia, and, to a lesser extent, the upwelling regions of the NEPO and southeast Pacific off the coasts of North and South America. The complementary use of ecological tracers alongside biogeochemical and physical variables showed 1) a marine biogeochemistry control, likely via variations in the concentrations of MeHg and the depth at which peak concentrations exist within the water column, which drives the spatial variability of skipjack Hg across the Pacific Ocean, and 2) a cumulated local anthropogenic Hg release effect enhancing tuna Hg levels along the Asian coasts.

Although Hg concentrations in skipjack were low compared to other tropical tuna species and were below the upper threshold of food safety guidelines (1 µg ⋅ g−1, wet weight [i.e., 3.3 µg ⋅ g−1, dw]) (43), our results raise questions regarding human exposure to Hg through the consumption of tunas. Globally, skipjack is the third most consumed marine fish and the most consumed tuna species (20), and the Pacific Ocean remains the primary source of tunas, representing more than 65% of global catch (44). Skipjack Hg concentrations along the Asian coasts mirror atmospheric Hg0 concentrations, and therefore anthropogenic Hg emissions and releases, illustrating a link between recent Asian anthropogenic Hg emissions, marine pelagic fish Hg levels, and human MeHg exposure. In the eastern region, the observations of growing low-oxygen zones (45, 46) and their predicted climate-driven expansion in the coming decades (47, 48) may point to more favorable conditions for MeHg formation. However, parallel changes in primary production and organic matter export, which are key components driving MeHg production at depth, could alter this trend. Although global models tend to project declines in both the production of organic matter by primary producers and its export to the ocean interior in low-latitude oceans over the present century, there are significant uncertainties associated with a range of poorly constrained processes, especially for the EPO (49, 50). Further research is needed to better constrain the integrated impacts of climate change, ocean deoxygenation, and organic matter export on MeHg production and, by extension, MeHg concentrations in tunas.

Our results indicate that Hg levels in skipjack tuna reflect regional anthropogenic atmospheric Hg0 distribution, suggesting they can be used as sentinels of Hg pollution, complementary to direct observations of atmospheric deposition and seawater Hg concentrations, and could therefore provide insights for the future design and implementation of large-scale Hg biomonitoring efforts to evaluate the effectiveness of the Minamata Convention. Given the expected different lifetimes of Hg in the atmosphere (i.e., up to several months) versus in oceans and marine biota (i.e., decades or even thousands of years), these complementary measurements of Hg would be particularly valuable to better evaluate the effect of legacy and recent anthropogenic Hg releases. On the other hand, the importance of natural biogeochemical processes driving MeHg production and enhanced tuna Hg concentrations brings additional information on the potential response of the marine Hg cycle to climate change (51). Further large-spatial-scale studies, combined with Hg environmental data and biogeochemical variables, are needed for other common tuna species (e.g., yellowfin and bigeye) to investigate if a similar response is present in other highly consumed emblematic pelagic fishes. Similarly, complementary large-scale spatial patterns of Hg concentrations in tunas from the Atlantic and Indian Oceans would be valuable for understanding the control factors of Hg in marine predators at the global ocean scale.

## Materials and Methods

### Mercury Concentrations and Ecological Tracers in Skipjack Tuna.

A total of 661 samples of skipjack tuna with Hg concentration data were compiled from published and unpublished studies (Fig. 1A and SI Appendix, Table S1). Samples were collected onboard commercial (longline, purse seine, troll, and pole-and-line) and recreational fishing boats from 122°E to 73°W and 24°S to 45°N, and metadata (i.e., dates, position, and gear type) were provided for each sample. Fork length (FL) was measured and ranged from 29 to 90 cm (58 ± 12 cm; mean ± SD). For each fish, a white muscle tissue sample was collected and stored frozen prior to analyses. Data were classified by oceanic areas: NWPO, CNPO, NEPO, EPO, SWPO, and WCPO (Fig. 1). Sampling spanned from 1997 to 2018, with no temporal sampling bias between oceanic areas (SI Appendix, Fig. S6).

Total Hg concentrations were measured using homogenized freeze-dried samples in different study-specific laboratories, with laboratory-specific reference standards (SI Appendix, Table S1). For the majority of samples (n = 528, ∼80% of the total number of samples), total Hg concentrations were obtained by hot-plate acid digestion (HNO3-H2O2) followed by cold vapor atomic fluorescence spectroscopy. For the rest (n = 133, ∼20% of the total number of samples), total Hg concentrations were measured by thermal decomposition, gold amalgamation, and atomic absorption spectrometry (DMA-80, Milestone). Blanks and biological standard reference materials were routinely used in each analytical batch to check Hg measurement accuracy and traceability. Hg contents are expressed on a dw basis and are considered to reflect MeHg concentrations, as most of the total Hg (>91%) is in its methylated form in tuna white muscle (8).

To account for potential confounding ecological effects, in particular, changes in trophic position, habitat, and feeding pathway, bulk muscle carbon and nitrogen isotopic compositions were also determined when sufficient material was available (n = 475, SI Appendix, Table S1). Measured in the tissues of consumers, carbon isotopes (δ13C) provide information on feeding pathways and habitat associations, while nitrogen isotopes (δ15N) are used to estimate trophic position, as δ15N values increase predictably between prey and consumer. For each sample, δ13C and δ15N values were measured on homogenized freeze-dried samples packed in tin cups and analyzed with an elemental analyzer (online C-N analyzer or Costech elemental analyzer) coupled to an isotope ratio mass spectrometer. Results were reported in the δ unit notation and expressed as parts per thousand (‰) relative to international standards (atmospheric N2 for nitrogen and Vienna Pee Dee belemnite for carbon). Depending on the laboratory, a combination of certified and/or in-house reference materials were used to determine that uncertainty was below 0.2‰ (SI Appendix, Table S1). For samples with elevated lipid content (C:N > 3.5, n = 51), δ13C values were corrected using a mass balance equation with parameters derived from Atlantic bluefin tuna muscle (52).

### Measured Seawater Methylmercury Concentrations in the Pacific Ocean.

Seawater MeHg data from the Pacific Ocean were compiled from the literature from both filtered and unfiltered samples (12, 34, 38, 5357). In particular, we used the maximum total Hg and MeHg concentrations observed in the water column (picomolar) and the depth of maximum MeHg concentrations (meters). When MeHg concentrations were not measured directly from the water samples, we used the sum of monomethylmercury and dimethylmercury.

### Anthropogenic and Natural Environmental Model Estimates Driving Mercury Bioaccumulation.

The global chemical transport model GEOS-Chem, version v11-02 was used to simulate atmospheric Hg0 concentrations. The model couples a three-dimensional atmosphere, a two-dimensional (2D) subsurface-slab ocean, and a 2D terrestrial reservoir at a 6 × 7.5° horizontal resolution. Additional details on model design are given by Angot et al. (58) and Travnikov et al. (59). The simulation was performed using an anthropogenic emissions inventory developed for the year 2015 in the context of the latest Global Mercury Assessment (60). The simulation was driven with assimilated meteorological data from the NASA Goddard Earth Observing System (MERRA-2). In order to account for interannual variability, the simulation was performed with meteorological data for the years 2013 through 2015 following a 1-y spin-up. Annually averaged Hg0 concentrations were extracted from grids for which tuna samples were collected. In order to account for the spatial migration of skipjack tuna, Hg0 concentrations were averaged over adjacent grid boxes. Note that adjacent grid boxes with a fraction of ocean <50% were not taken into account in the average calculation. Model outputs were evaluated against available globally distributed observations (SI Appendix, Fig. S2).

As seawater MeHg data were not available at a large spatial scale and high resolution in the Pacific Ocean, we used a suite of physical and biogeochemical variables known to influence MeHg production in seawater (4, 57, 61, 62) in addition to Hg environmental data. Variables included sea surface temperature (SST, °C), depth of net heterotrophy (meters), depth of oxycline (meters), depth of the mixed layer (MLD, meters), total chlorophyll (Chl-a, milligrams ⋅ meters−3) averaged over the upper 100 m, and vertically integrated net primary production (NPP, milligrams C ⋅ meters−2 ⋅ days−1). Oxygen concentrations in subsurface waters (i.e., 375 m, mmol ⋅ m−3) were also included, as Hg methylation is known to occur in the subsurface layer and as this variable has been identified as the first best environmental predictor of tuna trophic position (63). Baseline estimates included δ15N of particulate organic matter (δ15NPOM, ‰) and δ13C of dissolved inorganic carbon (δ13CDIC, ‰) averaged over the upper 100 m. Selected parameters were extracted from the same hindcast simulation, guaranteeing internal consistency among all of them. Simulations were performed using the isotope-enabled version of the Pelagic Interactions Scheme for Carbon and Ecosystem Studies version 2 biogeochemical model attached to the Nucleus for European Modeling of the Ocean version 4.0 (NEMO) general ocean circulation model (64). The model was forced by the Japanese Reanalysis 55 from 1958 to the present (65). Ocean model resolution is nominally 2° but approaches 0.5° near the equator, while vertical resolution varies between 10- to 500-m thickness over 31 levels between the surface and the abyssal ocean. All outputs on the NEMO curvilinear, tripolar grid were regridded onto a regular 1° × 1° horizontal grid prior to data extraction using the remap function in Climate Data Operators (66). Simulated physical variables were extracted from grid cells nearest to the date and location of tuna samples using the samplexyt_nrst function of the Ferret program (https://ferret.pmel.noaa.gov/Ferret/). Given the long Hg isotopic turnover during bioaccumulation in tuna white muscle (∼700 d; ref. 67), we assumed that Hg values of a given individual captured at a single date and place are not exclusively explained by the environmental conditions at this date but also by prior conditions over previous months or years. As 15N turnover in tuna white muscle is approximately 6 mo (68), we extracted averages of the listed variables above over 6 mo prior to the actual sampling dates, as already done in previous studies on other tropical tuna species (8, 15, 63).

### Spatial Interpolation and Main Drivers of Mercury Concentrations in Skipjack.

Tuna Hg concentrations were first log-transformed to guarantee the homogeneity of variance (69). Furthermore, as Hg is known to bioaccumulate with length (and age), a power-law relationship (

$log(Hg)=a×FLb−c)$

was fit between log(observed Hg) and fish length (FL) to characterize the bioaccumulative processes in skipjack tuna and remove this length effect for further analyses. As the SWPO and WCPO display slow anthropogenic emissions and negligible loadings from anthropogenic sources (60), the power-law relationship was fitted on samples from thes regions only (Fig. 1B). This allowed for investigation of the influence of potential factors, other than length, governing Hg concentrations in tuna and for quantifying the potential Hg enrichment due to anthropogenic emissions. Residuals from the length-based Hg model (i.e., observed values − predicted values) were extracted and used to calculate length-standardized Hg concentrations (at mean skipjack length [i.e., FL = 60 cm]).

A GAM was used to generate smoothed spatial contour maps of observed and standardized Hg levels in skipjack tuna, following the equation

where Y is the expected value of the response variable (i.e., observed and standardized Hg concentrations), α is the model intercept,

$si(Xi)$

is a thin-plate-spline smooth function of the explanatory variable i, and ε is the error term. Here, spatially interpolated maps were generated by fitting 2D thin-plate regression splines on the catch sample location [i.e., s(longitude, latitude)]. Observed and standardized Hg concentrations were assumed to follow a Gamma distribution, and a log link was used in both cases. Spatial interpolation was applied at a 0.25° resolution, with a buffer of 5° around each tuna sample to avoid excessive interpolation where no data were available.

GAMs were also used to test the relative importance of potential predictors on spatial variations of Hg concentrations. The response variable, log(observed Hg), was assumed to follow a Gaussian distribution. Explanatory variables tested included surface (SST, Chl-a, and NPP) and deep (MLD, depth of net heterotrophy, depth of oxycline, and [O2] in subsurface waters) oceanographic variables. Fish length, baseline outputs (δ15NPOM and δ13CDIC), and atmospheric Hg0 estimates were also tested. To account for ecological processes, we also included tuna δ13C values and baseline-corrected tuna δ15N values (i.e., tuna bulk δ15N values − δ15NPOM values) as a proxy for tuna trophic position. Before performing model computation, variance inflation factors (VIF) were calculated between all explanatory variables to detect collinearity. Covariates with the highest VIF were subsequently removed until the highest VIF value was <5 (69). With this method, MLD, Chl-a, and SST were found to be collinear with other variables and were removed as explanatory variables. The remaining explanatory variables (i.e., fish length, δ13C and baseline-corrected tuna δ15N values, δ15NPOM and δ13CDIC, NPP, [O2] in subsurface waters, depth of net heterotrophy, depth of oxycline, and Hg0 estimates) were fitted in the GAM with a low spline complexity (k = 3) to reduce overfitting. A backward selection approach was used, and we chose the model with the lowest Akaike’s Information Criterion corrected for small sample sizes (70). Finally, for the best-fit GAM, assumptions of the residuals trend and autocorrelation were examined graphically with diagnostic plots. The percent DE of the model was compared to assess predictive capacity. All statistical analyses were performed with R 3.6.1 (71), and GAMs were fitted using the mgcv package (72).

## Acknowledgments

We are grateful to our sampling providers, including the many fishery observer programs and fishermen. In particular, we thank the Western and Central Pacific Fisheries Commission Tuna Tissue Bank, the Pacific Marine Specimen Bank managed by the Pacific Community (SPC), and the Environmental Specimen Bank in Ehime University, as well as Robert Olson and Felipe Galvan-Magaña. We thank Laure Laffont from Géosciences Environnement Toulouse (GET) and Jean-Louis Duprey and Stéphanie Berne from Laboratoire des Moyens Analytiques (LAMA) for Hg analyses, as well as the Union College Stable Isotope Laboratory team, Sarah Katz, and Madelyn Miller. Shelagh Zegers provided valuable analytical assistance at Stony Brook University. We thank Laura Tremblay-Boyer from Dragonfly Data Science for her contribution in designing smooth-contour maps in R. H.A. acknowledges Noelle E. Selin and the use of the Svante cluster provided by the Massachusetts Institute of Technology’s Joint Program on the Science and Policy of Global Change. We are grateful to two anonymous referees for providing thoughtful comments on the manuscript. This study was funded by the Pacific Fund VACOPA Project (spatial VAriations of COntaminants levels in PAcific ocean trophic webs) and ANR-17-CE34-0010 MERTOX (unravelling the origin of methylMERcury TOXin in marine ecosystems, 2017–2021) from the French Agence Nationale de la Recherche. The US NSF funded Union College’s isotope ratio mass spectrometer and peripherals (NSF-Major Research Instrumentation No. 1229258). This work was supported by the Interdisciplinary Graduate School for the Blue Planet (ANR-17-EURE-0015) and cofunded by a grant from the French government under the program “Investissements d’Avenir.” P.B. was also supported by a grant from the Regional Council of Brittany (SAD program, Stratégie d’Attractivité Durable). H.A. received financial support from the Swiss National Science Foundation (Grant No. 200021_188478).

## Footnotes

• Accepted November 2, 2021.
• Author contributions: A.M., D.P., and A.L. designed and performed research; A.M. analyzed data; A.M., D.P., T.I., H.A., P.J.B., V.A., L.F., S.G., D.P.G., J.E.S., L.-E.H.-B., M.-M.D., C.E.M., D.J.M., P.B., O.G., A.T., L.B., A.V., and A.L. wrote the paper; D.P. and A.L. supervised and funded the study; T.I., V.A., L.F., S.G., and D.J.M. collected and provided samples; H.A., P.J.B., C.E.M., A.T., and L.B. provided model outputs; A.M., D.P.G., and A.V. analyzed samples; J.E.S., L.-E.H.-B., and M.-M.D. provided compiled Hg concentrations in marine sediments and seawater; and P.B. provided R codes for statistical analyses.

• The authors declare no competing interest.