Phylogenetic and ecophysiological novelty of subsurface mercury methylators in mangrove sediments
Introduction
The production of neurotoxic methylmercury (MeHg) by microorganisms is mediated by the functional proteins encoded by the two-gene cluster hgcA and hgcB [1]. The presence of the hgcAB gene pair is a reliable indicator of a microbe’s ability to methylate mercury (Hg) [2], and has previously been linked to sulfate-reducing bacteria (SRB) [1], iron-reducing bacteria (IRB) [3], methanogens [4], and other anaerobes [2]. The discovery of the hgcAB gene pair has significantly expanded the known diversity of microorganisms capable of producing MeHg [5]. For instance, recent studies investigating the subsurface biosphere [6,7,8] have identified Hg-methylating microbes from several phyla, including Lokiarchaeota, Aminicenantes, Kiritimatiellaeota, Spirochaetes, Marinimicrobia, and Actinobacteria, that were previously not thought to be associated with MeHg production. The subsurface biosphere is a major reservoir of phylogenetically and metabolically diverse microorganisms that have largely yet to be characterized [9]. Despite the detection of bioaccumulated MeHg in this environment [2], the microbiota responsible for mediating MeHg cycling in the subsurface biosphere remains an elusive “dark matter”.
Mangrove sediment is a crucial component in the global Hg cycling, serving as ideal niches of microbial Hg methylation [10]. Previous research has primarily focused on MeHg biogeochemistry in the superficial sediments of mangroves, revealing the presence of well-known Hg methylators [11]. Recent findings have demonstrated that the subsurface microbial communities showed apparent acclimation to Hg, and Hg methylation may be favored by more reducing conditions (i.e., in the deeper layers) [10, 12]. However, whether the more anoxic and highly reducing condition in mangrove subsurface sediments can lead to intensive MeHg production remains unknown, as does the identity of the subsurface microorganisms responsible for Hg methylation. Our prior research, indicating that the more reducing conditions in the subsurface layer favor anaerobes (i.e., SRB and methanogens) [13], raises the possibility of a more diverse Hg-methylating microbial community in the mangrove subsurface sediments.
Mangrove subsurface sediments have long been recognized as energy-limited [14,15,16], since more recalcitrant organic matters are trapped in deeper sediments with physical forces (i.e., hydrodynamics and gravity) [17]. However, owing to the lack of cultured relatives, much remains unknown about the metabolic feature and adaptive strategy of subsurface Hg methylators under such energy-limited conditions. In recent years, high-throughput culture-independent sequencing has provided a window into the metabolic novelty of subsurface-dwelling microorganisms [18,19,20]. The discovery of novel methanogen genomes from the terrestrial subsurface has uncovered their diverse ecophysiological adaptations with unexpected metabolic features (i.e., the potential for extracellular electron transfer) [21]. In energy-limited conditions, syntrophic relationships have been proposed as a common strategy for mediating biogeochemical cycles [22, 23]. A recent study has provided evidence that the SRB-methanogen syntrophy via interspecies electron transfer dramatically stimulated Hg methylation in co-cultures [24]. However, the potential for syntrophic relationships in energy-limited mangrove subsurface sediments for Hg methylation remains uncertain. With the continued advancement of high-throughput sequencing, the integration of microbial genomic information into syntrophic models will enhance our understanding of the adaptive strategies of Hg methylators in the subsurface biosphere.
In this study, we used MeHg measurements, Hg isotope assays, and metagenomic sequencing to examine the levels of MeHg, phylogenetic diversity, and metabolic adaptations of Hg methylators in 100-cm mangrove sediment columns. Our findings indicate elevated diversity of Hg methylators in subsurface layers. This prompted us to investigate the regional-scale distribution of these subsurface-dwelling Hg methylators across southeastern mangrove regions of China, spanning over 1500 km. Our most significant discovery is the identification of the candidate phylum Zixibacteria as a currently understudied Hg methylator in anoxic subsurface environments, with unexpected metabolic and lifestyle features, such as Hg methylation independent of the Wood–Ljungdahl pathway but involved in methionine synthesis, benzoyl-CoA, and syntrophic pathways. This study significantly expands the diversity of subsurface-dwelling Hg methylators and sheds light on their unique metabolic adaptations in mangrove sediments.
Materials and methods
Mangrove sediment sampling
The sampling site is located at Qi’ao Mangrove Wetland Park (22° 26′ 12.28′′ N, 113° 38′ 26.12′′ E), Guangdong Province, China. The park was characterized by a mean annual temperature of 22.4 °C and annual precipitation range of 1700–2200 mm. The tidal regime was irregular semidiurnal, with average high and low tide levels of 0.17 m and −0.14 m, respectively [25]. We collected partially air-exposed sediments from the native mangrove species (Kandelia obovata) in the park, in August 2019. Five replicate sediment cores were obtained using a 1-m long PVC sampling column after ebb. The cores were sliced into 10 depth intervals (0–5, 5–10, 10–15, 15–20, 20–30, 30–40, 40–50, 50–60, 60–80, and 80–100 cm), yielding 50 samples in total. The sliced sediments were stored in a portable cooler at 4 °C and transported back to the laboratory within 24 h for further analysis. Each sample was divided into two sub-samples for different purposes: one for biogeochemical properties analysis, stored at 4 °C, and the other for DNA extraction, stored at −80 °C.
Biogeochemical properties analyses
Sulfate anion in the sediments was analyzed using a Thermo Scientific Dionex ICS-1100 Ion Chromatography System, which consisted of a guard column (IonPac AG14A), an analytical column (IonPac AS14A), a suppressor (Dionex ERS 500), and a column heater (Dionex ICS-1100 Column Heater). CH4 measurement was carried out by placing the sediment samples in 20 mL serum bottles containing 5 mL of 1 M NaOH-buffered deionized water [26]. The porewater CH4 concentrations were then measured using an online GC–Mass (GCMS-QP2010 SE). Fe (II) concentrations were measured directly, whereas Fe (III) concentrations were calculated as the difference between total Fe and Fe (II) concentrations. Total Fe and Fe (II) were determined using iron-phenanthroline spectrophotometry, which relies on the Fe2+-dependent optical absorbance of bathophenanthrolinedisulfonic acid (BPS) at 535 nm at pH 5.4 [27].
Mercury analysis
All analyses for Hg and MeHg concentration as well as total Hg (HgT) stable isotope were performed as described in previous studies [28, 29]. The homogenized and freeze-dried sediment samples were prepared for Hg concentration and isotope analysis using a double-stage tube furnace in conjunction with 40% anti-aqua regia trapping solution (HNO3/HCl = 2/1, v/v) [30]. All digests were heated to 95 °C for 8 h. Standard reference materials (GSS-5, Hg concentration: 290 ± 30 ng/g) and method blanks were processed similarly to the samples. The former yielded Hg recoveries of 109 ± 4% and the latter showed Hg concentrations lower than the detection limit, precluding lab contamination. The HgT and MeHg contents in the digested sample were measured by a MERX automated HgT analytical system (Brooks Rand Laboratories, USA) and inductively coupled plasma mass spectrometry (ICP-MS) according to USEPA standard methods 1631E and 1630, respectively [31, 32].
For HgT stable isotope analyses, the preconcentrated solutions were diluted to 0.5 ng/mL with an acid concentration of 10–20% prior to Hg isotope analysis using a Neptune Plus MC-ICP-MS. Hg isotope ratios were reported following the convention proposed by Blum and Bergquist [33]. To ensure that pre-concentration processes did not induce fractionation [34], the National Institute of Standards and Technology (NIST) 3133 was also run before and after each sample. The mass-dependent fractionation (MDF) is expressed in δ202Hg notation in units of ‰ referenced to the NIST-3133:
The mass-independent fractionation (MIF) is reported in Δ notation, which describes the difference between the measured δxxxHg and the theoretically predicted δxxxHg value in units of ‰:
β is 0.252 for 199Hg, 0.5024 for 200Hg, and 0.752 for 201Hg. NIST-3177 secondary standard solutions, diluted to 0.5 ng/mL Hg with 10% HCl, were measured before and after each sample. The overall average and uncertainty values of NIST-3177 (δ202Hg: −0.53 ± 0.06‰; Δ199Hg: −0.15 ± 0.01‰; Δ200Hg: 0.02 ± 0.03‰; Δ201Hg: 0.00 ± 0.03‰; 2 SD, n = 8) and GSS-5 (δ202Hg: −1.76 ± 0.03‰; Δ199Hg: −1.60 ± 0.02‰; Δ200Hg: 0.03 ± 0.03‰; Δ201Hg: −0.29 ± 0.03‰; 2 SD, n = 8) agree well with previous results [28, 29]. The larger values of 2-SD between NIST-3177 and GSS-5 are used to reflect maximum analytical uncertainties.
For Hg isotope fractionation, the previous study summarized the research published up until the year 2011, with the aim of providing a basic framework for understanding MDF and MIF in the environment [35]. Nearly all kinetic reactions involving Hg produce products with lower δ202Hg (i.e., isotopically lighter) and leave a residual pool of reactants with higher δ202Hg (i.e., isotopically heavier). In contrast, microbial reactions do not produce significant MIF (i.e., Δ199Hg does not differ between products and reactants). As these reactions proceed, the ratio of Δ199Hg to δ202Hg can be diagnostic of the type of reaction. The overview of the general patterns in isotope fractionation is detailed in Supplementary Table S1.
DNA extraction
The whole community DNA was extracted from 5.0 g mangrove sediment and purified using a combined protocol of sodium dodecyl sulfate extraction method [36] and Power Soil DNA Isolation Kit (Mo Bio Laboratories, Carlsbad, California, USA) according to the manual. The DNA purity was checked by NanoDrop ND-2000 Spectrophotometer (Thermo Fisher Scientific, MA, USA), and ratios of 260/280 and 260/230 were about 1.80 and 1.70, respectively. The DNA concentrations were quantified using a fluorescent method (Qubit 4 Fluorometer, Thermo Scientific, USA).
Metagenomic sequencing and processing
For metagenomic sequencing, 1.0 μg of DNA from each sample was prepared using the NEBNext DNA Library Prep Kit. Samples were individually sequenced with the NextSeq550 platform (Illumina, San Diego, CA, USA), resulting in 2 × 150 bp paired-end reads with an average of 12.47 Gbp per sample. Low-quality (quality score ≤ 38, base N > 10 bp, the overlap length between adapter and reads >15 bp) paired-end reads were filtered. After trimming the raw sequencing reads using BBDuk (https://jgi.doe.gov/data-and-tools/bbtools/, parameters: ktrim = r k = 28 mink = 12 hdist = 1 tbo = t tpe = t qtrim = rl trimq = 20 minlength = 70), genome assembly and binning were performed according to the MetaWRAP pipeline v1.3.0 [37]. The sequences were assembled with MEGAHIT (v1.1.3; parameters: -m 500) to generate contigs [38]. Genome binning of assembled contigs was done using MetaBAT2 (v2.12.1) [39] and MaxBin2 (v2.2.6) [40], and the resulting metagenome-assembled genomes (MAGs) were consolidated with the Bin_refinement module. The Reassemble_bins module further improved the consolidated bin sets to generate MAGs (parameters: -m 800). The quality of MAGs was evaluated with CheckM v1.0.12 [41]. Only MAGs with a completeness of over 50% and a contamination level below 10% were analyzed further. Representative MAGs were chosen using dRep (v3.4.0, parameters: dRep dereplicate -sa 0.99 -nc 0.1) [42]. The abundance of each MAG was expressed as genome copies per million reads and calculated using CoverM v0.6.1 genome mode (https://github.com/wwood/CoverM). Taxonomic assignments of the MAGs were performed using the GTDB-Tk v1.3.0 [43].
Search for hgcAB genes
In this study, we utilized the Hg-MATE-Db [44], an up-to-date database of Hg methylators from pure culture and environmental datasets, to search for hgcAB genes. The search was conducted using HMMER v3.3.2 [45], and the results were further confirmed by examining the presence of conserved motifs (N(V/I)WCA(A/G) in hgcA and (CX2CX2CX3C) in hgcB, respectively. The hgcAB sequences were clustered at 100% identity using CD-HIT to remove redundancy [46].
Relative and normalized hgcA abundance calculation
To quantify the relative abundance of hgcA gene, we calculated the average number of aligned reads that overlapped each position on the contig. The hgcA index was built using the Salmon v0.8.1 index mode [47]. Reads were mapped to the hgcA using Salmon in the quant mode. To analyze the normalized hgcA abundance, we needed to estimate total cell counts, whereas sequencing data can generate only relative abundance values (in the absence of an internal standard). We began with the general equation:
This approach has been discussed and tested in-depth previously [48], where DNA yield was utilized as an estimate of microbial load in a sample. DNA yield is not a perfect estimate of microbial load due to potential ploidy [49] and bias in sequencing representation depending on the DNA extraction method [50]. But in this study, all samples were extracted using the same DNA extraction kit, and DNA yield was measured using a fluorometry-based method (Qubit 4 fluorometer dsDNA HS Assay), making it the best available measurement of microbial load (Supplementary Fig. S1).
Phylogenetic analyses
Hg methylator genomes and hgcA genes recovered in this research were analyzed and compared with reference sequences from Hg-MATE-Db. MAFFT (–auto) [51], trimAl (v1.4.rev15) [52], and IQ-TREE 2 (parameters: -st AA -alrt 1000 -bb 1000) [53] were used for sequence alignment, trimming, and construction of phylogenetic trees, respectively. The phylogenetic trees for Hg methylator genomes were inferred by a concatenated set of specific single-copy marker genes in the GTDB (https://gtdb.ecogenomic.org/). The orthologs of these marker genes in the Hg methylator and the reference genomes were identified using GTDB-Tk based on HMMER. The results were edited and modified using iTOL [54] and Adobe Illustrator, respectively. The pairwise genomic ANI values were calculated by FastANI (v1.33, https://github.com/ParBLiSS/FastANI) with default parameters based on the assemblies [55]. FastANI was targeted to estimate ANI in the 70–100% identity range.
Functional characterization of genomes
Gene prediction for all MAGs was performed using Prodigal v2.6.3 (default settings) [56]. The predicted genes were characterized using McycDB [57], SCycDB [58], and FeGenie [59]. MAGs and reference genomes were also annotated using the KEGG-based annotation program METABOLIC v4.0 [60], and custom databases were searched using DIAMOND v2.0.15 [61] and HMMER. Hydrogenases were identified using METABOLIC, consisting of NiFe-, FeFe-, and Fe-hydrogenase catalytic subunits [62]. Putative hydrogenase sequences identified in the search were then uploaded to the HydDB webserver [62] to identify and remove non-hydrogenases. Anaerobic hydrocarbon degradation genes were identified using a HMMER database (hmmbuild, hmmpress, and hmmsearch modes) of reference sequences from the anaerobic hydrocarbon degradation database [63]. Prokaryotic Hg resistance genes (merAB) were identified using publicly available homologs of merAB [64] and the search custom-markers function using HMMER v3.3.2. The distribution of biogeochemically important genes identified using these methods is shown in Supplementary Data 1. To detect multi-heme cytochromes, we utilized the Python script “cytochrome_stats.py” available for download at https://github.com/bondlab/scripts [65]. This script counts the occurrence of the heme-binding motif, CXXCH, in protein sequences. Proteins with more than three CXXCH motifs were considered as multi-heme cytochromes, following the established criteria [65, 66]. The detailed information of multi-heme cytochromes is shown in Supplementary Data 2. The annotation of Wood-Ljungdahl pathway, methionine synthesis, and benzoyl-CoA reduction was based on previous studies [67,68,69,70]. The typical cobalamin-dependent methyltransferases were predicted using KofamKOALA (release 106.0, https://www.genome.jp/tools/kofamkoala/, parameters: E‐value ≤ 1e−5) [71, 72].
Protein AlphaFold2 modeling
The three-dimensional structures of the putative HgcAB and MetH protein sequences found in this study were built and refined upon the official release of the source code and monomer neural network models of AlphaFold2 (version 2.0.1) [73]. The transmembrane regions in HgcA were predicted by CCTOP [74]. Multiple structure alignment for different taxa was performed using mTM-Align [75]. The cofactor cobalamin was placed in the binding pocket by superposing the HgcA and CFeSP (PDB entry 4DJF) [76]. Consistent with the expected Cys coordination patterns from other dicluster ferredoxins, such as that from Clostridium acidurici (PDB entry 2FDN) [77], preliminary de novo models of HgcB with coevolution restraints suggest that one [4Fe-4S] cluster was bound to Cys20, Cys23, Cys26, and Cys60 and another was bound to Cys50, Cys53, Cys56, and Cys30. The preliminary models of the HgcAB complexes were generated by PyMOL 2.5.2 (www.pymol.org), and complexes were energy minimized using the Schrodinger Suite 2021-4 (Schrödinger, LLC, Portland, OR) energy minimization server. The C-terminal tail of HgcB was also introduced at this step. The alignment of HgcA AlphaFold2 models between Zixibacteria and six confirmed Hg-methylating bacterial strains was performed by PyMOL 2.5.2.
Regional-scale distribution of Hg methylators in mangrove subsurface sediments along southeastern coastal China
A total of 27 samples were collected in mangrove sediments along southeastern coastal China, with 7 samples from Hainan (5 replicates for 30–40 cm, 2 replicates for 80–100 cm), 10 samples from Fujian (5 replicates for both 30–40 cm and 80–100 cm), 10 samples from Zhejiang (5 replicates for both 30–40 cm and 80–100 cm). The sampling sites were located at Bamen Bay Mangrove National Wetland Park (19°33′25.2“N, 110°50′20.40″E) with a mean temperature of 26.5 °C in Hainan Province, Zhangjiang Estuary Mangrove National Nature Reserve (23°55′17.18″N, 117°25′11.04″E) with a mean temperature of 28.5 °C in Fujian Province, and Ximen Island (28°20′45.15′′N, 121°10′58.27′′E) with a mean temperature of 27.0 °C in Zhejiang Province, China, respectively. The methods for DNA extraction and initial metagenomic analysis to obtain MAGs were conducted as described above. Also, we performed a series of analyses, including prodigal, hgcAB genes annotation, and conserved motifs examination, on the MAGs obtained from these three sampling points to examine the presence of Hg methylators, following the methods described above.
Global distribution
The distribution of representatives of the Zixibacteria lineage across various environments was estimated following the methodology outlined in [78]. In brief, the longest 16 S rRNA gene sequence in Zixibacteria MAGs obtained from this study was utilized for a homology-based search against the 16 S rRNA Public Assembled Metagenomes database of the Integrated Microbial Genomes and Microbiomes (IMG/M v7) system [79], employing the IMG BLAST Tool. BLAST hits with sequence identity greater than 75% (phylum threshold) [80], and the associated metadata (longitude, latitude, and habitat type) were subsequently retrieved from IMG for further analysis.
HgcA-MetH protein interaction simulation and analysis
The HDOCK server (https://hdock.phys.hust.edu.cn/) was utilized to predict a binding conformation between HgcA and MetH [81]. The HDOCK algorithm incorporates both template-based and free approaches, providing cross-validation [81]. This method involves an attempt to sample all conceivable binding modes between two separate protein structures. Given the typical absence of binding site information, a global docking approach is employed to sample assumed binding modes, involving six degrees of freedom (three rotational and three translational).
Subsequently, a scoring function, comprising docking and confidence scores, was applied to rank the sampled binding modes based on the sampling process and the resulting binding conformations [82, 83]. The docking score signifies the likelihood of a binding model, with more negative scores indicating a higher probability of binding models. The confidence score was computed using the formula:
A confidence score exceeding 0.70 indicates a strong likelihood of binding between the two molecules. The confidence score between 0.50 and 0.70 suggests a potential binding interaction between the two molecules. A confidence score below 0.50 indicates a lower probability of binding between the two protein molecules. To generate a 2D residue-residue interaction plot depicting the entire interface of a protein-protein complex, we inputted the complex into LigPlot [84].
Co-occurrence analysis
The co-occurrence network was based on the full-length 16 S rRNA amplicon sequencing data. For full-length 16 S rRNA amplicon sequencing, both the forward and reverse primers were tailed with sample-specific PacBio barcode sequences to allow for multiplexed sequencing. The KAPA HiFi Hot Start DNA Polymerase (KAPA Biosystems) was used to perform 20 cycles of PCR amplification, with denaturing at 95 °C for 30 s, annealing at 57 °C for 30 s, and extension at 72 °C for 60 s. Post-amplification quality control was performed on a Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Amplified DNA from the sediment samples was then pooled in equimolar concentration. SMRTbell libraries were prepared from the amplified DNA by blunt ligation according to the manufacturer’s instructions (Pacific Biosciences). Purified SMRTbell libraries from the mangrove sediment bacterial communities were sequenced on PacBio Sequel cells. The DADA2 method is an algorithm for the inference of the amplicon sequence variants (ASVs) present in a sample from the library of noisy reads generated by amplicon sequencing [85]. The standard processing steps in the DADA2 workflow include quality filtering, dereplication, learning the dataset-specific error model, ASV inference, chimera removal, and taxonomic assignment against a SILVA138 database [86]. After processing the full-length 16 S rRNA amplicon sequencing data, we investigated the co-occurrence network using FastSpar (v1.0) algorithm [87]. All p values were corrected using the Benjamini–Hochberg multiple correction method. Afterward, all the statistically significant correlated taxa (p < 0.001) were used to infer the network and visualized using Gephi (version 0.10.1, https://gephi.org/) [88].
Statistical analysis
Statistical analysis and plotting were performed using R software version 4.2.1 (www.r-project.org) and associated packages, including vegan, spaa, ggtrendline, gggenes, ggmsa, mapchina, ggcor, and ggplot2 [89,90,91,92]. Richness was calculated using the spa package. The analysis of the non-metric multidimensional scaling (NMDS) was performed using the vegan and ggcor packages. Correlations were evaluated by a linear fit and Spearman test using the ggtrendline and vegan packages. Alignments of the hgcAB amino acid sequences were plotted using the ggmsa package, while the genomic context was visualized with the gggenes package. The map of China was produced using the mapchina package.
Results and discussion
Microbial Hg methylation in mangrove sediments
We investigated the variations in HgT and MeHg across five mangrove sediment columns (0–100 cm) collected in Guangdong Province, China. Across the mangrove sediment columns, the average HgT concentration was 420.06 ± 95.49 ng/g, and an increase in the concentration of HgT was observed from the surface to 20–30 cm (Fig. 1A). However, the concentration of MeHg was found to peak at 0–15 cm (3.27 ng/g), followed by a gradual decrease towards the bottom of the sediment columns. The MeHg/HgT ratio, a proxy for the potential HgT available for MeHg production [93], showed a constant decrease with increasing sediment depth, with an average value of 0.67 ± 0.09% (Fig. 1A). This falls within the typical range of 0.11−7.13% observed in natural mangrove ecosystems [10]. The in-depth profiles of sulfate, CH4, and iron (referring to the iron (III)/iron (II) ratio) were found to reflect that of MeHg across the mangrove sediment columns (Supplementary Fig. S2), which are key factors involved in Hg methylation [94].
A The variations of total dissolved Hg and MeHg in the mangrove sediments column. B The relationship between δ202Hg and Δ199Hg in mangrove sediments, with reference data [96,97,98,99,100,101,102] from soil, marine sediments, and seawater ecosystems. The gray shaded area represents the terrestrial reservoir (soil), the light purple shaded area represents the marine reservoir (marine sediments and seawater), and the light orange shaded area represents the mangrove sediments. The trend of δ202Hg versus Δ199Hg with a near-zero ratio (Ratio = −0.16) in mangrove sediments indicates a significant contribution of microbial actions to the Hg methylation. The shadow indicates the 95% confidence interval.
Through Hg isotopic analyses, we investigated Hg sources in mangrove sediments by analyzing surface (0–5 cm), medium (40–50 cm), and deep (60–80 cm) sediment layers (three replicates each). The Δ199Hg, Δ200Hg, Δ201Hg, and δ202Hg values were determined using mass-independent fractionations (MIF) and large mass-dependent fractionations (MDF) [95]. Our data showed that the sediment across different depths exhibited similar δ202Hg ranges (−1.52 to −0.67‰, Fig. 1B and Supplementary Table S2), which were ~14 times larger than the 2 standard deviations (±0.06‰) analytical uncertainty for δ202Hg. This distinctive negative MDF of 202Hg in mangrove sediments was likely due to the mixture of δ202Hg sources in the land-ocean connection [13]. To further support this finding, we compared our results with 247 reported isotope values [96,97,98,99,100,101,102] of Δ199Hg and δ202Hg from terrestrial soil (177), marine sediments, and seawater (70), and found that the isotope values of mangrove sediments were in between the values from terrestrial soil and seawater (Fig. 1B). This result indicates that mangrove sediments received terrestrial- and marine-sourced Hg probably due to complex hydrodynamic (either surface runoff or tides and waves) [103].
Our dataset showed a limited MIF of 200Hg (−0.09 to 0.07‰; Supplementary Table S2), which is close to the analytical uncertainty for Δ200Hg (±0.03‰). As microbial Hg methylation does not result in odd-MIF shifts [104], the near-zero odd-MIF values detected suggest that the Hg methylation in mangrove sediments is predominantly driven by microbial reactions. To verify this hypothesis, we determined the Δ199Hg/δ202Hg ratio, which serves as a diagnostic of Hg methylation reaction types [35], with values of ~1.2–2.4 and ~0 expected for photochemical and microbial reactions, respectively (Supplementary Table S1). Our analyses revealed a strong correlation between Δ199Hg and δ202Hg in the studied samples (R2 = 0.47, p < 0.05), and the sediment Δ199Hg/δ202Hg ratio was −0.16 (Fig. 1B), indicating a substantial contribution of microbial actions to Hg methylation in mangrove sediments.
Composition and distribution of Hg methylators in mangrove subsurface sediments
Through metagenomic analysis, we identified 432 distinct hgcA sequences with conserved motifs (N(V/I)WCA(A/G)) in mangrove sediment columns and quantified their normalized abundance. We observed a notable pattern where the total normalized abundance of hgcA sequences was comparatively elevated above 20 cm and decreased below this depth (as depicted in Fig. 2A). Correlation analyses confirmed a strong association between the total normalized abundance of hgcA sequences and the concentration of MeHg in the sediments (Supplementary Fig. S3; Spearman ρ = 0.79; p < 0.0001). This finding further supports the significant role of microbial Hg methylation in mangrove sediments.
A Normalized abundance and species richness of hgcA gene. B Occurrence frequency of representative hgcA genes from different taxonomic groups. Taxonomic classifications of the Hg methylators are labeled by different colors. C Regional-scale distribution of Hg methylators in mangrove sediments in southeastern coastal China. The samples were collected from Hainan (n = 7), Guangdong (n = 50), Fujian (n = 10), and Zhejiang (n = 10). The number of metagenomic samples is indicated by n. Taxonomic classifications of the Hg methylators are labeled by different colors, which are consistent with Fig. 2B. Among them, the gray color represents the Hg-methylating Dehalococcoidia solely present in mangrove sediment samples from Hainan. White indicates the absence of a taxonomic group at the sampling site. D Energy metabolism redox pairs of each Hg-methylating lineage.
We found that hgcA sequences in the subsurface sediments (below 20 cm) displayed a higher species richness than those in the surface (above 20 cm) (Fig. 2A, B), pointing to an elevated diversity of Hg methylators in the mangrove subsurface sediments. To verify this observation, we conducted a detailed analysis of the composition and variation of Hg methylators across the mangrove sediment columns. Our metagenomic study identified 65 medium- and high-quality MAGs with conserved hgcA (N(V/I)WCA(A/G), Supplementary Fig. S4) and hgcB (CX2CX2CX3C, Supplementary Fig. S5) motifs in the samples collected from all mangrove sediment layers. These MAGs represented a diverse range of taxa (Fig. 2B and Supplementary Data 3), including Acidobacteria (5), Anaerolineae (9), Bacteroidales (2), candidate division Zixibacteria (6; hereafter referred to as Zixibacteria throughout the manuscript), Desulfobacterales (6), GWC2-55-46 (4), Methanomicrobia (6), MBNT15 (7), Nitrospirae (8), Phycisphaerae (6), and Syntrophobacterales (6). We also observed a clear divergence in the distribution pattern of Hg methylators at the partition depth of 20 cm (Adonis test, R2 = 0.19, p < 0.001; Anosim test, R = 0.58, p < 0.001; Supplementary Fig. S6), with most of these MAGs exhibiting elevated occurrence frequencies in subsurface sediments (Fig. 2B).
To assess the regional-scale distribution of these Hg methylators, we collected and analyzed 27 additional mangrove subsurface sediment samples from three other provinces (Zhejiang, Fujian, and Hainan) along the southeastern coast of China, covering a distance of more than 1500 km (Fig. 2C). In addition to 65 MAGs recovered in Guangdong, metagenomic analyses of these 27 samples also generated 30, 15, and 38 hgcA-carrying MAGs in Zhejiang, Fujian, and Hainan, respectively. Three types of Hg methylators (Desulfobacterales, Anaerolineae, and Syntrophobacterales) were detected in all four sampling provinces (Fig. 2C). These members were the dominant Hg methylators in mangrove subsurface sediments, with the cumulative relative abundance of 26.5–91.2% (Supplementary Fig. S7). Owing to their dominance and environmental adaptability [105, 106], these Hg methylators may largely contribute to the ecological interception effect of mangrove sediments on terrestrial- or marine-sourced Hg by microbial methylation [10]. However, hgcA-carrying Zixibacteria, Dehalococcoidia, and GWC2-55-46 were only detected in mangrove subsurface sediments from Guangdong, Hainan, and Guangdong, respectively (Fig. 2C). Together with a recent discovery of Zixibacteria MAGs with hgcA in another mangrove sediment area of Guangdong [8], our findings point toward the geographic selectivity of specific Hg methylators across mangrove sediments.
Metabolic adaptation of Hg methylators in mangrove subsurface sediments
The adaptation of Hg methylators to mangrove subsurface sediments requires tolerance to intermittent low levels of O2 and energy exposure [16]. Our functional annotation revealed that the hgcA-carrying MAGs in mangrove sediments contained multiple terminal oxidases (Fig. 2D; Supplementary Data 1). Nearly all hgcA-carrying MAGs except Methanomicrobia contained the complete operon encoding the cytochrome bd complex, which has a high affinity for O2 and supports adaptation to low-oxygen environments [107]. Furthermore, most members in the GWC2-55-46 and MBNT15 groups contained genes encoding cytochrome c cbb3-type (ccoPQNO) oxidase (Supplementary Data 1), which could exhibit high affinity for O2 and enable bacteria to respire O2 under both oxic and suboxic conditions [108]. The presence of genes involved in anaerobic hydrocarbon degradation, such as those encoding putative anaerobic benzene carboxylase (abcA) and phenylphosphate carboxylase (ppcAB), provides additional evidence for the anaerobic lifestyle of Hg methylators in mangrove subsurface sediments.
Consistent with the previously documented feature of Hg methylators [109], the diverse Hg methylators prevailing in mangrove subsurface sediments can also use versatile electron acceptors (Fig. 2D). First, most MAGs belonging to Zixibacteria, Nitrospirae, Syntrophobacterales, and Desulfobacterales contained the key genes involved into the complete sulfate-reducing pathway (i.e., sat, aprAB, and dsrABD), which has been previously shown to support Hg methylation in mangrove sediments [93]. Second, consistent with the traditional view [4], the hgcA-carrying Methanomicrobia MAGs encoded functional genes responsible for CO2 reduction to CH4. Third, the GWC2-55-46, MBNT15, Nitrospirae, and Desulfobacterales MAGs harbored genes encoding nitrate/nitrite reductase to support Hg methylation potentially. Last, the Nitrospirae and Syntrophobacterales MAGs also contained porin-cytochrome modules (mtrAB), suggesting iron as an alternative terminal electron acceptor in these Hg methylators.
Phylogenetic diversity of Hg methylators in mangrove subsurface sediments
A meta-omic study has recently provided an updated catalog of Hg methylators, the Hg-MATE-Db, using both pure culture and environmental datasets [44]. This has expanded our understanding of the phylogenetic diversity of Hg methylators. Nine microbial groups were identified as having the highest diversity of Hg-methylating microorganisms across different environments [44], including paddy soils, brackish water, lake water, and reservoir and lake sediments. These groups belong to Acidobacteria, Bacteroidetes, Desulfobacterales, Chloroflexi, Euryarchaeota, Firmicutes, Nitrospirae, Spirochaetes, and PVC lineages.
To uncover the phylogenetic diversity of Hg methylators in mangrove sediments, a phylogenetic tree was constructed by integrating our detected hgcA sequences into Hg-MATE-Db (Fig. 3). Results showed that 65 hgcA-carrying MAGs in mangrove sediments were affiliated with seven microbial groups with the highest diversity (excluding Firmicutes and Spirochaetes lineages), as well as two additional groups (GWC2-55-46 and MBNT15 lineages). The number of microbial groups with the highest diversity in mangrove sediments was found to be second only to paddy soils, but higher than in brackish water, lake water, and sediments from reservoirs and lakes (Fig. 3 and Supplementary Table S3). These findings, together with more than 90% of these hgcA-carrying MAGs with a high occurrence frequency below 20 cm (Fig. 2B), point to a high phylogenetic diversity of Hg methylators in mangrove subsurface sediments.
Maximum-likelihood tree (IQ-TREE 2, LG + F + I + I + R9) was constructed using an alignment of HgcA proteins from hgcA-carrying MAGs and Hg-MATE-Db. The representative hgcA sequences recovered in this study are indicated by stars. The Hg-methylating lineages are labeled with colored circles, indicating their distribution in different environments. The scale bar represents one amino acid substitution per sequence position.
To assess the taxonomic novelty of these Hg methylators, we performed comparative genome analyses of average nucleotide identity (ANI) using our hgcA-carrying MAGs and the reference sequences in the Hg-MATE-Db (Fig. 3 and Supplementary Figs. S8–S17). Our results showed that Nitrospirae MAGs in our dataset and Hg-MATE-Db were largely syntenic and shared up to 99.4% ANI, suggesting the conservation of Hg-methylating Nitrospirae across different environments. On the other hand, our Methanomicrobia, GWC2-55-46, MBNT15, and Acidobacteria MAGs shared around 80% ANI with reference genomes, implying that these MAGs could be previously unreported species with Hg-methylating capacity. Notably, the Anaerolineae, Phycisphaerae, Desulfobacterales, and Bacteroidales MAGs shared less than 70.0% ANI against their reference genomes, indicating that they form distinct subfamilies and differ from previously known ones when following the proposed guidelines [110, 111]. These findings suggest that mangrove subsurface sediment is a promising source for identifying novel microorganisms and broadening the repertoire of Hg methylators.
Candidate phylum Zixibacteria is an understudied Hg methylator in diverse anoxic environments
Based on our phylogenetic analyses against the Hg-MATE-Db, we discovered a surprising finding: the Zixibacteria lineage, located between Bacteroidales and Phycisphaerae in the phylogenic tree (Fig. 3), has not been included in the updated catalog of Hg methylators. However, we found that the relative abundance of Zixibacteria MAGs ranked among the top 9 Hg-methylator taxa across all depths, and even ranked 5th at 30–50 cm depths (up to 9.1%, Supplementary Fig. S7) by mapping all the metagenomic reads against the genomic assemblies. Previous study has also identified subsurface-dwelling Zixibacteria lineage as a crucial hub member in mangrove sediment microbial communities, not only interacting frequently with Deltaproteobacteria but also providing metabolites for Bacteroidetes and Chloroflexi [112]. The latter three phyla have been often classified as the microbial taxa of putative Hg methylators in mangrove sediments [8]. Especially for the Syntrophobacteraceae members in Deltaproteobacteria, a recent cultivation study showed that they could enhance Hg methylation by interacting with their syntrophic partners [24]. Thus, we endeavored to decipher whether and how the “hub” Zixibacteria members methylate Hg in mangrove sediments, which is crucial for expanding the diversity of Hg methylators and understanding the microbial mechanism of Hg methylation in the subsurface biosphere.
As no cultivated representative of Zixibacteria is currently available for functional assays, we firstly investigated its potential in Hg methylation by comparing the HgcAB protein sequences from six Zixibacteria MAGs (average completeness 86.7%, average GC counts 53.2%; Supplementary data 4) against six confirmed Hg-methylating bacterial strains isolated from anaerobic environments (Supplementary Table S4) [113,114,115,116,117]. Our results showed that the genetic context of HgcAB sequences and the conservation of respective orthologous gene clusters in six Zixibacteria MAGs were all consistent with the confirmed Hg-methylating bacterial strains (Fig. 4A). The alignment of the HgcA protein sequences has revealed that HgcA is a cobalamin-dependent protein with a highly conserved motif in the CfsA subunit of CFeSP. Aside from HgcA, HgcB is also required for methylation activity [1], and its protein sequence alignment revealed the presence of two CXXCXXCXXXCP motifs, which are known to bind [4Fe-4S] clusters [1].
A Multiple sequence alignments of 12 HgcA and HgcB orthologs from candidate phylum Zixibacteria and six confirmed Hg-methylating bacterial strains. Red boxes indicate highly conserved regions including the putative cap helix [consensus sequence motif, N(V/I)WCA(A/G)GK] in HgcA and two strictly conserved CX2CX2CX3C motifs characteristic of [4Fe-4S] clusters in HgcB. Abbreviations for amino acid residues: A, Ala; C, Cys; D, Asp; E, Glu; F, Phe; G, Gly; H, His; I, Ile; K, Lys; L, Leu; M, Met; N, Asn; P, Pro; Q, Gln; R, Arg; S, Ser; T, Thr; V, Val; W, Trp; Y, Tyr; X, any amino acid. B AlphaFold2 model of Zixibacteria HgcA protein, with an enlarged view of the functional domain bound to cobalamin. The ionic interaction between cysteine and cobalt is represented through a dotted box. C AlphaFold2 model of Zixibacteria HgcB protein in complex with two [4Fe4S] clusters.
We constructed AlphaFold2 models of the putative HgcA and HgcB proteins in Zixibacteria. Our analysis of the putative HgcA protein structure in Zixibacteria revealed a globular domain at the N-terminus and five transmembrane spanning helices at the C-terminus (Fig. 4B). By superposing the cobalamin onto the HgcA AlphaFold2 model in Zixibacteria, we found a strong H-bonding network between HgcA protein and cobalamin (Fig. 4B), which was thought to be bound with HgcA protein for catalytic methylation activity [6]. Importantly, the HgcA AlphaFold2 models in Zixibacteria and six confirmed Hg-methylating bacteria were aligned with a striking level of structural concordance (Supplementary Fig. S18). In addition to the HgcA, our analysis of the HgcB protein structure in Zixibacteria also revealed the binding of its protein with 2[4Fe-4S] clusters at the N-terminal by molecular stacking (Fig. 4C), showing a high similarity with the functionally validated Hg-methylating strain [1]. Together, our results provide genomic and structural support for the functional potential of Zixibacteria lineage in Hg methylation.
In an effort to understand the extent of Hg-methylating capacity in Zixibacteria lineage, we conducted a broader search of the NCBI SRA metagenomic datasets and obtained 135 medium- and high-quality Zixibacteria MAGs belonging to eight described orders (Fig. 5, Supplementary Data 4, 5). Among all these currently available Zixibacteria MAGs, 18 MAGs were identified to carry hgcA with conserved motifs and were distributed within six different orders, with the exception of DG-27 and MSB-5A5 orders. Our functional annotation showed that more than 70% of hgcA-carrying MAGs contained genes encoding sulfate reduction, while DG-27 and MSB-5A5 orders contained no MAGs with the potential for sulfate reduction. These results suggest that the Zixibacteria members primarily perform sulfate reduction to support Hg methylation in anoxic environments, as hgcA-carrying MAGs in Zixibacteria did not possess the capability to utilize other commonly reported electron acceptors such as iron, NO3−, NO2−, or CO2 (Supplementary Data 6).
The maximum-likelihood tree was constructed using an alignment of maker genes from Zixibacteria genomes obtained in this study and from the NCBI SRA database. The genomes with hgcA sequences are highlighted in red. The diverse sources of Zixibacteria genomes are depicted in the outer circle with different colors representing the distinct environments. The capacity of sulfate reduction is indicated by brown circles, and the scale bar indicates substitutions per site. The tree was constructed using 1000 ultrafast bootstrap replicates, with black dots representing values greater than 10% at the nodes.
The environmental distribution of Zixibacteria was thoroughly investigated by screening homologous 16 S rRNA gene sequences against all available metagenomes within the IMG/M database. This search generated a total of 710 unique occurrences. The categorization of sampling locations revealed a broad environmental distribution of the Zixibacteria lineage, spanning 17 different types of environments and including natural, host-associated, and engineered systems (Supplementary Fig. S19). Of them, the hgcA-carrying Zixibacteria were mainly recovered from natural environments, including mangrove sediment, marine sediments, cold seep, hydrothermal vent, estuarine sediment, groundwater, microbial mat, and deep mine (Fig. 5). In addition to mangrove sediment (this study and Zhang et al.) [8], the Zixibacteria lineage in other environments have not been reported to carry hgcA gene, highlighting its currently understudied roles in Hg methylation across most anoxic environments. However, no clear environment-specific distribution was observed across the six different orders (Fig. 5), suggesting that hgcA-carrying Zixibacteria are well adapted to diverse anoxic environments.
Unique metabolic and lifestyle features of candidate phylum Zixibacteria in Hg methylation
During Hg methylation, the methyl group (as CH3+) in MeHg is widely thought to be transferred to HgcA by a cobalamin-dependent methyltransferase (i.e., CFeSP methyltransferase from the Wood-Ljungdahl pathway) [1, 118, 119]. However, we found all Zixibacteria MAGs lacking genes for the entire Wood-Ljungdahl pathway (Fig. 6A and Supplementary Data 7). To investigate the putative CH3+ transfer within the Zixibacteria lineage, we conducted a search for typical cobalamin-dependent methyltransferases [72, 119,120,121,122,123,124,125,126,127]. Our results revealed that the Hg-methylating Zixibacteria lineage only encoded the cobalamin-dependent methyltransferase involved in the methionine synthesis pathway, rather than the Wood-Ljungdahl pathway, coenzyme M pathway, or one carbon folate cycle (Supplementary Data 8). In particular, all six different orders in Zixibacteria consistently encoded the MetH (Fig. 6B), a cobalamin-dependent methyltransferase catalyzing CH3+ transfer from 5-methyl-tetrahydrofolate to homocysteine during the methionine synthesis pathway [128]. Considering that CH3+ in MeHg is traditionally obtained from 5-methyl-tetrahydrofolate by cobalamin-dependent methyltransferase [1, 118, 119], we hypothesized that the MetH in Hg-methylating Zixibacteria may have the potential to transfer CH3+ from 5-methyl-tetrahydrofolate to HgcA, forming CH3-cob(III)alamin-HgcA (Fig. 6B). To validate this hypothesis, we performed docking predictions for HgcA and MetH. For the most reliable docking model of HgcA and MetH (confidence score: 0.93), the interaction interface was defined as those amino acid residues within 5.1 Å of the surface volume encapsulating the atom positions (Fig. 6C and Supplementary Data 9). The docking score of −280.66 kcal/mol indicates a strong interaction of HgcA and MetH, which is mediated by 26 hydrogen bonds, 7 hydrophobic interactions, and 1 salt bridge (Fig. 6C and Supplementary Fig. S20). Among them, Arg in MetH exhibited the highest binding activity to HgcA, and the closest atom distance between MetH and HgcA was 1.4 Å in a hydrogen bond between Arg156 and Ser83 (Supplementary Data 9). These observations suggest that hydrogen bonding interactions associated with Arg side chains are likely important for HgcA and MetH protein binding.
A Metagenomics-based mechanism reconstruction of Hg-methylating Zixibacteria. Abbreviation: Fd, ferredoxin. The blue balls with “e” indicate electrons. B Potential mechanism as used for the methionine biosynthesis pathway to donate methyl groups to Hg methylation across all orders within Hg-methylating Zixibacteria. The dashed arrows signify proposed steps that require experimental validation. C Schematic representation of HgcA-MetH interaction in the 3D structure. All of the top 5 candidate models show confidence scores greater than 0.9. The visualization was accomplished using PyMOL.
Furthermore, Zixibacteria MAGs encoded the benzoyl-CoA reductase, which was reported to be structurally similar to the CODH/acetyl-CoA synthases in Wood-Ljungdahl pathway and showed a strong correlation with HgT and MeHg levels [129, 130]. The benzoyl-CoA reduction enables the Hg methylators to adapt to environments contaminated by hydrocarbons [131], which is of particular significance as hydrocarbons are frequently found in mangrove sediments, as a result of industrial discharge, shipping activities, and urban runoff [132]. These findings indicate that the Zixibacteria lineage has the potential to play a role in Hg methylation, with the added benefit of removing hydrocarbons from contaminated mangrove sediments.
Energy limitation of mangrove subsurface sediments has been widely acknowledged [14,15,16], but the provision of energy is essential for Hg methylation [109]. To address this issue, the researchers have proposed syntrophic microbial interaction as a major energy source for Hg methylation in energy-limited environments [24]. Through an examination of microbial syntrophic features, our study found evidence for the presence of these partnerships in mangrove subsurface sediments. The Zixibacteria MAGs, which carry the hgcA gene and have the sulfate-reducing capacity, were found to contain genes involved in interspecies electron transfer (including multi-heme cytochromes, conductive filaments and molecular hydrogen; Supplementary Data 10), the significant features for a syntrophic lifestyle [22, 133]. Also, we here recovered several Syntrophobacterales MAGs, which not only possessed genes involved in interspecies electron transfer (Supplementary Data 10), but also were dominant in mangrove subsurface sediments (up to 13.1% in the Hg-methylating genomic assemblies, Supplementary Fig. S7). In addition, our 16 S rRNA gene-based microbial network inference and metagenomic-based linear regression revealed the positive co-occurrence between the Syntrophobacterales and Zixibacteria lineages (Supplementary Figs. 21, 22), further supporting the potential for syntrophic partnerships in the mangrove subsurface sediments.
Further evidence for their potential syntrophic relationship was the absence of cobalamin biosynthetic genes in the Hg-methylating Zixibacteria. This was inconsistent with our AlphaFold2 model (Fig. 4A; indicating the participation of cobalamin in Hg methylation) and metabolic predictions (Fig. 6; with cobalamin-dependent methyltransferase involved in methionine synthase). However, the syntrophic Syntrophobacterales lineage contained all relevant genes for cobalamin biosynthesis (Supplementary Data 11), potentially leading to its interactions with Hg-methylating Zixibacteria [134]. Finally, both the co-occurrence pattern and the paired genomic analyses suggest that the Zixibacteria lineage has a unique metabolic strategy for Hg methylation in mangrove subsurface sediments, dependent on a syntrophic relationship with another organism in energy-limited conditions. Future cultivation experiments with Hg methylators isolated from mangrove sediments will help to verify this hypothesis and further clarify the nature of their syntrophic interactions.
Responses