Mercury stocks and methylmercury production remain unaffected by 17 years of experimental permafrost thaw

Haugk, Charlotte; Azaroff, Alyssa; Bertilsson, Stefan; Johansson, Margareta; Liu, Tong; Thompson, Lauren; Jonsson, Sofi

doi:10.3389/fenvc.2025.1715512

ORIGINAL RESEARCH article

Front. Environ. Chem., 02 February 2026

Sec. Inorganic Pollutants

Volume 6 - 2025 | https://doi.org/10.3389/fenvc.2025.1715512

This article is part of the Research TopicEnvironmental Chemistry of Mercury: Sources, Pathways, Transformations and Impact Vol IIView all 6 articles

Mercury stocks and methylmercury production remain unaffected by 17 years of experimental permafrost thaw

Charlotte Haugk^1,2*

Alyssa Azaroff¹

Stefan Bertilsson³

Margareta Johansson⁴

Tong Liu³

Lauren Thompson^5,6

Sofi Jonsson¹*

¹Department of Environmental Science and Analytical Chemistry, Stockholm University, Stockholm, Sweden
²Bolin Centre for Climate Research, Stockholm University, Stockholm, Sweden
³Department of Aquatic Sciences and Assessment, Swedish University of Agricultural Sciences, Uppsala, Sweden
⁴Department of Physical Geography and Ecosystem Science, Lund University, Lund, Sweden
⁵Department of Renewable Resources, University of Alberta, Edmonton, AB, Canada
⁶Hatfield Consultants, Calgary, AB, Canada

Mercury (Hg) has been stored in permafrost peatlands for millennia. As permafrost thaw is predicted to increase with ongoing climate warming, Hg is at risk to be remobilized from those peatlands and hotspots for Hg methylation could potentially form. Monomethylmercury (MeHg) is a known neurotoxin and a health concern to northern communities if Hg is remobilized, transformed to MeHg and subsequently bioaccumulated in the food chain. It is uncertain how Hg cycles in thawing permafrost systems and how much of it could potentially be remobilized by thaw processes. In this study, we have investigated Hg dynamics in a permafrost peatland of northern Sweden, where a snow fence field experiment was set up in 2005 to simulate accelerated permafrost thaw through winter warming. We compared total mercury (THg) and MeHg concentrations in soil plots representing intact and thawed permafrost conditions, investigated seasonal variations and examined the coupling between microbial community composition and MeHg concentration. Similar stocks of both total THg and MeHg were observed in intact and thawed permafrost conditions, suggesting that 17 years of winter warming manipulation and accelerated permafrost thaw had not led to substantial Hg loss from the peat nor extensive MeHg production. The apparent stability of the Hg stocks contrasted with our hypothesis and with many previous studies. While there was no difference in microbial communities between treatments with or without accelerated thaw, putative methylators were more abundant in thaw plots in the fall. This indicates that permafrost thaw has increased the potential for Hg methylation, although these shifts have not yet been strong enough to measurably affect MeHg stocks. Our study emphasizes the complexity of Hg dynamics in thaw-affected permafrost landscapes and the need to consider thaw-related perturbations of the Hg cycle on various timescales.

1 Introduction

Accelerated permafrost thaw could remobilize currently locked mercury (Hg) and may worsen the Hg pollution in the Arctic. Current estimates suggest between 600 and 1,600 Gg of Hg to be stored in the top 3 m of permafrost soils (Schuster et al., 2018; Lim et al., 2020). Various chemical forms of inorganic Hg (Hg^II) account for most of this Hg, though Hg^II is typically not found in nature at harmful levels. Microbial processes, however, can lead to the formation of monomethylmercury (MeHg), a neurotoxic and organic form of Hg that can bioaccumulate in food webs to levels of severe health concern. In the Arctic, high levels of Hg in human and wildlife are already evident (AMAP, 2021). Despite the potential impact that thawing permafrost systems may have on future levels of Hg exposure, large uncertainties remain in our understanding of Hg cycling in these complex and rapidly changing environments.

The risks associated with Hg in thawing permafrost landscapes are closely linked to the timing and extent of physical changes that occur in the landscape. In past decades, permafrost thaw has accelerated due to enhanced Arctic warming (Rantanen et al., 2022). Thaw processes occur on an annual to decadal scale and can cause varying degrees of permafrost degradation (Grosse et al., 2011) that result in diverse thermokarst landforms. The thaw may be gradual, with permafrost slowly degrading by centimeters over decades, or abrupt involving rapid erosion that occurs within days to years (Turetsky et al., 2020). As these thaw processes reshape the landscape, they can also impact the Hg cycle in permafrost regions, shifting current Hg sinks into sources that release Hg to the atmosphere and downstream waters. St. Pierre et al. (2018) for example, observed increased export of Hg^II and MeHg, mainly associated to remobilized soil particles from an actively eroding thaw slump in Canada. Hg mobilization has also been observed along less dramatic thaw features, either as a loss of Hg from the soil matrix (Klaminder et al., 2008), through hydrological transport to downstream systems (Rydberg et al., 2010a; Schaefer et al., 2020), or via Hg reduction and subsequent Hg(0) emissions to the atmosphere (Fahnestock et al., 2019). In addition to the concerns of Hg mobilization, formation of waterlogged and organic-rich thermokarst environments can result in increased production of MeHg (Fahnestock et al., 2019; Tarbier et al., 2021). The production of MeHg in these environments is depending on a variety of geochemical factors including redox conditions, dissolved organic matter (DOM) quantity and composition, the bioavailability of Hg and microbial community composition (Ullrich et al., 2001; Mitchell et al., 2008; Gordon et al., 2016; Poulin et al., 2019; Gindorf et al., 2025). As collapsed and saturated wetland environments are typically anoxic, they provide ideal conditions for the microbial transformation of inorganic Hg^II to MeHg by bacteria and archaea containing the hgcAB⁺ gene pair (Parks et al., 2013). Dominant methylators in anoxic environments include sulfate reducing bacteria (SRB), iron reducing bacteria (IRB) and methanogenic archaea (King et al., 2000; Bravo et al., 2018). Microbial activity has been shown to increase along a permafrost thaw gradient (Fahnestock et al., 2019), likely promoting the production of MeHg (Podar et al., 2015). These findings highlight the impact of permafrost thaw on Hg dynamics, yet the timeline of thermokarst development and its connection to Hg mobilization upon thaw remains poorly understood.

Existing studies of Hg cycling in permafrost systems have mainly relied on laboratory incubation experiments, natural thaw gradients, and/or modelling. Field experiments, including the use of snow fences, fill an important gap by providing opportunities to test the effects of permafrost thaw within a known timeframe and with replicatication while preserving the complexity of natural environments, but they have not yet been utilized to study Hg dynamics. Snow fences are geoengineering structures that have widely been used to simulate winter warming events and study accelerated permafrost thaw and its environmental consequences in cold regions (Johansson et al., 2013; O’Neill and Burn, 2015; D’Imperio et al., 2018). By trapping and accumulating snow, these structures can alter the ground thermal regime leading to thicker snow cover and greater ground insulation, which over time results in permafrost thaw. Especially in ice-rich permafrost soils, this warming effect causes thermokarst erosion, including ground subsidence, active layer thickening and increased soil moisture followed by accumulation of standing water and pond formation (Hinkel and Hurd, 2006; Johansson et al., 2013). Winter warming events are predicted to become more frequent (Pascual and Johansson, 2022), particularly in the Abisko region, where climate predictions indicate increasing precipitation with snow depths expected to rise by 2–3 cm per decade (Kohler et al., 2006). With both continued warming and increasing snow accumulation, research on the effects of permafrost thaw on the Hg cycle becomes increasingly important. In particular, a better understanding of the extent and the timeframe of permafrost thaw is needed to assess how it will alter the biogeochemical cycle of Hg.

