Methane cycling microbes are important predictors of methylmercury accumulation in rice paddies

TL;DR

Methane-cycling microbes in rice paddies are key to predicting methylmercury levels, which pose health risks, and could help reduce both pollution and greenhouse gases.

Contribution

Methanogen and methanotroph abundances are shown to be better predictors of methylmercury than mercury methylation genes.

Findings

Methane cycling genes and microbial beta-diversity are strongly associated with methylmercury levels.

Methanogen abundance correlates with higher methylmercury under high total mercury.

Methanotroph mbnT gene abundance unexpectedly correlates positively with methylmercury.

Abstract

Microbial production of methylmercury from inorganic mercury in rice paddies poses health risks to consumers of this essential dietary staple. Although mercury-methylating communities are well characterized, the microbial guilds contributing to methylmercury accumulation in rice paddies remain unclear. Here, we collected paddy soils across a mercury concentration gradient throughout the rice-growing season to identify microbial and environmental factors influencing methylmercury dynamics. We show that hgcA gene abundance, the key gene required for methylation, was not a significant predictor of methylmercury concentration in paddy soils. We also show that the merB gene abundance correlated with methylmercury in mercury-polluted rhizosphere samples. Methane cycling genes were actively expressed, and their beta-diversity was significantly associated with methylmercury levels. Methanogen abundance correlated with higher methylmercury under elevated total mercury concentrations. Analysis of the methanotroph-associated mbnT gene, implicated in demethylation, revealed an unexpected positive correlation with methylmercury. Multiple regression and machine learning models converged on mercury bioavailability and methanogen/methanotroph abundances as key predictors of methylmercury, with methanogen-associated hgcA gene abundance and methanogen-methanotroph interactions highlighted under flooded, low-redox conditions. These findings suggest that methane-cycling microbes play key roles in methylmercury cycling dynamics and point to management strategies that could simultaneously mitigate mercury pollution and greenhouse gas emissions.

Methylmercury is a microbially derived neurotoxin that accumulates in the food staple rice (Oryza sativa). Mitigating the health effects of methylmercury exposure requires predicting mercury cycling dynamics in rice paddies. This task is challenging because of the complex interplay of microbial and environmental factors. Our study coupled genomic and geochemical measurements with machine learning models to identify the key biological indicators of methylmercury accumulation. We demonstrated that the abundance of methanogens and methanotrophs is a major microbial predictor of methylmercury variability. This predictive framework, which considers the interactions between these coupled microbial guilds, offers greater power than methods relying only on mercury methylation genes. These findings inform better management and remediation strategies for rice paddies, offering a path to reduce methylmercury exposure and mitigate greenhouse gas emissions.

Figures

Hg dynamics in sampled rice paddies over the growing season. (A) Total mercury (THg), (B) MeHg, and (C) water-soluble Hg fraction (F1) in surface soils (left panels) and rhizosphere soils (right panels). Data points represent individual measurements from replicate soil samples collected at three sites (HX, GX, and SK), with smoothed trend lines (±95% CI). Note that surface soils were sampled throughout the entire study period (including post-harvest), while rhizosphere soils were only collected during the rice-growing season when plants were present, resulting in different x-axis scales. THg is shown in mg/kg, while MeHg and F1-fraction are in µg/kg. Day 0 corresponds to the beginning of the sampling period.

Analysis of hgcA genes and their correlation with MeHg. (A) Temporal dynamics of hgcA gene abundance across four major metabolic guilds: IRB (blue), SRB (green), methanogens (orange), and syntrophs (purple). Unclassified hgcA sequences were excluded from visualization but retained in total abundance calculations. Taxon-level hgcA abundances (family level) for all samples are provided in Data S7. (B) Spearman correlation coefficients between MeHg concentrations and gene abundances for different metabolic guilds. Syntrophs are not shown due to insufficient data. Asterisks indicate statistical significance levels (*P < 0.05; ns: not significant).

Community structure and environmental drivers of methanogens and aerobic methanotrophs. (A) Taxonomic distribution and relative abundance of mcrA genes across sampling days. (B) Taxonomic distribution and relative abundance of pmoA/mmoX genes across sampling days. (C and D) db-RDA of beta-diversity in relation to geochemical variables, with panel C showing methanogens and panel D showing methanotrophs. Points represent individual samples colored by site. Arrows indicate environmental variables with significance levels (*P < 0.001, P < 0.01, and *P < 0.05).

Presence/absence of methanobactin-related genes in methanotroph genomes. This heatmap illustrates the presence and absence of methanobactin biosynthesis and transport genes across experimentally studied methanotrophs (24) and MAGs recovered in this study, along with their respective quality metrics. The x-axis denotes the methanobactin genes, while the left y-axis represents the taxonomic classification of each genome, and the right y-axis lists the RefSeq or MAG IDs. All RefSeq genomes lack methanobactin biosynthesis, except OB3B (Note mbnB is present in OB3B as a pseudogene, which is not recovered through the HMM-based analyses, but confirmed with RefSeq annotation).

Variable importance and model performance for predicting MeHg concentrations using RF analysis (full model). (A) Variable importance plot showing the mean percent increase in mean squared error (%IncMSE) for each predictor variable when permuted, with 95% confidence intervals across 100 different data splits. Larger point size and green coloration indicate higher significance rates. Values show the percentage of data splits where the variable was significant (P < 0.05), followed by the mean P-value. (B) Model performance metrics displaying mean values with 95% confidence intervals across all data splits. The model was built using 2,001 trees with mtry = 4 and validated through 100 permutation tests with random 70/30 train/test splits.

Equations

Full text

INTRODUCTION

Mercury (Hg) is a metal primarily emitted as gaseous elemental Hg (Hg[0]) from natural and anthropogenic sources and deposits in ecosystems worldwide, where it undergoes oxidation to Hg(II) and subsequent microbial conversion to methylmercury (MeHg), a potent neurotoxin (1). MeHg accumulates in rice, a global food staple, with particularly elevated concentrations in paddies near mining, smelting, and industrial areas with elevated Hg(0) emissions (2). The global atmospheric transportation of Hg could also affect rice paddies further away from pollution sources (3). Since paddy soil represents a major source of MeHg in rice grains (4), predicting MeHg accumulation in rice paddies is critical for mitigating Hg exposure through diet.

Net MeHg accumulation in rice paddies depends on the interplay between methylation and demethylation processes (5). Hg methylation is facilitated by anaerobic microbes harboring the hgcAB gene cluster (6–8), while demethylation occurs via the Hg resistance (mer) operon-dependent pathways (e.g., merB/merA) and mer-independent microbial pathways (5). These processes are also influenced by environmental variables such as Hg bioavailability, redox potential, dissolved organic matter (DOM), sulfide (5, 9–12), and other microbial players indirectly (13).

While Hg methylating communities have been extensively characterized in rice paddies and other ecosystems by surveying the genetic determinants of Hg cycling (14–19), MeHg poorly correlates with Hg-cycling genes (e.g., hgcA, merA, and merB) and hgcA transcripts (20). This disconnect highlights the need for integrated approaches that consider multiple variables simultaneously when predicting what controls MeHg accumulation in rice.

Recent evidence suggests that methane-cycling microbes play underappreciated roles in mediating MeHg dynamics. In rice paddies near Hg-mining areas, methanogens emerged as critical players controlling MeHg levels (14, 15, 21), potentially contributing to both methylation and demethylation simultaneously (14, 22, 23). Furthermore, methanotrophs were experimentally shown to degrade MeHg at picomolar to nanomolar concentrations under near-neutral pH conditions typical of rice paddies, using a methanobactin-mediated mechanism distinct from the conventional mer-dependent pathway (24, 25). However, the ecosystem-level significance of methanotrophic MeHg demethylation and its interaction with methanogenic processes in determining net MeHg accumulation has not been evaluated in rice paddy systems.

To address these gaps, we aim to (i) identify and rank the relative importance of microbial communities and environmental variables relevant to MeHg turnover as predictors of MeHg accumulation patterns in mining-impacted rice paddies; (ii) characterize the abundance and distribution of methane-cycling microorganisms in relation to MeHg concentrations; and (iii) develop predictive models to inform future mechanistic studies on Hg cycling.

By integrating metagenomic, geochemical, and statistical approaches, we show that methanogen and methanotroph abundances are leading biological predictors of MeHg content, while Hg bioavailability is the primary environmental determinant of MeHg variability. Machine learning analyses further highlighted the importance of methanogen methylator abundance and a methanogen–methanotroph interaction, particularly under flooded conditions, where model performance improved. Together, these findings provide environmental-scale data that can guide hypothesis development for broader geographical testing and inform strategies to jointly mitigate greenhouse gas emissions and Hg exposure.

MATERIALS AND METHODS

Soil collection, geochemical analyses, and sequencing

Rice paddy soil samples were collected from three sites in Guizhou province, Southwest China: Huaxi (HX, background site), Gouxi (GX, artisanal Hg smelting area), and Sikeng (SK, abandoned Hg mining area). Sampling occurred from May 2021 to April 2022, spanning pre-plantation through post-harvest periods. Surface and rhizosphere soils were collected and processed as follows: Individual replicate samples were analyzed separately for geochemical parameters, while replicate samples were combined into composite samples for DNA extraction and subsequent metagenomic analyses. Soil geochemical parameters were analyzed, including Hg species (THg by Cold Vapor Atomic Fluorescence Spectrometry [CVAFS], MeHg by gas chromatography-CVAFS followingUSEPA method 1630, and the bioavailable fraction Hg by the Mg[NO₃]₂ extraction method), anions (by ion chromatography), and elemental sulfur (by high-performance liquid chromatography). Dissolved organic carbon and soil organic matter were analyzed from water-extractable constituents from the soil. Redox-sensitive species (S^2−^, Fe²^+^, and Fe³^+^) were analyzed from porewater collected during the flooded period. Water-soluble and exchangeable Hg fractions were measured following an established protocol (26) based on earlier works (27, 28). Briefly, freeze-dried soil samples (1 g) were placed in 50 mL polytetrafluoroethylene tubes and extracted with 8 mL of 1 M Mg(NO₃)₂ solution (pH 7). Following agitation at 150 rpm for 1 h, samples were centrifuged at 3,000 rpm for 15 min. The resulting supernatant was filtered through 0.45 μm nitrocellulose filters and transferred to clean 25 mL borosilicate containers. Hg concentrations in the filtered extracts were determined using the same analytical protocol applied for total Hg analysis in aqueous samples. Recovery rates of the standard reference materials for THg and MeHg analyses are provided in Data S1. DNA and RNA extractions were performed using the PowerSoil DNA Isolation Kit and the RNeasy PowerSoil Total RNA Kit (QIAGEN) following the manufacturer’s instructions. All geochemical analyses and nucleic acid extractions were performed at the Institute of Geochemistry, Chinese Academy of Sciences (Guiyang, China). Both metagenomic and metatranscriptomic sequencing were performed on an Illumina NovaSeq 6000 platform, producing 150 bp paired-end reads. Sequencing was conducted at Azenta (Suzhou, China). The sequencing aimed for a read depth of approximately 15 Gbp per sample (range: 10.9–20 Gbp). Quality metrics are provided in Data S2. In total, 50 metagenomic and 6 metatranscriptomic samples were obtained. Details on sample collection, geochemical analyses, and sequencing are provided in Method S1.

Genome assembly and binning

Metagenomic data were processed using fastp and FastQC for quality control (Data S3) (29, 30). Individual assemblies were performed with MetaSPAdes, while site-specific co-assemblies used MEGAHIT (31). Contigs >2,000 bp were retained and mapped with BWA-MEM to generate coverage data. Key assembly metrics, such as N50, total contigs, length distribution, and mapping rates, are provided in Data S4. Four binning algorithms (MetaBAT 2, MaxBin 2, CONCOCT, and VAMB) were applied (32–35), yielding 5,307 initial bins. Bin quality was assessed with CheckM2 (36), and DAS-Tool (37) was used to select the highest quality bins, followed by dereplication with dRep (38) (ANI threshold 98%). We retained 267 metagenome-assembled genomes (MAGs) with >50% completeness and <10% redundancy following manual refinement in Anvi'o (39). Taxonomic classification was performed using GTDB-Tk (40) based on GTDB (release 214; Data S5). The co-assembly approach was chosen over individual assembly due to higher MAG yield and greater taxonomic diversity. Details on genome assembly, binning, and software versions are provided in Method S2.

Metatranscriptomic sequence processing