Here, we investigated Hg dynamics in a snow fence field experiment where, at the time of our study, the snow depth had been manipulated for 17 consecutive years. We analyzed and evaluated the concentrations of Hg and MeHg in soil and water collected from manipulated and control plots in both spring and fall, along with microbial communities (including putative methylators), field observations of thaw depth, surface moisture conditions and other ancillary parameters. We hypothesized that increased snow depth, and the associated permafrost thaw, had resulted in the remobilization of Hg from the peat matrix, as well as higher net production of MeHg, to such a degree that the stocks of Hg and MeHg in the peat had been altered. We further hypothesized that changes in net production of MeHg could be linked to the diversity and abundance of putative Hg methylators. By utilizing the now 17-year-old snow fence field experiment to examine Hg cycling in permafrost peatlands, this study provides valuable insight into the timing of permafrost thaw related perturbation of the Hg cycling in northern peatlands.

2 Materials and methods

2.1 Experimental set-up and sample collection

This study was conducted at a snow fence experimental set up in the ombrotrophic mire complex “Storflaket” (68°20′48″N, 18°58′16″E), 5 km east of Abisko Scientific research station in northern Sweden (Figure 1). The mire has a peat cover reaching a thickness of 0.5–0.9 m underlain by silty sediments of glacial origin (Åkerman and Johansson, 2008; Klaminder et al., 2008). The dominant vegetation is consist of peat moss (e.g., Sphagnum spp.), dwarf shrubs (Andromeda polifolia L.,Vaccinium vitis-idaea L., Betula nana L., Empetrum nigrum L.), lichens, graminoids (Eriophorum vaginatum L.) and Rubus chamaemorus (Johansson et al., 2013; Olid et al., 2020). The permafrost peatlands around Abisko are characterized by isolated or sporadic permafrost (Gisnås et al., 2017) and have mean annual ground temperatures close to 0 °C (Johansson et al., 2011). These conditions make the Abisko region especially sensitive to permafrost thaw in comparison to colder and more stable permafrost zones, as a relatively small change in ground thermal regime could lead to accelerated permafrost degradation and thermo-erosion (Johansson et al., 2011).

Figure 1

Map depicting Abisko in Sweden near the Arctic Circle. Abisko is marked with a black dot, adjacent to Torneträsk Lake. Storflaket is marked with a red dot. The map includes a scale bar indicating distances in kilometers and a north arrow for orientation.

Figure 1. Map of sampling area in Abisko, northern Sweden. Left panel shows the location of Abisko in the north of Sweden, and right panel shows Storflaket mire in relation to Abisko south of Torneträsk Lake. Map created by the author using ArcGIS Pro 3.3.

The field experiment was set up in 2005 to simulate ongoing permafrost thaw with increased snow thickness and has been described thoroughly in Johansson et al., 2013. Briefly, snow fences were erected on 6 out of 12 randomly chosen plots in the western part of the mire (Supplementary Figure S1A). By accumulating snow on top, the ground is insulated and experiences warming, simulating permafrost thaw. These six plots are further referred to as “manipulated”; the remaining six plots without snow fences serve as “control” plots. A short description of studies connected to this experiment can be found in the Supplementary Material. We collected one core from each plot in spring (11th-12th of June 2022, after snow melt) and in fall (3rd September 2022, within the period when thickest active layer is expected). Sampling was performed at 5 cm fixed depth intervals and conducted manually to minimize disturbance of the experimental site. The surface layer was collected by sawing peat blocks with a knife or scissors to avoid compressing the upper peat layer. For the rest of the profiles in the wet plots, we used a Russian peat corer with a 50 cm long chamber. For dry sites, we used a shovel to continue the soil profile and sampled from the side of the pit. To extract permafrost samples, we used a steel pipe by hammering it into the ground in 5–10 cm depth increments. Total length of cores was reached once we sampled the first 10–20 cm of permafrost table (in spring) or reached the mineral layer (in fall). Soil samples for mercury speciation and elemental analysis were subdivided in the field, transferred to zip-lock bags and kept cold during transport to the Abisko Research station where all peat samples were frozen. Prior to chemical analysis, samples were freeze-dried in a Heto Drywinner 6.55 freeze-drier and then homogenized using either an agate mortar and/or a scientific grade coffee mill (DeLonghi KG40). Bulk density (g cm^-3) was calculated by dividing the dry weight (d.w. in g) by the volume of sample (cm³).

For microbiological analysis, aliquots of soil samples were taken from the same cores (spring n = 34; fall n = 27) from two or three depths in each plot. In all plots, a surface sample (“top”, at 5–10 cm depth) was taken alongside an active layer sample (at 20–25 or 45–50 cm), later on called “above permafrost table”. If permafrost was reached in a plot, the first permafrost depth below the permafrost table was sampled as well (“below PFT”, at ca 45–50 cm or 60–65 cm).

Water samples for chemical analysis were collected from ponds that formed on the plots (spring: n = 6; fall: n = 7) using a sterile and pre-washed syringe. Water samples were filtered (<0.45 µm polyethersulfone (PES) syringe filter) and acidified (with 0.5% v:v with trace metal grade hydrochloric acid) in the field and collected in burned 60 mL amber glass bottles, 120 mL PET bottles and 50 mL PE falcon tubes for further analysis of THg, MeHg and dissolved organic carbon (DOC), respectively. Unfiltered water samples were collected in 50 mL PE falcon tubes for determining pH and total organic carbon (TOC). All water samples were stored cold during transport until further analysis.

2.2 Soil analysis

2.2.1 Total mercury

Approximately 20 – 50 mg of freeze-dried and homogenized sample material was weighted and measured for total Hg concentration by direct thermal decomposition, amalgamation and atomic absorption spectrophotometry using a direct mercury analyzer (DMA-80, Milestone) and following US EPA Method 7473. The instrument detection limit was 0.01 ng HgT, with a working range of 0.05–600 ng. Accuracy and precision of the method were determined using an internal reference material of estuarine sediment SE-3 (n = 35, 210 ± 10 ng g^-1) after every 10 samples analysis. To track the risk of memory effects between samples an empty sample boat was also added every 10 samples. Recoveries and standard deviation (SD) of the internal reference material (210 ng g^-1) were 205.7 ± 18.2, with a relative standard deviation (RSD) of 9%. Triplicate analysis of samples (n = 24) showed a SD of 2.3 %RSD at 7% on average (Supplementary Table S1).

2.2.2 Methylmercury

MeHg was first extracted via solid-liquid extraction following Qian et al. (2000). Briefly, 0.2 – 1 g of sample material was weighted into 50 mL falcon tubes and 100–300 µL of an isotopically enriched Me²⁰¹Hg standard solution (concentration at 1.308 μL L⁻¹, in Milli-Q (MQ) water) was added as an internal standard and left to equilibrate for 1 hour. The amount of spike was adjusted for each sample according to the expected concentration. Subsequently, 10 mL of potassium bromide (1.4M), 2 mL copper sulfate (2M) and 10 mL dichloromethane (DCM) were added to each tube. After 45 min, samples were rotated at 65 rotations per minute (RPM) for 45 min and centrifuged for 5 min at 3000 RPM. The supernatant was then manually transferred to a 50 mL falcon tube using a glass Pasteur pipette. After adding 10 mL of MQ water, the DCM was purged with nitrogen gas while in controlled water bath at 40 °C (Lambertsson et al., 2001). The extracts were kept at 4 °C (if analyzed wihtin 24 h) or kept at −20 °C until analysis.