The RNA sequences were trimmed and filtered to a minimum mean quality score of 30 using fastp (29). Forward and reverse sequences were overlapped with a minimum and maximum of 10 and 150 bp, respectively, using FLASH (41). Sequences were subsequently sorted into SSU ribosomal RNA (rRNA), LSU rRNA, and non-rRNA using SortMeRNA (42) against the default reference database (smr_v4.3_default_db). The non-rRNA fraction was mapped to the open reading frames (ORFs) identified on the contigs from the corresponding sample and linked with featureCounts results to determine the transcript abundance of mcrA, pmoA, mmoX, hgcA, and merB. Taxon-specific transcript abundance was calculated by dividing the RNA reads mapped to transcripts by the total RNA reads in the sample, normalized to hidden Markov model (HMM) length.

Functional gene abundance analyses

For microbial abundance estimation of various functional guilds, including methanogens, methanotrophs, Hg-methylating microbes, and demethylating microbes, we employed a gene-centric approach at the contig level using individual assemblies. We selected individual assemblies over co-assembly to preserve sample-specific information and enable more accurate quantification of functional genes within each sample. While co-assembly is valuable for generating high-quality MAGs by aggregating data across samples, it inherently homogenizes information from multiple samples, potentially obscuring subtle variations in gene abundance at the individual sample level.