In a second step, the extract of soil samples (as well as pondwater samples) was prepared in 40 mL amber glass vials for MeHg analysis. Depending on the expected concentration, between 0.5 and 2 mL of soil extract or 15 mL of pondwater was added into MQ water and then adjusted to a final volume of 30 mL. The extract was ethylated via NaBEt₄ (sodium tetraethyl-borate) at an optimal pH 4.9 (using 225 µL 2 M acetate buffer) (Munson et al., 2014). MeHg was then quantified using a Tekran Model 2700 Automated Methylmercury Analyzer coupled with Teflon tubing to an Inductively Coupled Plasma Mass Spectrometer (ICPMS, Thermofisher X-series 2). The concentration of MeHg for each sample was calculated through signal deconvolution using a mass-bias corrected signal derived approach after (Qvarnström and Frech, 2002). Accuracy and precision of the method were determined with certified reference material (n = 16; ERM-CC580, 75 ± 4 ng g^-1). Replicates of the reference material had a mean and SD of 73.1 ± 6.2 ng g^-1 with RSD of 8%. Replicates of soil samples (n = 16) had a SD of 0.03 and %RSD of 11% on average (Supplementary Table S2). Blanks were extracted each new day alongside samples to ensure no contamination of MeHg during the extraction process (n = 10; 0.015 ± 0.0045 ng g⁻¹).

2.2.3 Elemental analysis

All samples were analyzed for carbon and nitrogen with a Thermo FlashSmart Elemental Analyzer. 2–3 mg of soil were packed into tin capsules and then flash combusted at 950 °C. Compound peaks were first integrated by EagerSmart software and adjusted manually after peak inspection. The instrument was calibrated using 0.1 – 3 mg of 2.5-bis (5-tetra-butyl-benzoxazol-2-yl) thiophene (BBOT) at the start of every run and every 10 samples 0.5 mg of Methionine and 2 mg of Acetanilide for quality control were tested (Supplementary Table S3).

2.2.4 Microbial community measurements

DNA extraction was performed using the DNeasy PowerSoil Pro Kit (Qiagen), and DNA concentrations were measured using a Qubit fluorometer (Invitrogen). Bacterial 16S rRNA gene libraries were constructed by amplifying the V3–V4 hypervariable region with polymerase chain reaction (PCR) using primers 341F (CCTACGGGNGGCWGCAG) and 805NR (GACTACNVGGGTATCTAATCC) (Herlemann et al., 2011). The PCR and library preparation protocols were adapted from Sinclair et al. (2015) with slight modifications. For the first PCR, 1 µL of 10-fold diluted DNA templates (average concentration 25.4 ng μL⁻¹, range: 8.4–36.2 ng μL⁻¹) were each combined with 0.2 µL of Q5 High-Fidelity DNA Polymerase (2 U µL⁻¹; New England Biolabs), 2 µL of 2 mM dNTPs, 0.5 µL each of forward and reverse primers (5 µM), and 15.8 µL of nuclease-free water in a 20 µL reaction. The thermal cycling program included an initial denaturation at 98 °C for 3 min, followed by 20 cycles of denaturation at 98 °C for 10 s, annealing at 48 °C for 30 s, and extension at 72 °C for 30 s, with a final 2 min extension at 72 °C. PCR products were verified for size and quality via agarose gel electrophoresis and purified using Seramag Select magnetic beads as recommended by the manufacturer (Cytiva).

For the second PCR, Illumina sequencing adapters and unique index primers were added to the amplicons. A reaction mix containing 4 µL of the primary PCR product, 0.2 µL of Q5 High-Fidelity DNA Polymerase, 2 µL of 2 mM dNTPs, 1 µL each of forward and reverse index primers (5 µM), and 11.8 µL of nuclease-free water was used in a 20 µL reaction. The thermal cycling protocol involved an initial denaturation at 98 °C for 30 s, followed by 15 cycles of denaturation at 98 °C for 10 s, annealing at 66 °C for 30 s, and extension at 72 °C for 30 s, with a final extension at 72 °C for 2 min, then held at 6 °C. The final products were confirmed by agarose gel electrophoresis and purified using the same magnetic bead protocol. DNA concentrations of purified products were quantified, and equal amounts (20 ng) of each sample were pooled to create the sequencing library. The library was sequenced using the Illumina MiSeq platform with Reagent Kit v3 (2 × 300 cycles). Raw sequencing data were processed using the DADA2 pipeline (version 1.16) in R (version 4.0.2) (Callahan et al., 2016). Forward and reverse primers were removed, and sequences were trimmed to lengths of 277 bp and 245 bp, respectively, with quality control thresholds set at maxEE=(2, 2) and truncQ = 2, based on FIGARO predictions (Weinstein et al., 2019). Forward and reverse reads were merged with minimum overlap = 10 bp, maximum mismatch = 0. Taxonomic assignment of amplicon sequence variants (ASVs) was performed using the SILVA rRNA database (release 138.1) (Quast et al., 2012). Samples with sequencing depths below 1,000 reads were excluded from further analyses.

2.3 Pondwater analysis

2.3.1 Total mercury and methylmercury

Total Hg analysis in water was performed following EPA-method 1631 by oxidation, purge and trap and cold vapor atomic fluorescence spectrometry (CV-AFS) using a Tekran 2600 Total Hg analyzer. Water samples were analyzed in burned 40 mL glass vials, Briefly, 0.125 mL of bromine chloride (BrCl, 0.5% v:v) was added to each sample, to allow Hg compounds to be digested and oxidized overnight. Remaining halogens were then removed by adding 0.03 mL of hydroxalamine hydrochloride (H.H, 0.2% v:v) and after 10 min 0.06 mL of tin chloride (SnCl, 30%) was added to reduce Hg to Hg⁰. On a gold trap Hg⁰ was collected and desorbed 550 °C for detection by atomic fluorescence. MeHg in pondwater was analyzed in the same manner as for the soil extracts (Tekran 2700 coupled to ICPMS), as described above. For MeHg measurements limit of detection (LOD) and limit of quantification (LOQ) were 0.05 and 0.16 in peat and 0.005 and 0.02 in pondwater respectively.

2.3.2 Ancillary parameters

Both DOC and TOC of the pondwater were measured with a total organic carbon analyzer (TOC-L, Shimadzu, Japan). Briefly, 3 mL of the sample was diluted in Milli-Q water to a total volume of 30 mL (dilution factor of 10) and analyzed in 40 mL glass vials. Samples and blanks that were not previously acidified in the field, were acidified with HCL to 0.5% v-v and measured in triplicates. pH was measured using a pH electrode (Metrohm AG, Switzerland).

2.4 Statistical analysis

Normality of the data distribution was tested by using the Shapiro−Wilk test prior to further analysis using Anderson-Darling test for goodness of fit. If data were not normally distributed (e.g., THg, MeHg and %MeHg data for peat) it was adjusted by using log10(x+1) transformation. Differences between treatments (manipulated, control), surface moisture condition (dry, inundated) or season (spring, fall) were tested for significance using One-way analysis of variance (ANOVA) with Student’s t-test. All statistical analyses were done with JMP Pro (version 17.2.9) at a significance level of p < 0.05. The ratio of THg to C (R_HgC) was calculated to estimate how much Hg relative to C is bound in the soil matrix. Concentrations of THg, MeHg, C and N are reported as per the dry weight of soil, unless otherwise stated.

2.4.1 Calculation of stocks

Mercury and carbon stocks were calculated for each core by multiplying the THg (ng Hg g⁻¹), MeHg (ng MeHg g⁻¹) and C (g C cm⁻³) concentrations with the dry bulk density (DBD, g cm⁻³) for the respective depth (cm), respectively. Values were integrated to the underlying mineral layer acting as baseline for comparison of peat profiles between cores. If the mineral layer was not reached, the peat thickness was assumed to be 0.90 m according to previous observation of the mire (Åkerman and Johansson, 2008). The stocks of THg, MeHg and C were expressed in mg THg m^-2, mg MeHg m^-2 and kg C m⁻².

2.4.2 Microbiological multivariate analysis

Prior to microbiological multivariate analysis, all samples were rarefied to the lowest read count (3,216) using the rrarefy function from the vegan package (version 2.6-6.1). Permutational multivariate analysis of variance (PERMANOVA) was conducted using the adonis2 function in the vegan R package (version 2.6-6.1) to assess which environmental variables significantly explained variation in microbial community composition (based on Bray–Curtis dissimilarities, 999 permutations). Microbial community dissimilarities were visualized with non-metric multidimensional scaling (NMDS) based on Bray-Curtis distances of rarefied ASV read counts, using the Phyloseq package (version 1.44.0).

To identify taxa contributing most to between-group differences, we applied SIMPER (similarity percentages) as a descriptive follow-up to PERMANOVA for contrasts of interest, season (spring vs. fall) and treatment (manipulated vs. control) evaluated overall and within each season (spring, fall). SIMPER was run on Bray–Curtis dissimilarities computed from phylum-level relative abundances aggregated from the rarefied ASV table (vegan 2.6-6.1; simper). For each contrast, between-group dissimilarity was partitioned among phyla and summarized as mean ± SD percent contribution with a consistency index (mean/SD); the direction of change was inferred from group means (e.g., fall−spring; manipulated−control), and contribution significance was evaluated with 999 label permutations. In addition, for phyla highlighted by SIMPER we quantified between-group differences in mean relative abundance (in percentage points) and reported percentile bootstrap 95% CIs (5,000 resamples) and two-sided permutation p-values (20,000 label permutations); when multiple phyla were assessed simultaneously, p-values were adjusted using Benjamini–Hochberg false discovery rate (FDR).

The dataset was then partitioned into subsets for further analysis of microbial population effects using linear discriminant analysis (LDA) implemented in the MicrobiotaProcess package (version 1.12.3) (Xu et al., 2023), to detect differences in microbial relative abundances across seasons and treatments (including within-season contrasts). The LDA workflow consisted of two stages: (i) Kruskal–Wallis to identify taxa differing across groups (p ≤ 0.05), followed by (ii) Wilcoxon pairwise tests with Benjamini–Hochberg FDR correction. Correlation analysis among environmental variables was conducted using corrplot (version 0.92) and phyloseq (version 1.44.0) (McMurdie and Holmes, 2013). All 16S rRNA data analyses were performed in R 4.4.2 (R Core Team, 2024).

3 Results

3.1 Field observations

Measured thaw depths in the spring ranged from 0.22 to 0.36 m and 0.28 – 0.72 m in control and manipulated plots, respectively (Figure 2A; Table 1). In the fall, the permafrost table was not detected in the upper 1 m soil in the manipulated plots, while thaw depth in control plots reached between 0.59 and 1.05 m. For five of the six manipulated plots, ponds (ca. 6–10 m in diameter) were present around the length of snow fences both in the spring and fall. Ponds were also present in two of the control plots (one in both fall and spring and one in fall only; Figure 2B). The remaining plots (one manipulated and four control plots) were dry at both sampling occasions. At least one pond had clearly expanded from spring to fall and started to expand into another plot (Supplementary Figure S1B). The ponds were marked by sedges and/or open water with bryophyte (e.g., Sphagnum spp.) cover (Supplementary Figure S2). Although we did not measure the water depth at the centre of the ponds, we consider all of them shallow ponds with maximum depth of around 0.5 m. In spring, we noticed shallower thaw depths (0.23 – 0.28 m) in plots without ponds, while inundated plots with ponds had deeper thaw depths (0.22 – 0.72 m; Supplementary Figure S3). As the conditions in manipulated plots (formation of thaw ponds) and control plots (dry conditions) were not consistently presented in all plots, we compared Hg dynamics both between the treatments (manipulated, control) as well as between surface moisture conditions, (plots with and without thaw ponds present, hereon referred to as inundated and dry, respectively).

Figure 2

Line and bar graphs display thaw depth data. In graph A, lines for control (blue) and manipulated (orange) show depth from 2005 to 2022, with data sources noted. Graph B compares control and manipulated thaw depths by season, using bars representing peat, mineral, and pond formation. Thaw depth and core loss are marked, with asterisked notes for thaw depth over one hundred centimeters.

Figure 2. (A) Active layer thaw thickness (cm) averages and standard deviation over time since the start of the experiment in 2005 until this study in 2022. Data from respective publications show fall measurements only at the time of maximum active layer thickness. Empty diamonds in 2022 represent measurements below 100 cm (control: n = 1, manipulated: n = 6), while active layer thickness was reached in five control plots. (B) Sampling schematic of all 12 experimental plots where samples in control (left panel) and manipulated plots (right panel) are marked in orange for peat and grey for mineral layers. Cores from the same plots are shown next to each other for sampling in spring and fall 2022. Black solid lines are thaw depths. If > 100 cm, the active layer was thicker than the length of the probe (*). Ponds that formed on plots were subsequently sampled for chemical analysis of total mercury (THg), methylmercury (MeHg), dissolved organic carbon (DOC) and total organic carbon (TOC).

Table 1

Table 1. Concentration of total mercury (Hg, ng g⁻¹), methylmercury (MeHg, ng g⁻¹) and carbon (C, %), %MeHg, R_HgC (μg Hg g C⁻¹) and carbon to nitrogen ratio (C/N)as well as stocks for THg (mg Hg m⁻²), MeHg (mg MeHg m⁻²) and C (g C m⁻²) of peat profiles in intact permafrost conditions (control) versus thawed conditions (manipulated) presented as mean with standard deviation (mean ± SD). Plot ID derived from original publication (Johansson et al., 2013). Each plot is described by treatment type, season, surface moisture condition, active layer thickness (ALT, m) and peat depth (m). Note that plot “C2” was dry in spring but had formed a pond in the fall (Figure 2B). If we could not detect the permafrost table with our equipment, active layer thickness was > 1 m.

3.2 Mercury and ancillary parameters

3.2.1 Peat

Concentrations of THg and MeHg in peat samples ranged from 1.5 to 270 ng g⁻¹ and from 0.030 to 2.6 ng g⁻¹, respectively (Table 1). No consistent trend with peat depth was observed among the cores for either THg or MeHg; however, THg often peaked within the top 5–30 cm of both control and manipulated cores, which is also shown by the higher average concentration in this section (Figures 3A,B; Supplementary Figure S4). We found no significant difference between treatments (manipulated vs. control), surface moisture condition (inundated vs. dry) or season (fall vs. spring) for THg, MeHg concentration, %MeHg or R_HgC in the top 30 cm of the peat or the full peat profile (Supplementary Tables S4–S6). THg concentrations positively correlated with %C (R² = 0.13, p < 0.001, Supplementary Figure S6A) and the R_HgC ranged from 0.0032 to 0.60 μg Hg g C⁻¹ in peat samples. The concentration of MeHg correlated with THg and %C (p < 0.05, Supplementary Figures S6B,C). The stocks of THg and MeHg were on average 7.9 ± 3.4 mg Hg m^-2 and 0.029 ± 0.012 mg MeHg m^-2, respectively, in control plots and 7.4 ± 2.7 mg Hg m^-2 and 0.059 ± 0.048 mg MeHg m^-2, respectively, in manipulated plots (Table 1). We could not observe any statistically significant difference of THg and MeHg stocks between treatments, surface moisture condition or seasons when comparing calculated stocks of THg and MeHg (Supplementary Table S5). The same was true for both stocks of C and N. Only when comparing stocks of THg and C in the fall was a statistical difference between control and manipulated treatment observed.