ORFs were identified using Prodigal (43), short reads were mapped using BWA-MEM (44), and reads mapped to ORFs were quantified using featureCounts from Subread package (45). Functional genes were retrieved using HMMs through hmmsearch of HMMER software (46). Specifically, the McrA HMM (PF02249) was sourced from PFAM, while PmoA and MmoX HMMs were obtained from the GraftM gene packages (https://data.ace.uq.edu.au/public/graftm/7) with an E-value threshold of 1 × 10^-5^. We took an additional filtering step to maximize the correct recovery of methane-cycling gene sequences. We conducted protein homology searches using BLASTp against the NCBI non-redundant database (as of 5 December 2023) to confirm their identity and refine search results. Queries consisted of HMM hits of the corresponding genes, each limited to the single best hit (-max_target_seqs 1) to focus on the most relevant homolog. McrA, PmoA, and MmoX hits were filtered based on sequence titles indicative of their respective functions (e.g., McrA hits containing “methyl coenzyme M reductase”; PmoA and MmoX hits classified as “methane monooxygenase”). The PmoA sequence is homologous to other monooxygenases, like AmoA of ammonia oxidizers, so we only retained sequences classified to known methanotrophic taxa to ensure accuracy.

Methanobactin biosynthesis and transport genes from reference genomes, MAGs, and metagenomes were identified using HMMs listed from a previous study with predefined threshold values (47).

Hg-cycling genes required specialized validation procedures. hgcAB genes were identified using Hg-MATE-Db HMMs (ver. 1.01142021) (48) and filtered for conserved motifs [N(V/I)WCA(A/G)GK] and [C(M/I)ECGA] (6). MerB sequences were retrieved using pre-compiled HMMs (49) and validated for characteristic cysteine and aspartic acid signatures (50).

For additional functional guilds, detection of iron-reducing bacteria (IRB) employed HMM profiles from FeGenie (51) to identify genes related to iron reduction, with gene-specific domain E-value cutoffs as defined in the database configuration files. Abundance for IRB was calculated by summing coverage counts of all identified genes and normalizing against total coding sequence counts in each metagenome. For quantification of sulfate-reducing bacteria (SRB), we targeted the dsrA gene using the TIGR02064.1 HMM profile with a sequence cutoff of 392.9.

Taxonomic classification was conducted using MMseqs2 (52) against GTDB (release 214) (53) with lowest common ancestor strategy. To facilitate cross-sample comparison, we converted taxon-specific coverage values to the percentage of reads normalized to HMM lengths as shown in the equation below. Total reads numbers in each sample were referred to as effective library size values calculated using the edgeR package (54). Details on functional gene analyses and software versions are provided in Method S3.

$$[eqn]$$

Complementary metagenomic approaches

Our study employed two complementary metagenomic strategies. Gene-centric analysis quantifies functional potential across all community members, capturing the full breadth of metabolic capacities, including genes from low-abundance taxa that may not assemble into complete genomes. Genome-resolved analysis using MAGs provides taxonomic context, linking functional genes to specific organisms and enabling network analyses of population-level responses to environmental changes. Together, these approaches address different scales: community-wide functional potential (gene-centric) and organism-specific contributions (genome-resolved), both essential for understanding Hg cycling processes.

Statistical analyses

Principal component analysis of geochemical variables was performed with standardized data, and distance-based redundancy analysis (db-RDA) was used to assess relationships between methane-cycling communities and environmental factors based on Bray-Curtis dissimilarity. Permutational multivariate analysis of variance (PERMANOVA) evaluated the effects of site and soil compartment on community composition (i.e., beta-diversity), with pairwise comparisons using Bonferroni correction. Multiple linear regression with Akaike information criterion corrected (AICc) model selection was employed to evaluate the relative importance of methylation potential (hgcA abundance) and Hg bioavailability (F1-Hg) in predicting MeHg concentrations. Seven candidate models representing different hypotheses were tested, and model assumptions were verified through standard diagnostic tests.

MAG relative abundances were calculated by aligning reads to dereplicated genome sets using CoverM (v0.7.0) (55). Time-dependent correlations between MAG abundance and MeHg concentrations were identified using extended Local Similarity Analysis across 11 time points for each site (Data S6) (56, 57). False discovery rate correction was applied using the Benjamini-Hochberg procedure, with significant associations defined as FDR-adjusted P < 0.05. Random forest (RF) regression models were developed to predict MeHg concentrations from 12 predictor variables (nine microbial functional guilds and three environmental parameters) using 100-fold cross-validation with 70:30 train-test splits. Model performance was evaluated using root mean square error (RMSE), mean absolute error (MAE), R², and Pearson correlation coefficients. All analyses were conducted in R (v4.3.0) with statistical significance set at α = 0.05. Detailed methodological parameters are provided in Method S4.

RESULTS AND DISCUSSION

Spatial and temporal dynamics of MeHg concentration across the rice paddies

THg and MeHg levels varied significantly across rice paddies (Fig. 1; Fig. S1). Average THg concentrations across the entire sampling period were significantly elevated at the mining-impacted sites SK (42.03 ± 21.5 mg/kg) and GX (21.41 ± 4.04 mg/kg), compared to the background site HX (0.17 ± 0.02 mg/kg). However, average MeHg concentrations over the sampling period were significantly higher at GX (3.44 ± 0.70 µg/kg), followed by SK (2.42 ± 0.76 µg/kg) and HX (0.78 ± 0.21 µg/kg; Fig. S1). This supports prior studies indicating that THg is a poor predictor of the net MeHg production (58).

Hg dynamics in sampled rice paddies over the growing season. (A) Total mercury (THg), (B) MeHg, and (C) water-soluble Hg fraction (F1) in surface soils (left panels) and rhizosphere soils (right panels). Data points represent individual measurements from replicate soil samples collected at three sites (HX, GX, and SK), with smoothed trend lines (±95% CI). Note that surface soils were sampled throughout the entire study period (including post-harvest), while rhizosphere soils were only collected during the rice-growing season when plants were present, resulting in different x-axis scales. THg is shown in mg/kg, while MeHg and F1-fraction are in µg/kg. Day 0 corresponds to the beginning of the sampling period.

Distinct temporal trends were observed for MeHg at the mining-impacted sites, which exhibited different patterns compared to the relatively stable levels of THg (Fig. 1A). At GX and SK, MeHg peaked at cultivation onset (day 0) and declined throughout the rice-growing season (day 0–80), with further reductions observed post-harvest after drainage. MeHg concentrations were typically higher during the water-saturated season (day 0–100) than in the post-harvest period (day 100–260). A notable increase in MeHg occurred at GX on day 290, coinciding with pre-planting and re-drainage for the next growing season. In contrast, at HX, MeHg levels remained relatively stable overall. Throughout the three sampling sites, no significant differences in MeHg were detected between surface soil and rhizosphere across the three sampling sites (Mann-Whitney U test, P > 0.05). This pattern of fluctuating MeHg in mining-affected sites during cultivation and stability at the background site is previously undocumented, due to the lack of extensive longitudinal sampling. Subsequent sections consider the potential microbial and geochemical drivers contributing to these trends in Hg concentrations.

Geochemical factors partially explain site-specific MeHg patterns

Geochemical factors were assessed to contextualize their potential impacts on MeHg accumulation (please refer to Fig. S2; Result S1 throughout this section). Redox potential measured in overlying water revealed that HX is more oxidizing (109.29 ± 19.12 mV) compared to GX (57.41 ± 37.09 mV) and SK (60.96 ± 20.3 mV). Regardless of soil compartments, DOM concentrations differed significantly between sites, with HX showing the highest levels (0.72 ± 0.31 mg/g). MeHg correlated positively with DOM in surface soil at GX (Spearman’s ρ = 0.33) and SK (ρ = 0.58), while correlating negatively with DOM SUVA_254_ in SK surface (ρ = −0.56) and rhizosphere (ρ = −0.60; Fig. S3). The bioavailable fraction of Hg (F1-Hg) exhibited temporal fluctuations that partially paralleled MeHg changes, primarily in GX surface soil (ρ = 0.54, P = 0.0034, n = 27) and SK rhizosphere (ρ = 0.62, P = 0.014, n = 15; Fig. 1C; Fig. S3). Noticeably, GX initially showed higher F1-Hg than SK, but these differences diminished over the growing season with both sites reaching comparable levels after day 100 (Fig. 1C). This decline aligns with the post-harvest period when rice fields are dried and aerated, promoting immobilization of dissolved Hg species through redox-driven precipitation and adsorption processes (e.g., Fe-oxyhydroxide precipitation and adsorption) (12, 59).

Several key factors showed no significant differences between sites. Soil pH displayed similar temporal trends at HX and GX, while SK showed more apparent fluctuation but with no significant site differences overall (Kruskal-Wallis test, P > 0.05, Fig. S1). Terminal electron acceptors, including sulfate (avg. ~100 mg/kg or 100 μM) and nitrate (avg. ~4 mg/kg or 6.45 μM), demonstrated similar temporal patterns across all sites. Hydrogen sulfide remained consistently low (0–5 μM) with no significant site differences (Kruskal-Wallis test, P > 0.05, Fig. S1). This suggests similar rates of microbial sulfate reduction and constraints on Hg methylation or bioavailability, as sulfide can readily react with Hg(II), forming HgS (60, 61). Fe²^+^/total Fe ratios showed similar temporal trends except in HX subsurface porewater, which displayed lower ratios, suggesting reduced iron reduction at the depth of rice roots.

The three sites differ primarily in bioavailable Hg content (i.e., F1-Hg) and DOM, which are both variables known to control Hg methylation (62–64). Since DOM affects Hg(II) bioavailability at low sulfide conditions (≲30 μM) (65, 66), the consistently low sulfide levels and differences in F1-Hg suggest that DOM potentially contributes to the MeHg fluctuations observed. Interestingly, DOM aromaticity (i.e., SUVA_254_) correlated negatively with MeHg at SK, where lower DOM SUVA_254_ (indicating higher DOM bioavailability [67]) correlated with higher MeHg in both soil compartments.

Taken together, our geochemical data suggest that iron and sulfur cycling processes have limited influence on net MeHg accumulation and are insufficient to explain the observed spatiotemporal patterns of MeHg fluctuation across our study sites. While oxidation-reduction potential (ORP) measurements (Fig. S1) indicate that HX conditions favor microaerophilic to anaerobic processes (denitrification and manganese reduction), and GX and SK are more conducive to iron reduction and sulfate reduction, these differences in redox environments do not adequately account for the MeHg distribution patterns observed. Rather, our findings suggest that net MeHg accumulation is likely controlled by methylating microbes beyond SRB and IRB and active demethylation processes. These results highlight the need to assess the potential contributions of alternative microbial methylators and demethylators to MeHg dynamics in rice paddy environments under the geochemical conditions observed in this study.

MeHg levels show limited association with hgcA abundance across metabolic guilds

As the key marker gene of Hg methylation, we searched for and classified hgcA in the assembled metagenomes (Fig. 2) to assess methylation potential and identify the microbial players involved. Total hgcA abundance was, on average, higher at mining-impacted sites GX (mean = 8.4 × 10^−8^ ± 3.0 × 10^−8^; n = 16) and SK (mean = 8.4 × 10^−8^ ± 2.7 × 10^−8^; n = 17) compared to the background site (mean = 5.5 × 10^−8^ ± 1.8 × 10^−8^; n = 16), aligning with the higher MeHg levels observed at contaminated sites ([Fig. 1 and 2](#F1 F2)). Metatranscriptomic analysis confirmed active hgcA expression across all sites (n = 6), with similar expression in surface soils (∼1.0 × 10^−8^) but substantially higher expression in GX and SK rhizospheres (∼1.4–1.7 × 10^−8^), supporting the functional relevance of the metabolic potential characterized using metagenomics (Fig. S4A).

*Analysis of hgcA genes and their correlation with MeHg. (A) Temporal dynamics of hgcA gene abundance across four major metabolic guilds: IRB (blue), SRB (green), methanogens (orange), and syntrophs (purple). Unclassified hgcA sequences were excluded from visualization but retained in total abundance calculations. Taxon-level hgcA abundances (family level) for all samples are provided in Data S7. (B) Spearman correlation coefficients between MeHg concentrations and gene abundances for different metabolic guilds. Syntrophs are not shown due to insufficient data. Asterisks indicate statistical significance levels (P < 0.05; ns: not significant).

Temporally, hgcA abundance at HX remained relatively stable, mirroring MeHg trends (Fig. 1B). Across the growing season, correlations between hgcA abundance and soil MeHg were highly inconsistent (Fig. 2B). A significant positive correlation between MeHg and hgcA abundance was restricted solely to the SK rhizosphere, observed for both total hgcA and IRB-associated hgcA abundance, while correlations for other sites and guilds were non-significant. Taxonomically, hgcA was primarily assigned to Deferrimicrobiaceae at HX, whereas Methanoregulaceae-related hgcA increased at GX and SK, particularly in rhizosphere samples (Data S7).

As an additional test of drivers of MeHg accumulation, we fitted multiple linear regression models with model selection by AICc (Result S3). The best-supported model included additive effects of bioavailable Hg(II) (F1-Hg) and overall hgcA, while controlling for site and soil compartment. In this model, F1-Hg was a major positive predictor of MeHg concentrations (β = 0.73, P < 0.001), whereas total hgcA abundance was not significant. Adding an interaction between F1-Hg and hgcA did not improve fit (ΔAICc = 3). Consistent with this, the explained variance (adjusted R^2^) increased from ~60% (additive model without site/compartment) to ~84% when site and compartment effects were included, underscoring the influence of other environmental on MeHg.

These results suggest that net MeHg accumulation is not primarily controlled by methylation potential in Hg-contaminated rice paddies, but rather by Hg substrate availability and other environmental factors. This contrasts with findings from diverse aquatic systems (e.g., lakes, reservoirs, and brackish systems) where direct correlations between hgcA abundance and MeHg production have been observed under more uniform geochemical conditions (68–70). In other cases, studies have demonstrated that MeHg formation is governed by a synergy of Hg bioavailability and methylation capacity, suggesting that either factor can act as a rate-limiting step (71, 72). Specifically, microbial productivity can limit MeHg production regardless of bioavailability below a certain threshold (71), while community shifts caused by environmental constraints may outweigh subtle changes in substrate availability (69).

In rice paddies, the relationship between hgcA and MeHg content is complicated by several factors: (i) pronounced geochemical gradients in DOM characteristics, redox conditions, and Hg bioavailability; (ii) the presence of multiple MeHg sinks including biotic and abiotic demethylation processes (20, 73); and (iii) MeHg uptake by rice plants, which demonstrate efficient MeHg bioaccumulation from paddy soil (i.e., bioaccumulation factors of 5.5–5.6) (74, 75). Although the quantitative effect of plant uptake on soil MeHg concentrations was not considered in this study, plant-mediated removal may partially contribute to the soil MeHg content we observed. The interplay between MeHg production, degradation, and plant uptake processes operating simultaneously in rice paddies underscores the need to consider the genetic potential of methylating microbes alongside broader ecosystem variables that operate at larger spatial scales to better predict the fate of MeHg in rice paddies.

merB-mediated demethylation and MeHg dynamics

Following hgcA analyses, we investigated demethylation potential represented by merB, which encodes the organomercurial lyase that performs reductive demethylation (76, 77). Total merB gene abundance in the rhizosphere was generally higher compared to surface soil, particularly at HX and SK (HX: 5.1 × 10^−8^ ± 5.5 × 10^−8^ surface vs. 8.2 × 10^−8^ ± 5.6 × 10^−8^ rhizosphere; SK: 5.8 × 10^−8^ ± 3.9 × 10^−8^ surface vs. 1.6 × 10^−7^ ± 2.2 × 10^−7^ rhizosphere), though GX showed similar abundances between compartments (3.5 × 10^−8^ ± 5.2 × 10^−8^ surface vs 3.9 × 10^−8^ ± 3.5 × 10^−8^ rhizosphere; Fig. S5). The higher merB abundance in rhizosphere samples is likely attributed to radial oxygen loss (ROL) (78), creating more oxidizing conditions conducive to aerobic *merB-*carrying microbes (79), although direct soil ORP measurements were not obtained in this study. The highest rhizosphere merB was observed at day 0, followed by a gradual decline, aligning with previous findings that ROL rates decrease over the rice growth period (79). This temporal trend further supports the potential role of plant-mediated oxygen dynamics in structuring the merB-carrying community.

Multiple regression analysis revealed that merB gene abundance has a significant positive association with MeHg levels (β = 2.51 × 10^6^, P = 0.0113) after accounting for site-specific effects. The positive association between gene abundance and absolute MeHg concentrations suggests adaptation of microbial communities to Hg-contaminated environments. In highly polluted aquatic systems, increased THg load drives higher absolute MeHg concentrations (76). The microbial communities respond by enriching for Hg-resistant bacteria carrying the mer-operon, including merB genes (76).

The correlation between merB and MeHg is particularly strong in GX (Spearman’s ρ = 0.6; n = 5) and SK (Spearman’s ρ = 0.83; n = 6) rhizosphere samples. This suggests that the ROL-driven oxidizing conditions select for the aerobic merB-carrying organisms, and where THg loads are high (GX and SK), these adapted populations are most enriched, leading to the strong co-variance between merB gene abundance and MeHg concentration in this specific micro-environment.

At the transcriptomic level, merB transcripts were detected only in HX surface soils (n = 6, Fig. S4B). This disconnect from gene abundance likely reflects: (i) merB-mediated demethylation is not active as MerB was favorably expressed at oxic and high-Hg environments; (ii) the transient nature of merB transcription and the short half-lives of mer transcripts, making them easy to miss in single time-point sampling; (iii) temporal fluctuations in Hg(II) bioavailability that narrow the window for detectable merB expression (76, 80). These findings suggest a lack of MerB activity in this system and highlight the need for temporal transcriptomic sampling to establish relationships between merB expression and MeHg dynamics in situ.

Methane-cycling communities show coupled dynamics and shared responses to Hg contamination

Recently, the role of methane-cycling microbes in mediating MeHg dynamics has become increasingly evident through controlled laboratory studies (14, 22, 81–83). However, the environmental relevance and field-scale implications of these laboratory-derived mechanisms remain largely unexplored in complex agricultural systems.

We show that despite contrasting Hg contamination sources (GX: atmospheric deposition; SK: historical smelting residues), methane-cycling communities at GX and SK exhibited similar beta-diversity (PERMANOVA, P < 0.001; Fig. 3C and D). Constrained ordinations showed that geochemical variables explained 27.67% of the beta-diversity in methanogens and 22.14% in methanotrophs (P = 0.001 and 0.004, respectively; Fig. 3C and D). THg and MeHg were significant predictors for both groups, and DOM concentration showed an additional significant association with methanogen beta-diversity (permutation tests, P < 0.05; Result S3).

*Community structure and environmental drivers of methanogens and aerobic methanotrophs. (A) Taxonomic distribution and relative abundance of mcrA genes across sampling days. (B) Taxonomic distribution and relative abundance of pmoA/mmoX genes across sampling days. (C and D) db-RDA of beta-diversity in relation to geochemical variables, with panel C showing methanogens and panel D showing methanotrophs. Points represent individual samples colored by site. Arrows indicate environmental variables with significance levels (***P < 0.001, **P < 0.01, and P < 0.05).

Methanogen and methanotroph abundances were significantly and positively correlated (β = 0.51, P = 0.003, n = 50; Fig. S6), indicating that their populations tended to co-vary across samples. However, methanotrophs were consistently present at roughly an order of magnitude lower abundance than methanogens (Fig. 3A and B). This observation aligns with studies demonstrating positive co-occurrence networks between these groups (84) and synchronized responses to environmental factors in rice paddies (85), likely reflecting their metabolic interdependence. Methanotrophs rely on methane produced by methanogens and may, in turn, help regulate local oxygen concentrations (84, 85). This interdependence is further supported by the notably lower abundances of both groups at HX, the site with the highest average ORP, where conditions appeared too oxidizing for robust methanogenesis and suboptimal for aerobic methanotrophy, resulting in reduced populations of both metabolically linked functional groups.

The similarity of methane-cycling community structures at the Hg-impacted sites, in addition to their significant association with THg and MeHg in db-RDA analysis, suggests two possibilities. The presence of a core methane-cycling community that has either adapted similarly to the Hg stress or is functionally robust under the high Hg impacts, regardless of contamination source. However, the selective pressures exerted by Hg exposure on methane-cycling guild assembly remain poorly characterized, as do potential adaptations within different methanogenic and methanotrophic lineages to Hg toxicity. These mechanisms warrant further investigation to determine whether community similarity reflects true adaptive convergence responding to Hg stress.

Context-dependent methanogen-MeHg relationships link to DOM quality and compartment effects

We subsequently examined how methanogen abundance (via mcrA) links to MeHg concentrations across sites and soil compartments. Overall methanogen abundance showed marginal association with MeHg (β = −1.30 × 10^7^, P = 0.064), but this relationship varied with THg concentration (β = 4.57 × 10^2^, P = 0.049), indicating context-dependent relationships where higher methanogen abundance amplifies MeHg production under elevated THg conditions.

Our field data reveal a significant correlation between methanogen community composition and DOM characteristics (Fig. 3C), suggesting a potential linkage between DOM, methanogens, and MeHg production. Specifically, methanogen abundance was negatively correlated with both DOM quantity (β = −5.94 × 10^-8^, P = 0.007) and SUVA_254_ (β = −2.88 × 10^-8^, P = 0.003), indicating that methanogen populations may thrive in conditions with lower DOM concentrations and less aromatic (more labile) DOM (10). Interestingly, while aromatic DOM (indicated by higher SUVA254 values) is generally understood to increase Hg(II) bioavailability and methylation by stabilizing HgS complexes (66, 72), we observed a negative correlation between SUVA_254_ and MeHg (β = −0.48, P = 0.005). This contradictory relationship suggests a potential role for methanogens in mediating MeHg production, whereby DOM with lower SUVA254 promotes methanogen activity, leading to enhanced MeHg production. However, this hypothesis requires further investigation.

Specific methanogen taxa showed compartment-dependent patterns. Methanoregulaceae, the dominant family across sites with diverse Hg methylators (7) (Fig. 3A; Result S4), mirrored overall methanogen trends with positive MeHg correlations in SK surface soil (ρ = 0.745, P = 0.011; n = 11) but negative correlations in SK rhizosphere (ρ = −0.943, P = 0.017; n = 6). These family-specific mcrA patterns differed from the methanogen-associated hgcA relationships reported in Fig. 2. This discrepancy likely reflects that methanogens can contribute to both methylation and demethylation processes through direct enzymatic activity or indirect metabolic effects. Consequently, mcrA abundance, which reflects overall methanogen abundance, may not necessarily correlate with hgcA abundance or MeHg production in the same way.

Minor taxa showed different compartment preferences. Methanospirillaceae (ρ = 0.974, P = 0.005; n = 6) and Methanosarcinaceae (ρ = 0.941, P = 0.005; n = 6) correlated positively with MeHg in SK rhizosphere, with weaker patterns at HX and GX. Metatranscriptomic analysis (n = 6) detected active mcrA expression across sites, with higher levels at Hg-contaminated sites, particularly by Methanoregulaceae, providing functional support for gene abundance data (Fig. S4C).

Interestingly, we recovered sequences from Methanoperedenaceae, an archaeal group (ANME), which performs anaerobic oxidation of methane via reverse methanogenesis (86). This finding corroborates the detection of Methanoperedenaceae-associated hgcA from our rice paddy metagenomes (Fig. 2A; Result S7). Although mcrA and hgcA transcripts associated with Methanoperedenaceae were detected at the Hg-impacted sites (Fig. S4C), we exclude this group from subsequent analyses due to the lack of empirical evidence for their methylation capacity. Their potential role in Hg cycling is discussed in Discussion S1.

Flooding links higher methanotroph abundance to lower MeHg

Methanotroph abundance relationships with MeHg followed similar site-specific patterns to those observed for methanogens, consistent with the two guilds co-varying across samples (Fig. 3A and B). However, when focusing specifically on the flooded period, multiple regression revealed a significant negative association between methanotroph abundance and MeHg concentrations (β = –3.21 × 10⁷, P = 0.009), after adjusting for site and soil compartment, an effect not observed in the methanogen data.

Metagenomic analyses revealed that Methylomonadaceae, Beijerinckiaceae, and Methylococcaceae were the most abundant methanotroph taxa (Fig. 3B; Result S5). Metatranscriptomic data further confirmed active methanotrophy at all sites, with particularly high pmoA expression at SK, where Methylomonadaceae was the most highly expressed methanotroph family (Fig. S4D). Methylomonadaceae showed the strongest negative correlations with MeHg at GX (ρ = −0.8, P = 0.133; n = 5) and SK (ρ = −0.9, P = 0.083; n = 5) rhizosphere, while other methanotroph taxa showed weaker correlations (Result S6).

It is likely that the negative methanotroph–MeHg relationship reflects two non-mutually exclusive mechanisms. First, flooding-induced oxygen depletion may suppress aerobic methanotrophs while simultaneously stimulating anaerobic Hg methylation, creating conditions favorable for MeHg accumulation. Alternatively, taxa such as Methylomonadaceae, which are adapted to low-oxygen conditions through their capacity for denitrification (87–89), may persist in flooded, microaerophilic rhizospheres and actively demethylate MeHg via a methanobactin-related pathway, which we explore in the next section.

Exploratory analysis of mbnT as a potential demethylation biomarker

Recent studies indicate a methanotroph-specific MeHg demethylation pathway facilitated by methanobactins (24, 25). Methanobactins are copper-sequestering chalkophores essential for particulate methane monooxygenase activity in some methanotrophs (90). While primarily binding copper with high affinity, methanobactins also bind other metals, including Hg(II), albeit with lower affinity (90).

The cellular uptake of methanobactins is mediated by MbnT, a TonB-dependent transporter (TBDT), which is an outer membrane protein (91). This uptake is crucial as methanobactin appears to deliver MeHg inside the cell, and a study shows that deleting mbnT impairs methanobactin uptake and, consequently, MeHg degradation (24, 25). The degradation of MeHg is thought to be carried out by the periplasmic methanol dehydrogenase, targeting the carbon-Hg bond, and is distinct from the merB-dependent organomercurial lyase pathway, as these methanotrophs generally lack the merB gene (24, 25).

Given the prevalence of methanotrophs in rice paddies, we investigated whether their observed association with MeHg dynamics is linked to MbnT-mediated demethylation. We examined the genomic context of previously studied demethylating methanotrophs and methanotroph-MAGs recovered from rice paddies to determine if they encode genes related to methanobactin biosynthesis (i.e., mbnA, mbnB, and mbnC as core components, and additional genes like mbnD, mbnF, mbnN, mbnS, and mbnX) and transport (mbnE, mbnM, and mbnT) (47).

We confirmed that the demethylating Methylosinus trichosporium OB3b encodes the full suite of genes for methanobactin biosynthesis and uptake/transport in line with previous research (Fig. 4) (24). Methylomicrobium album BG8 and Methylocystis sp. strain Rockwell, which are also able to demethylate, encode mbnT but lack genes necessary for methanobactin biosynthesis. Conversely, Methylococcus capsulatus Bath lacks genes for both methanobactin biosynthesis and uptake (mbnT) and was shown to not demethylate (Fig. 4) (24). Subsequently, we show that none of the methanotroph MAGs possessed mbnT genes sufficiently similar to those found in experimentally confirmed demethylating methanotrophs, and these MAGs also lacked all essential genes for methanobactin biosynthesis. Therefore, the rice paddy MAGs are unlikely to exhibit methanobactin-dependent MeHg demethylation activity, although some of these MAGs showed significant correlations with MeHg via network analyses (Result S8 and Discussion S2).

Presence/absence of methanobactin-related genes in methanotroph genomes. This heatmap illustrates the presence and absence of methanobactin biosynthesis and transport genes across experimentally studied methanotrophs (24) and MAGs recovered in this study, along with their respective quality metrics. The x-axis denotes the methanobactin genes, while the left y-axis represents the taxonomic classification of each genome, and the right y-axis lists the RefSeq or MAG IDs. All RefSeq genomes lack methanobactin biosynthesis, except OB3B (Note mbnB is present in OB3B as a pseudogene, which is not recovered through the HMM-based analyses, but confirmed with RefSeq annotation).

Based on experimental evidence, we explored whether mbnT could serve as a potential indicator for MeHg demethylation capacity, as the demethylating methanotrophs are shown to possess the gene. We recovered mbnT that resembles those found in known demethylating methanotrophs using an HMM (TIGR01783, sequence cutoff 223.5). However, recognizing that TBDTs, including mbnT, have homologs in diverse organisms beyond methanotrophs, we employed a conservative approach, retrieving sequences classified to known demethylating methanotroph genera: Methylomicrobium, Methylocystis, and Methylosinus (Fig. S7). Transcriptomic analyses show that mbnT transcripts from these genera were present in Hg-contaminated rice paddies, though to a limited extent (n = 6, Fig. S4E). However, multiple regression revealed that total mbnT abundance from these genera was a significant positive predictor of MeHg concentration (β = 4.579 × 10⁷, P = 0.0136) after accounting for site and compartment specific effects (Result S9), contrary to the expected negative correlation we would expect from an effective demethylation marker.

Our findings suggest that while methanotroph demethylation capacity appears to be present in rice paddies, utilizing mbnT as a marker to assess net MeHg accumulation may not be reliable without consideration of several limitations. First, mbnT is not a single-copy gene (92), and it belongs to the family of TBDTs, which are responsible for the uptake of diverse substrates, including metal chelates, and have numerous homologs in organisms other than methanotrophs (47). Despite using a stringent cutoff, the presence of mbnT homologs with high similarity could lead to biased abundance estimations. Furthermore, gene abundance does not directly equate to activity. Second, methanotrophs exhibit varied demethylation rates (24). These rates depend on various biological factors, including the structure of methanobactin (24), the efficiency of the MbnT transport system (90), and the presence of additional copper-handling proteins, like MbnH and MbnP (91), as well as environmental factors such as redox conditions, pH, and the presence of competing metals (93–95). Methanobactins can influence Hg cycling through secondary mechanisms. This includes altering microbial community composition (e.g., through copper competition impacting denitrifying bacteria [90]) or affecting the bioavailability of Hg species via competitive binding (95, 96). These broader effects complicate the interpretation of *mbnT’*s role in MeHg demethylation.

Future studies should prioritize establishing direct relationships between mbnT expression and demethylation rates in situ through quantitative approaches with adequate sample sizes, before linking mbnT abundance to environmental MeHg dynamics. Additionally, broader experimental validation across diverse methanotrophic taxa is needed to confirm the universal applicability of this biomarker. Given the complexity revealed by our environmental data, controlled laboratory studies examining mbnT expression under varying Hg, redox, and competing metal conditions would help define when mbnT serves as a reliable indicator of methanotroph-mediated MeHg demethylation.

Methane-cycling microbes as key biological predictors of MeHg variability

Our analyses demonstrate that relying on single factors has limited explanatory power for net MeHg accumulation, despite significant associations found for individual variables. The complexity expected from concurrent methylation and demethylation processes is reflected in the increase in model R² when accounting for site and compartment effects (i.e., from 60% to 84% in hgcA models). To address the limitations of single-factor analyses and identify the key, potentially non-linear, variables driving MeHg patterns across diverse functional genes and varying environmental conditions, we employed RF analysis. RF allowed us to identify non-linear relationships and interactions among variables by building multiple decision trees on training data and assessing their predictive power on held-out test samples (97).

The RF models built using all variables achieved moderate predictive performance (R²=0.70 ± 0.13, range: 0.40–0.95, Fig. 5). Bioavailable Hg(II) (F1-Hg) consistently emerged as the dominant predictor across all data splits (mean IncMSE = 30.9, 100% significance rate, P = 0.010). Methanogen abundance ranked as the second most important variable (mean IncMSE = 20.8, 100% significance rate, P = 0.013), followed by methanogen-hgcA abundance (mean IncMSE = 9.3, 49% significance rate), substantially outperforming other microbial indicators like total Hg methylators (hgcA), SRB, and IRB abundances.

Variable importance and model performance for predicting MeHg concentrations using RF analysis (full model). (A) Variable importance plot showing the mean percent increase in mean squared error (%IncMSE) for each predictor variable when permuted, with 95% confidence intervals across 100 different data splits. Larger point size and green coloration indicate higher significance rates. Values show the percentage of data splits where the variable was significant (P < 0.05), followed by the mean P-value. (B) Model performance metrics displaying mean values with 95% confidence intervals across all data splits. The model was built using 2,001 trees with mtry = 4 and validated through 100 permutation tests with random 70/30 train/test splits.

Methanotroph abundance was identified as the third most important biological variable overall (mean IncMSE = 8.8; 31% significance rate). An explicit interaction term between methanogen and methanotroph abundances (i.e., their product) was a critical predictor of MeHg (mean IncMSE = 18.3; 80% significance rate; P = 0.038). This result suggests that the methanotrophic influence is primarily expressed through community-level coupling with methane-producing populations. Given that methanogens are key methylators that amplify MeHg under elevated THg, and methanotrophs showed a significant negative association with MeHg during the flooded period (suggesting a role in MeHg consumption/demethylation), the interaction term highlights that the net MeHg outcome is highly sensitive to the synchronous activity of these two interdependent microbial guilds.

Constraining the RF model solely to the flooded period substantially improved predictive accuracy (R^2^ = 0.86 ± 0.08; Fig. S8B). This improved performance was accompanied by a modest increase in the relative importance of methanogens (mean %IncMSE rose from approximately 21 to 23) and a decrease for F1-Hg (mean %IncMSE fell from approximately 32 to 26), suggesting that the predictive power of methanogens is comparatively greater under flooded, reduced conditions than that of Hg bioavailability alone. Integrating genomic and geochemical data through RF modeling revealed that methanogen, methanogen methylator, and methanotroph abundance were the top biological factors contributing to net MeHg accumulation, while Hg bioavailability and DOM aromaticity were the dominant geochemical factors. Prediction accuracy further improved when removing less influential variables, reinforcing the importance of these key factors.

While previous studies have associated methanogenesis with MeHg content in rice paddies (14, 21), our work advances this understanding by identifying methanotrophy as an important biological correlate of net MeHg accumulation and characterizing specific methane-cycling populations associated with MeHg variability. Methanotrophic involvement is particularly noteworthy, as it suggests a previously underappreciated potential demethylation pathway in rice paddies. Our findings indicate that studying the relationship between microbial methane metabolism and Hg cycling provides critical insights into the key biogeochemical factors associated with net MeHg accumulation in Hg-contaminated rice paddies.

While these findings advance our understanding of microbial associations with MeHg dynamics, several considerations should guide their interpretation and application. The observed coupling between methane and Hg cycling, though mechanistically plausible, requires validation across diverse environmental contexts before broader generalization to other rice paddy systems. Our RF approach, while effective at identifying statistical associations, cannot establish causal links between predictors and MeHg accumulation, highlighting the need for controlled experimental validation of these relationships. We recommend validation across diverse rice paddy systems with (i) varying flooding patterns and water management practices (12, 98), (ii) different rice cultivars (99), and (iii) diverse geographical span (3) and contamination sources (100–102) to benchmark model transferability.

Our results underscore that effective assessment of MeHg risk requires integrating microbial community data with physicochemical gradients that govern both Hg bioavailability and microbial niche occupation. The importance of methane-cycling communities in our predictive models highlights promising avenues for management strategies that could simultaneously address climate change mitigation and Hg pollution through targeted manipulation of shared microbial pathways (103).