Figure 3

Four graphs showing the depth (in centimeters) versus concentration of total mercury (THg) and methylmercury (MeHg) in nanograms per gram, comparing spring and fall seasons. Panels (A) and (B) display THg, while panels (C) and (D) show MeHg. Spring data is marked in blue, while fall data is in red, with shaded areas representing variability.

Figure 3. Vertical depth distribution of total mercury (THg) concentration (ng g⁻¹) in (A) control and (B) manipulated plots and for methylmercury (MeHg) concentrations (ng g⁻¹) in (C) control and (D) manipulated plots. Lines represent averages and shaded areas standard deviations for each season respectively (spring in blue, fall in red).

3.2.2 Pondwater

The concentration of dissolved THg, MeHg and DOC ranged from 4.4 to 23 ng Hg L^-1, 0.025 and 2.6 ng MeHg L^-1 and 30–150 mg C L^-1, respectively (Table 2). Pondwater concentrations did not differ between treatments, noted that ponds mainly formed on manipulated plots (n_manipulated = 10, n_control = 3). We did not find any difference between filtered and unfiltered waters for concentrations of THg or MeHg, implying that THg and MeHg exists predominantly in the dissolved form in pondwater. The %MeHg positively correlated with the concentration of MeHg (p < 0.05), and %MeHg showed a weak negative correlation with THg (p = 0.07) and DOC (p < 0.05) (Supplementary Figure S7). Dissolved MeHg and THg concentrations did not correlate. With values ranging from 0.27% to 31%. %, MeHg did not show significant seasonal variation. Notably, three ponds had exceptionally high %MeHg – one in spring and two in fall – that were 130–150 times higher than the lowest observed value (Table 2). Two ponds measured in fall had exceptionally low %MeHg (<1%, C2 and C3, both control plots, Table 2). The mean pH was higher in spring compared to fall (p<=0.001), while we observed higher concentrations of TOC (p < 0.05) as well as DOC (p < 0.05) in fall than in spring. When concentrations of Hg, MeHg and %C of the surface peat (0–30 cm) were tested against pondwater concentrations, we did not find any correlation (Supplementary Figure S7).

Table 2

Table 2. pH, THg (ng L⁻¹) unfiltered and filtered, MeHg (ng L⁻¹) unfiltered and filtered, dissolved organic carbon (DOC, mg L⁻¹), total organic carbon (TOC, mg L⁻¹) and %MeHg (for filtered waters) in pondwater samples of each plot for spring and fall. Values are shown as means with standard deviation (mean ± SD).

3.3 Microbial community structure

3.3.1 Microbial community composition in the peat profiles

Illumina sequencing of 16S rRNA gene libraries from 54 peat samples produced a total of 3,080,318 reads. After quality filtering, paired-end merging, and chimera removal, 1,520,967 high-quality sequences remained. Nine samples containing fewer than 1,000 reads after quality control were excluded from further analysis, leaving an average of 34,567 reads per sample, ranging from 3,126 to 128,717. No sample had more than 6% of merged reads that could not be classified taxonomically at the kingdom level. Additional details on read filtering, sequencing quality, and overall sequencing depth are provided in the Supplementary Table S7.

The 16S rRNA gene profile revealed that the peat sample microbial communities were predominantly composed of the bacterial phyla Acidobacteriota, Proteobacteria, and Actinobacteriota, which accounted for 34.9% ± 14.9%, 17.8% ± 9.1%, and 15.7% ± 9.8% of the total reads, respectively, followed by Verrucomicrobiota, Patescibacteria, Bacteroidota, Planctomycetota, WPS-2, Firmicutes, and Chloroflexi, which average relative abundances ranging from 7.1% to 1.1% (Supplementary Table 1). Among the environmental variables, season and sample depth were the most influential factors for community composition (PERMANOVA, p < 0.01), while no clear differences were seen by the different treatments (PERMANOVA, p > 0.05). Further alpha diversity analysis of the 16S rRNA data revealed that microbial richness (both estimated and observed), as well as Shannon diversity and the Simpson index, were significantly higher in the fall compared to the spring (p < 0.01, Supplementary Figure S8). Beta diversity difference between seasons was evident in NMDS ordinations of ASV-level relative abundances (Figure 4). Consistent with the NMDS separation, SIMPER partitioning of Bray–Curtis dissimilarities indicates that seasonal divergence is driven primarily by shifts in a few phyla (∼73% of spring–fall dissimilarity): Acidobacteriota (−8.07 percentage points, 95% CI [−16.53, 0.54], p > 0.05), Actinobacteriota (−5.96 pp, [−11.26, −0.76], p < 0.05), and Proteobacteria (−2.57 pp, [−8.61, 2.93], p > 0.05) are higher in spring, whereas Verrucomicrobiota (+3.83 pp, [0.75, 6.59], p < 0.01), Patescibacteria (+5.15 pp, [2.46, 8.06], p < 0.01), and Bacteroidota (+5.20 pp, [3.10, 7.55], p < 0.01) increase in fall (Figure 5; Supplementary Table S2). At lower taxonomic ranks, the pattern persists with within-phylum heterogeneity: Acidimicrobiia (Actinobacteriota) and Acetobacteraceae (Proteobacteria) show higher relative abundance in spring, whereas Pedosphaeraceae (Verrucomicrobiota), Parcubacteria (Patescibacteria), and Burkholderiales lineages such as Rhodoferax are higher in fall (Supplementary Table S3).

Figure 4

NMDS scatter plot showing two dimensions (NMDS1 and NMDS2) with points representing data from fall (red) and spring (blue) across control (circles) and manipulated (triangles) treatments. Vectors indicate relationships with MeHg ng/g, THg ng/g, MeHg%, C%, and N%. Stress value is 0.17.

Figure 4. Non-metric Multidimensional Scaling (NMDS) of amplicon sequence variants (ASVs) of the 16S rRNA genes in the core samples. Dots colored according to the sample seasons and shaped according to the treatment. Environmental variables associated with changes in microbial community structure are plotted as vectors, where the length and direction indicate the degree of the correlation with the data.

Figure 5

Heatmap comparing microbial phylum abundance in control and manipulated samples during spring and fall. Darker shades indicate higher abundance. Phyla like Proteobacteriota and Actinobacteriota show variable distribution across samples.

Figure 5. Heatmap visualizes the relative abundance of microbial community in phyla level in different samples. Spring samples in the left and fall samples in the right panel. From left to right the samples are ordered by sampling depth. Treatments (control and manipulated) are marked with curved brackets. See Supplementary Table 1 for details on microbial samples.

3.3.2 Potential effects of the treatment on the microbial community

Although NMDS ordinations indicated that separation was primarily structured by season, SIMPER partitioning of Bray–Curtis dissimilarities showed that between-treatment differences were largely accounted for by a limited set of phyla. Across all samples, the top six phyla explained 72.8% of the manipulated–control dissimilarity (Supplementary Data Sheet 1): Acidobacteriota (+1.32 pp, 95% CI [−7.56, 9.86], p > 0.05), Verrucomicrobiota (+1.96 pp, [−1.08, 4.77], p > 0.05), Patescibacteria (+4.57 pp, [1.10, 8.44], p < 0.01), and Bacteroidota (+1.26 pp, [−1.81, 4.56], p > 0.05) were higher in manipulated plots, whereas Actinobacteriota (−5.80 pp, [−11.11, −0.41], p < 0.05) and Proteobacteria (−3.23 pp, [−8.36, 1.48], p > 0.05) were higher in controls.

Within-season analyses indicated that the treatment contrast was strongest in fall, where the top six phyla explained 71.9% of the manipulated–control dissimilarity (Supplementary Data Sheet 2). Actinobacteriota (−8.77 pp, [−15.14, −2.49], p < 0.05) and Proteobacteria (−3.64 pp, [−8.78, 1.63], p > 0.05) were higher in controls, whereas Acidobacteriota (+0.75 pp, [−10.24, 11.43], p > 0.05), Bacteroidota (+2.32 pp, [−1.76, 6.85], p > 0.05), Verrucomicrobiota (+3.28 pp, [0.21, 6.23], p > 0.05), and Patescibacteria (+7.84 pp, [3.25, 12.50], p < 0.01) were higher in manipulated plots. In spring, Acidobacteriota (+2.85 pp, [−9.96, 14.21], p > 0.05), WPS-2 (+1.09 pp, [−5.04, 7.56], p > 0.05), and Planctomycetota (+3.38 pp, [−0.43, 7.82], p > 0.05) tended to be higher in manipulated plots, whereas Proteobacteria (−2.42 pp, [−11.87, 6.34], p > 0.05) and Actinobacteriota (−0.92 pp, [−7.15, 6.20], p > 0.05) were higher in controls; the top six phyla explained 83.8% of the between-treatment dissimilarity (Supplementary Data Sheet 2).

At lower taxonomic ranks, within-season differential-abundance results were concordant with these phylum-level patterns. In spring, Subgroup 2 and Bryobacter (each ∼5% of the community) had higher relative abundances in manipulated plots than in controls (Figure 6A). In fall, Occallatibacter and Parcubacteria were higher in manipulated plots (each ∼10% of the community), and Pedosphaeraceae and Firmicutes were also elevated, albeit at lower levels (<10%) (Figure 6B). Collectively, these analyses indicate that while season explains the dominant community structure, treatment-associated differences are detectable, particularly in fall, and are attributable to shifts in a small number of phyla and their subordinate lineages.

Figure 6

Bar and scatter plots illustrate microbial abundance and LDA scores under control and manipulated treatments. Graph (A) shows g__Bryobacter and o__Subgroup2; (B) includes taxa like p__Firmicutes and g__Occallatibacter. Control is in blue and manipulated in orange.

Figure 6. Amplicon sequence variants (ASVs) show significant differences in relative abundance between control and manipulated plots in spring (A) and fall (B). The left panel shows relative abundances, and the right panel displays the effect sizes as Linear Discriminant Analysis (LDA) scores on a log₁₀ scale. Taxa were considered significantly different based on Kruskal–Wallis tests (p < 0.05), followed by pairwise Wilcoxon tests (p < 0.05), and an LDA score threshold >4.0. Taxa are ordered by descending LDA score. The closest classified taxonomic level is shown: g (genus), f (family), o (order), c (class), p (phylum).

3.3.3 Mercury-methylating taxa and seasonal dynamics

Mercury-methylating taxa were detected in both seasons, which included the confirmed methylator Geobacter and putative hgcAB-bearing lineages such as Desulfobacca, Holophagaceae, and CPR clades (Parcubacteria, Nomurabacteria, Woesebacteria, Kaiserbacteria). Their relative abundances were low to moderate and varied among cores and seasons (Supplementary Figure S10). Comparing seasons, the summed relative abundance of confirmed/putative methylators was higher in fall than in spring (difference in means: 3.37 percentage points (pp), permutation p < 0.01; difference in medians: 3.49 pp, p < 0.05). At the taxon level, this seasonal increase was most apparent for Desulfobacca and CPR lineages (notably Parcubacteria; in some cores also Nomurabacteria/Woesebacteria), whereas spring generally showed lower values (Supplementary Figure S10; Supplementary Table 4). Across treatments, results depended on season. Considering all samples together, manipulated>control for the mean (difference: 3.26 pp, p < 0.01) but not for the median (2.76 pp, p > 0.05). Within seasons, there was no treatment difference in spring, whereas in fall the snow-manipulated plots showed higher methylator abundance than controls (difference in means: 5.12 pp, p < 0.01; difference in medians: 4.24 pp, p < 0.001).

3.3.4 Correlations between MeHg and taxa associated with methylation-related processes

Pearson correlation analysis was performed to examine relationships between key microbial guilds, specifically those known to associate with mercury methylation, sulfate reduction, iron reduction and methanogenesis, and measured environmental variables, including MeHg concentrations and net methylation potential (%MeHg) (Supplementary Table 5). These functional groups are known to include microorganisms that either directly methylate mercury or create geochemical conditions favorable for methylation (Bravo and Cosio, 2020).

A larger number of such taxa were identified in the fall samples compared to the spring samples (Supplementary Figure S11). In the fall samples, MeHg concentrations in control plots were positively correlated with members of Methylacidiphilaceae and Methylomonadaceae, while %MeHg showed positive correlations with Desulfobacca, Geobacteraceae, Sideroxydans, Methylocystis, Methyloligellaceae, Methylotenera, and Methanobacterium (Supplementary Figure S12B, left). In manipulated fall plots, both MeHg and %MeHg showed stronger positive correlations with Geobacter and MM1 (a subgroup of Methylophilaceae) compared to controls, where these associations were weaker or absent (Supplementary Figure S12B, right). In manipulated plots in the spring samples, positive correlations with MeHg levels were observed only for members of Methylacidiphilaceae (Supplementary Figure S12A, right).

4 Discussion

4.1 Permafrost thaw after 17 years of increased snow accumulation

Since the field experiment was initiated in 2005, active layer thickening and surface water accumulation have been monitored to track permafrost thaw and the development of thermokarst ponds (Johansson et al., 2013). In 2022, when we conducted our sampling, a thicker active layer was observed in the manipulated plots, and most of the manipulated plots had ponds, while the majority of controls remained dry (Figure 2B). These observations are in line with earlier evaluation from the same mire (Dobiński, 2010) and align well with the effects observed in similar snow fence experiments across the Arctic (Lafrenière et al., 2013; O’Neill and Burn, 2015; D’Imperio et al., 2018).

Although the overall conditions aligned with expectations, some variability in thaw depth and surface water accumulation was observed. In the fall, when thaw depth reaches its maximum, the active layer was similarly thick in the dry control plots (thaw depth of 0.68 ± 0.17 m, Table 1) as previously recorded from these plots (Figure 2A) and other parts of the mire (∼50 cm; Dobiński, 2010; Johansson et al., 2011). These similarities support minimal permafrost thaw in five of six control plots since the start of the experiment, suggesting that they likely still represent current permafrost conditions at the location. In contrast, the submerged control plot showed a deeper thaw depth of 1.05 m suggesting plot-specific conditions such as ground subsidence and pond formation had occurred. This was a “wet” site from the outset of the manipulated experiment and thus more prone to form ponds. In field experiments like this, replicate plots can develop differently (e.g., Rodenhizer et al., 2023) given the natural site heterogeneity. For the six manipulated plots, no permafrost was recorded in the upper 1 m (Figure 2A; Table 1). Drilling next to the experimental site indicated a permafrost thickness of at least 15 m at Storflaket (Johansson et al., 2011). Permafrost is thus still assumed to be present in the manipulated plots, but at greater depths. The continuous increase in thaw depth in the manipulated plots (Figure 2A) demonstrates that artificially increased snow accumulation has continued to accelerate permafrost thaw since these plots were last described (Pascual and Johansson, 2022). Thermokarst ponds that developed during the first 7 years of manipulation (Johansson et al., 2013) only lasted a few weeks at first and had disappeared by late fall (Johansson et al., 2013), but more recently they formed as early as spring and persisted well into summer and fall (Olid et al., 2020; Pascual and Johansson, 2022). Altogether this demonstrates that the manipulated plots exhibit formation of thermokarst as a result of permafrost thaw, while the control plots still resembled the initial permafrost conditions.

4.2 Stability of total Hg stocks

This study is the first to investigate Hg dynamics in a field experiment simulating accelerated permafrost thaw by increased snow accumulation. Mean THg concentrations observed (57 ± 35 ng g⁻¹) are consistent with values from other studies of peatlands across the Arctic (Lim et al., 2020; Kirkwood, 2024) and R_HgC compares well to previously reported values from organic rich peatlands of northern Scandinavia (Tarbier et al., 2021). Low THg concentrations in the near-surface peat (0–5 cm) likely reflect a decline in atmospheric Hg levels (Sonke et al., 2023) or post-depositional processes that may alter peat Hg concentrations (Biester et al., 2007; Haynes et al., 2019). Elevated THg concentrations at 5–30 cm (Figures 3A,B) most likely reflect increased atmospheric Hg deposition since the industrial revolution, a trend that has been observed previously (Klaminder et al., 2008; Tarbier et al., 2021). The link between THg and %C in our data has been well documented in earlier studies (e.g., Grigal, 2003; Hararuk et al., 2013; Lim et al., 2020) and supports that soil organic matter, in addition to historic levels of Hg in the atmosphere, drives Hg distribution in peatlands.

We did not observe changes in Hg stocks in response to permafrost thaw, neither between treatments nor surface moisture conditions, when analyzing seasons combined or the spring data alone. However, in the fall, there was higher THg in the manipulated plots (p < 0.05, Supplementary Table S5). We argue that this likely stems from heterogeneity among collected cores, rather than a true treatment effect. Given the stability of Hg when bound to peat (Khwaja et al., 2006; Vasilevich et al., 2023) and the large amounts of Hg stored in the soil matrix (e.g., Lim et al., 2020) in comparison to expected inputs (Dastoor et al., 2022), treatment differences would not be expected for an isolated season. Furthermore, if differences between the two treatments developed during the summer period, we would have expected an accumulative effect over time and the difference between manipulated and controls to become detectable. We thus conclude that, in contrast to our original hypothesis, that the peat Hg stocks have not been altered by the accelerated permafrost thaw and thermokarst pond formation within the timeframe of this field experiment.

The stability of Hg stocks observed in our study is in contrast to earlier studies, where loss of Hg from the soil matrix upon permafrost thaw has been demonstrated along natural permafrost thaw gradients at a nearby mire (Klaminder et al., 2008; Fahnestock et al., 2019). A possible explanation could be that Hg was primarily lost from the surface peat, and that this loss was not sufficient to alter the stock of Hg in the full peat profile. However, the vertical profile of Hg concentrations (and R_HgC) in manipulated plots (Figure 3B) resembled the profile found in the control plots (Figure 3A), suggesting that Hg in the surface soils of manipulated plots has also remained stable. While the wet and anoxic conditions following thermokarst formation may have promoted Hg remobilization through desorption of particulate Hg and abiotic or microbial reduction of Hg^II to Hg⁰ (Peters et al., 2007; Fahnestock et al., 2019; Li et al., 2023), organic matter (OM) may have had a protective function. In the soil matrix, Hg is tightly bound to thiol groups (Skyllberg et al., 2000) and likely embedded within soil particles, limiting its availability for particle desorption and reduction (Baptista-Salazar et al., 2022). Olid et al. (2020) examined C cycling in the same field experiment (after 10 years) and found slower decomposition rates of young C in the upper 20 cm of manipulated plots. They also noted slower C sequestration due to surface vegetation shifts, while deeper peat showed enhanced decomposition as thaw progressed–resulting in similar net C loss across treatments, as also seen in other peatlands (Blume-Werry et al., 2019; Heffernan et al., 2020). It is possible that Hg has been affected in a similar manner: Hg levels in surface peat remained stable in thaw affected plots because the lower Hg sequestration rates were offset by slower remobilization of existing Hg in the peat resulting from the slower C decomposition rates. This is supported by no observed difference of R_HgC between treatments (Supplementary Table S5). We thus argue that although the Hg fluxes may have been impacted, this has not yet resulted in altered Hg stocks. The discrepancy between our findings and previous studies reporting Hg loss is likely due to heterogenic responses between permafrost peatlands upon thaw or differences in the age of the thermokarst system (which typically is not known in studies covering natural thaw gradients). Further investigations are needed to investigate Hg remobilization within the peat profile of different permafrost thaw stages and make informed assumptions about short- and long-term Hg remobilization potential.

4.3 Net production of methylmercury

Permafrost thaw did not influence MeHg stocks across the peat profiles, nor were surface MeHg concentration higher in thaw affected plots. These results deviates from several earlier studies of thaw gradients in northern Scandinavia reporting elevated MeHg concentrations or production (Klaminder et al., 2008; Fahnestock et al., 2019; Tarbier et al., 2021). MeHg concentrations (0.030–2.6 ng g⁻¹) and %MeHg (0.058%–14%) are within the same range as other Scandinavian sites (0.01–28 ng g⁻¹ MeHg, %MeHg: 0.02%–17%, Tarbier et al., 2021). We also did not observe a strong vertical trend of MeHg concentrations, contrasting with reports of elevated surface MeHg concentrations in response to thaw (e.g., Fahnestock et al., 2019; Tarbier et al., 2021). We attribute increased thaw rates over the past 17 years to the snow fences, which altered the ground thermal regime. Previous studies investigating Hg dynamics in natural thaw gradients, although comparable to our study in experimental approach, provide limited information on the time of thaw feature formation or rather the age of thermokarst. Signs of thermokarst, including ground subsidence and active layer thickening, were evident as early as 2012 (Johansson et al., 2013) and have continued to progress with variable responses to water table, soil wetness and pond formation. The gradual thaw in the form of ground subsidence at Storflaket mire, as opposed to more abrupt thaw, including, e.g., block erosion observed at other sites (e.g., in Tarbier et al., 2021), likely explains the absence of elevated MeHg levels in our thaw plots, as both abrupt thaw and longer thermokarst development history can enhance the methylation observed in other studies.

The degradation of peat OM with permafrost thaw leads to the release of DOC into thaw ponds (Abbott et al., 2014). Pondwater THg concentrations (4.4–23 ng L⁻¹) exceeded values of those reported from other thaw ponds (0.75–4.4 ng L⁻¹, MacMillan et al., 2015) and THg correlated positively with DOC (Supplementary Figure S7), suggesting DOC release from peat may facilitate THg mobilization to the dissolved phase. The dissolved THg did not, however, reflect the Hg content of the surface peat, indicating that peat Hg composition does not directly drive current pondwater composition (Supplementary Figure S7). There was no correlation between dissolved MeHg and DOC, nor between THg and MeHg concentrations, suggesting that the MeHg pool may be governed by factors other than DOC or THg in the thaw ponds. The negative correlation observed between DOC and %MeHg (Supplementary Figure S7) is similar to other boreal peatland porewaters (Lehnherr et al., 2012; Thompson et al., 2025), implying that MeHg production might be limited at higher DOM concentrations.

The large variation in MeHg concentration and %MeHg (Table 2) is similar to other Arctic wetland ponds (St. Louis et al., 2005; Lehnherr et al., 2012). As %MeHg can reflect long-term net MeHg production (Varty et al., 2021), recently developed ponds may have lower %MeHg, while perennial ponds exhibit much higher %MeHg. This could explain why the majority of the rather mature thaw ponds in our study have %MeHg values >1%, but the two control plot ponds measured in fall have very low %MeHg (C2 and C3, <1% MeHg, Table 2). The three exceptionally high %MeHg values of 27%–31% were observed in ponds that likely represent the sites with most advanced thaw (manipulated plots, largest ponds, active layer depth >1 m, Table 2). Similar observations have been made by Fahnestock et al. (2019) in a nearby permafrost peatland, where matured thaw sites (fens) featured increased MeHg production. Rather rapid abiotic processes (e.g., photo-demethylation) could be involved in the net production of MeHg and contribute to the variable dissolved MeHg concentrations. Additionally, the formation of refractory Hg pools could lead to reduced reactivity of MeHg in the water column as they control the availability of Hg to enter the dissolved pool (Baptista-Salazar et al., 2022) and thus limit its availability for further biogeochemical transformation. Other sources such as methylation in Sphagnum moss mats (e.g., in pond 6, Supplementary Figure S2B) or in littoral biofilms (Yu et al., 2010; Branfireun et al., 2020; Burke et al., 2020) could further influence the dissolved MeHg pool. Combined, these findings highlight the complexity of the Hg cycle in permafrost affected systems. Site-specific factors such as thaw dynamics, timeline of thaw and the resulting bioavailability of DOC released appear to play a bigger role in regulating the MeHg pool and distribution rather than the mere presence of thaw.

4.4 The interaction of putative methylators and MeHg

With respect to mercury cycling, manipulated plots with accelerated permafrost thaw exhibited microbial shifts relevant to methylation. In fall, thaw plots showed higher summed relative abundance of putative methylators than intact permafrost conditions of control plots, together with increases in Occallatibacter and Parcubacteria, and elevations of Pedosphaeraceae and Firmicutes. In spring, Subgroup 2 and Bryobacter were higher in manipulated plots. Across both seasons, communities contained the confirmed methylator Geobacter and low–moderate abundances of putative hgcAB-bearing lineages indicating potential for MeHg production.

Taxa that co-varied with treatments and seasons but are not themselves Hg methylators, for example, aerobic methanotrophs (Methylocapsa, Methylocystis, Methylomonas), methylotrophs (Methylophilaceae/MM1), the iron oxidizer Sideroxydans, and acidophilic methanotrophs (Methylacidiphilaceae), likely form microhabitats that favor Hg methylation, rather than catalyze hgcAB directly. Methanotrophs consume oxygen and generate low-molecular-weight organics that steepen oxic–anoxic microgradients and fuel anaerobes. More specifically, methylotrophs recycle C₁ substrates and iron oxidizers produce fresh Fe^III oxyhydroxides that sorb Hg/sulfide and subsequently serve as electron acceptors for Fe^III-reducing taxa (e.g., Geobacter), collectively creating conditions favorable to Hg methylation (Hanson and Hanson, 1996; Hou et al., 2008; Chistoserdova, 2011; Zhou et al., 2022; Bravo et al., 2018). Thus these associated taxa may indirectly increase Hg^II bioavailability and sustain anoxia even if they lack hgcAB or do not directly methylate Hg (Bravo et al., 2018; Fahnestock et al., 2019; Bravo and Cosio, 2020; Baptista-Salazar et al., 2022). The elevated abundance of putative methylators and anaerobic lineages in fall (e.g., Desulfobacca among sulfate reducers; Geobacter/Geobacteraceae among iron reducers) is consistent with wetter and more reduced conditions with deeper thaw depths. Despite these treatment-linked community shifts, bulk peat MeHg remained unchanged, indicating that methylator presence and abundance alone is insufficient to increase MeHg production in thaw plots. Geochemical controls, including redox state, electron-acceptor availability, DOM quality, Hg bioavailability and MeHg loss pathways (e.g., photodemethylation, sorption), likely limited realized methylation (Davidson and Janssens, 2006; Fahnestock et al., 2019; Jansson and Hofmockel, 2020; Baptista-Salazar et al., 2022; Thompson et al., 2025). In line with this, Desulfobacca correlated with THg and DOC, but not with MeHg or %MeHg in fall pondwater (Supplementary Figure S13), suggesting activity in organic-rich, sulfate-reducing niches that primarily supports C turnover while net MeHg formation remains unaffected. Overall, treatment effects were most evident in the fall, when deeper thaw and inundation amplified methylation potential (microbial capacity and habitat conditions), but did not translate into higher peat MeHg or THg stocks. As discussed above, this likely reflects (i) time lags between community change and geochemical pools, (ii) gradual (rather than abrupt) thermokarst development at Storflaket (Johansson et al., 2013; Olid et al., 2020), and (iii) strong sorptive control of Hg in peat (Skyllberg et al., 2000; Baptista-Salazar et al., 2022). Continued monitoring that combines community data (including hgcAB screening), MeHg and Hg isotope measurements, and redox and DOM metrics will be will be critical to determine whether ongoing thaw eventually shifts the balance from potential to realized MeHg production at Storflaket mire (Fahnestock et al., 2019; Tarbier et al., 2021; Wang et al., 2023).

5 Conclusion

Our results demonstrate that accelerated permafrost thaw over 17 years had no net effect on the net Hg and MeHg pools. Although permafrost thaw has clearly progressed through active layer thickening and thermokarst pond formation, these physical changes have not translated into greater Hg mobilization post thaw. Mercury appears to remain largely bound within the peat matrix, likely stabilized by strong association with OM and limited exchange between solid and dissolved phases. Microbial data revealed noticeable seasonal variations of putative Hg-methylators in thaw-affected plots, yet these shifts did not result in higher MeHg concentration. Nevertheless, given that changes in this system are slow and marked by gradual rather than abrupt thaw, peat Hg levels may need several decades to be altered and bulk MeHg levels could potentially follow the rise of putative Hg methylators. The extent of these future changes in the Hg cycle will likely depend on the rate and extent of thaw and associated hydrological responses. These findings are highly relevant from the perspective of ongoing climate warming in regions that are already vulnerable to increased permafrost thaw. Long-term, process-based studies are needed to address how progressive permafrost degradation will influence Hg mobilization from the peat matrix, ecosystem exposure and potential health risks in Arctic regions.

Data availability statement

The sequencing data are publicly available through the National Center for Biotechnology Information (NCBI) database under BioProject PRJNA1303211, BioSample accession SAMN50501911.

Author contributions

CH: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review and editing. AA: Conceptualization, Writing – review and editing, Investigation, Methodology. SB: Methodology, Writing – review and editing, Resources, Supervision. MJ: Methodology, Supervision, Writing – review and editing. TL: Data curation, Investigation, Visualization, Writing – original draft, Writing – review and editing. LT: Writing – review and editing, Investigation. SJ: Writing – review and editing, Conceptualization, Funding acquisition, Resources, Supervision.

Funding

The author(s) declared that financial support was received for this work and/or its publication. CH and SJ have received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 860497.

Acknowledgements

We thank the Swedish Polar Research Secretariat and SITES for the support of the work done at the Abisko Scientific Research Station. SITES is supported by the Swedish Research Council. We thank the staff at the Abisko Scientific Research Station and Emily Pickering Pedersen for infrastructure and logistical support during fieldwork. We are grateful to Anna Hägglund for laboratory support with THg analyses, Arild Gustafsson for elemental analysis, Mingyue Li for laboratory help and Mitsy Polderman and Prune Leroy for microbial 16S rRNA sequencing.

Conflict of interest

Author LT was employed by Hatfield Consultants.

The remaining author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvc.2025.1715512/full#supplementary-material

References

Abbott, B. W., Larouche, J. R., Jones, J. B., Bowden, W. B., and Balser, A. W. (2014). Elevated dissolved organic carbon biodegradability from thawing and collapsing permafrost. J. Geophys. Res. Biogeosciences 119 (10), 2049–2063. doi:10.1002/2014JG002678