Mercury methylating microbial communities of boreal forest soils

The formation of the potent neurotoxic methylmercury (MeHg) is a microbially mediated process that has raised much concern because MeHg poses threats to wildlife and human health. Since boreal forest soils can be a source of MeHg in aquatic networks, it is crucial to understand the biogeochemical processes involved in the formation of this pollutant. High-throughput sequencing of 16S rRNA and the mercury methyltransferase, hgcA, combined with geochemical characterisation of soils, were used to determine the microbial populations contributing to MeHg formation in forest soils across Sweden. The hgcA sequences obtained were distributed among diverse clades, including Proteobacteria, Firmicutes, and Methanomicrobia, with Deltaproteobacteria, particularly Geobacteraceae, dominating the libraries across all soils examined. Our results also suggest that MeHg formation is also linked to the composition of non-mercury methylating bacterial communities, likely providing growth substrate (e.g. acetate) for the hgcA-carrying microorganisms responsible for the actual methylation process. While previous research focused on mercury methylating microbial communities of wetlands, this study provides some first insights into the diversity of mercury methylating microorganisms in boreal forest soils.

Mercury (Hg) is a potent toxin that might cause severe negative effects on wildlife and human health 1 . The toxicity of Hg is of such concern that 128 countries have signed the Minamata Convention, a global treaty that entered into force in August 2017 with the explicit objective to reduce Hg emissions and protect human health and the environment. High Hg emissions in the past have led to high present-day Hg levels in different parts of the atmosphere, oceans and terrestrial ecosystems 2,3 . Fish consumption is the main pathway for human Hg exposure 4 . High Hg levels in freshwater fish, not seldom at concentrations unsafe for human consumption, have in recent decades raised much concern in many boreal regions [5][6][7] . In these freshwater boreal catchments, the link between amount of atmospheric Hg deposited and the Hg accumulated in food-webs is complex 8 . Attention has therefore been directed to understand processes within the catchment soils that may contribute to the formation and mobilisation of readily bioavailable methylmercury (MeHg) in runoff water. Local biogeochemistry in catchment soils may for example influence the magnitude and timing of the response in fish Hg concentrations following decreased deposition loadings. One example of this is seen in a 40-year monitoring dataset from Canada showing flat or increasing trends of Hg in freshwater fish up to 2012 although Hg deposition has decreased during recent decades 7 .
Because Hg has a strong affinity for reduced sulphur or thiol (RSH) functional groups of soil organic matter (OM) 9,10 , the increased atmospheric deposition of Hg during the industrialisation period has resulted in high Hg concentrations in organic-rich soils 11 . As a consequence, the OM-rich soils, characteristic of the boreal biome, have retained Hg deposition from both natural and anthropogenic emissions, and now represent an important global Hg stock 9,12 . For example, recent atmospheric deposition has increased the stock of Hg in the organic-rich upper layer of Swedish forest soils by a factor of three to four 11,13 . This is of special concern because soil OM has been identified as an important vector of Hg and methylmercury (MeHg) transport from catchments to surface waters in boreal areas 14,15 . Indeed, the mobilisation of inorganic Hg (Hg (II) ) and, the more harmful, MeHg from soils by means of OM-mediated transport has been linked to MeHg accumulation in lake sediments within catchments 15,16 and in fish 17 . The role of forest soils is thus evident from the increase in MeHg export from forests to SCIENtIFIC REPORTS | (2019) 9:518 | DOI: 10.1038/s41598-018-37383-z aquatic ecosystems 16,18,19 , and the subsequent bioaccumulation in downstream fish 17,20 . Since forest soils are an important site for MeHg formation 21 , it is crucial to understand the processes and the organisms involved in MeHg formation in boreal soils. The methylation of Hg (II) to MeHg is biologically mediated 22 and takes place under oxygen deficient conditions such as those of flooded soils, sediments 15 , anoxic water columns 23 and suspended particles of aquatic systems 24,25 . Typical forest environments with high MeHg formation are wetlands 26 , organic-rich riparian soils 27 , and soils that have been recently become water logged soils after forest harvest 21 . Specific strains of sulphate-reducing bacteria 28,29 , iron reducing bacteria (FeRB) 30,31 , methanogens 32 and Firmicutes 33 have the capability to methylate Hg (II) . However, a number of factors controlling microbial activity and/or the geochemical speciation of inorganic Hg (II) will govern MeHg formation in the environment 21,[34][35][36] . For example, increases in temperature might lead to increases in biological activity and subsequently also higher Hg (II) methylation rates 37 . Redox potential also seems to be a key factor as suboxic and mildly reducing conditions seem to promote high Hg (II) methylation rates, whereas anoxic and strongly reducing conditions might lead to elevated sulphide concentrations that eventually prevent Hg (II) from being available for methylation 37 . Sulphur plays a major role in influencing Hg (II) methylation by directly affecting the activity of some methylating bacteria (e.g. sulphate reducing bacteria, SRB) and/or control the availability of Hg (II) for methylation 10 . Specific organic matter (OM) compounds can promote Hg (II) methylation by enhancing bacterial activity 15 , but also by defining Hg (II) speciation 38 and Hg (II) availability 39,40 . OM can also facilitate Hg (II) methylation by inhibiting mercury sulphide (HgS(s)) precipitation or enhance HgS(s) dissolution thereby providing available Hg (II) for methylating microorganisms 41 . High OM concentrations might also decrease Hg methylation by formation of high mass molecular mass complexes that hamper Hg (II) availability 39 . Recently it has been concluded that the availability of Hg (II) depends heavily on the S (−II) concentration in porewater and the RSH(aq)/RSH(ads) molar ratio of dissolved OM 38 . Besides all the geochemical factors that might directly or indirectly affect Hg (II) availability and methylation, a recent study suggest that the composition of the combined bacterial community may also influence the structure of Hg (II) methylating communities 42 . Together these studies, mainly performed in aquatic systems, highlight the importance of geochemical conditions for determining the availability of Hg (II) and the activity and composition of the microbial communities involved, directly or indirectly, in MeHg formation.
The identification of two functional genes, hgcA and hgcB, which play essential roles in Hg (II) methylation 22 , provided the means to more directly characterise the complexity of microbial communities involved in the formation of MeHg in natural ecosystems. This approach has been applied to marshes, sediments and swamps in several geographic regions [42][43][44][45][46] ; rice paddies in China 47 , and water conservation areas of the northern Everglades, USA 48 . However, very little work to date has been conducted to reveal the distribution of microbial groups responsible for Hg (II) methylation in forest soils within the vast boreal biome. To the best of our knowledge, no studies have directly described the composition and the spatial variation in Hg (II) methylating microbial communities in such forests. Therefore, the primary goal of this paper was to describe Hg (II) methylating microbial communities in various boreal forest soils and identify characteristics important for shaping these communities. High-throughput next generation sequencing of amplified 16S rRNA and hgcA genes combined with molecular barcoding and detailed soil geochemical characterisations were performed to study the Hg (II) methylating microbial communities in 200 soil samples from three different boreal forest regions ( Fig. 1) in order to shed light on the biogeography of microorganisms responsible for MeHg formation in the boreal landscape.

Results
Bacterial community composition in boreal forest soils. We collected 200 boreal forest soil samples distributed across eight catchments in Sweden in 2012 (Tables S1 and S2). A total of 3 321 197 high quality 16S rRNA sequences remained after quality control and chimera removal (7-72 911 reads per sample). The sample with only 7 reads was removed, and we then rarefied the rest of the data to the remaining sample with the fewest reads (1692 reads). The final rarefied sequence dataset (329 940 reads) clustered into 33 158 operational taxonomic units (OTUs) using a similarity threshold of 97%. In the rarefied dataset, 35 taxa at phyla level, 69 taxa at class level, 119 taxa at order level, and 187 taxa at family level were detected from all the soil samples across three regions. The overall coverage of the forest bacterial community is reflected in the combined richness detected for random subsets of analysed samples. The logarithmic shape indicated that most of the considerable OTU richness occurring in the forest soils was accounted for in the combined dataset (Fig. S1). Among the dominant phyla across all regions (>5% relative abundance), Acidobacteria was the most abundant, followed by Proteobacteria, Planctomycetes, Bacteroidetes, Parcubacteria and Verrucomicrobia (Table 1). Combined, these phyla accounted for 77.5% of the total sequences ( Table 1). Most of the previously identified clades known to contain Hg (II) methylators 33,49 were detected in the present study, including Deltaproteobacteria (3.31% of the total reads), Chloroflexi (2.60% of the total reads), Firmicutes (0.77% of the total reads) and Euryarchaeota (0.66% of the total reads) ( Table 1). Microbial community composition based on 16S rRNA sequences in the 34 studied MeHg hotspots (%MeHg >1%) showed a similar pattern in terms of the dominant phyla (>5% relative abundance), with Acidobacteria and Proteobacteria being the most abundant ones. However, Bacteroidetes and Chloroflexi contributed much more to the total communities at these hotspots compared to the combined dataset across all 200 samples (Table 1).
A non-metric multidimensional scaling (nMDS) plot based on 16S sequences was used to visualise the composition of the bacterial community among samples. Unclassified Acidobacteriales, Unclassified Ignavibacteriales, Spirochaetaceae, Holophagaceae, Anaerolineaceae, Betaproteobacteria and Tepisiphaeraceae were important contributing families for shaping the differences in bacterial community composition among samples (Fig. 2). Geochemical factors that were correlated (correlation coefficients > 0.5) with the bacterial composition were projected on top with longer vectors implying stronger correlations (Fig. 2). %MeHg, reflected by bubble sizes, presented a strong coupling to the bacterial community composition, which was further confirmed by %MeHg SCIENtIFIC REPORTS | (2019) 9:518 | DOI:10.1038/s41598-018-37383-z presenting a long vector among all the geochemical factors (Fig. 2). Water content, C%, S% and N% were all found to be the factors that affected the composition of soil bacterial community (Fig. 2), indicating that a supply of organic matter and nutrients in the moist soil shapes the bacterial community. This is in agreement with previous research that pointed out the contribution of nutrients and organic matter to bacterial activities and Hg (II) methylation 15,37 . Also, S was well correlated with both C and N (Table S3), suggesting that most of the measured sulphur in the sampled soils has likely an organic origin. This has been found as a common feature in boreal soils 27,36,50 . Unclassified Fibrobacterales, Methanosaetaceae, unclassified Ignavibacteriales, Spirochaetaceae, Holophagaceae and Anaerolineaceae exhibited the highest correlations with %MeHg (Table 2). Syntrophobacteraceae, Methanosarcinaceae, Methanoregulaceae, Desulfobulbaceae, Syntrophaceae, Desulfobacteraceae and Dehalococcoidaceae, were also found relevant to the bacterial community composition in high-%MeHg sites ( Table 2).

Distribution of Hg (II) methylators. Among the screened 200 soils samples, we selected those with high
MeHg concentrations and %MeHg (>1%), and defined them as "MeHg hotspots" (see "MeHg hotspots" soils geochemistry descriptors in Table S4, n = 34). In 34 of these MeHg hotspots, the relative abundance of microbial families carrying representatives known to methylate Hg (II) was assessed based on hgcA sequences 33,49 . A total of 1 257 577 hgcA sequences remained after quality control and chimera removal (11 404-55 461 reads per sample). The hgcA dataset was rarefied to the remaining sample with the fewest reads (11 404 reads). The rarefied sequence dataset accounted a total of 387 736 reads that clustered into 573 operational taxonomic units (OTUs) using a similarity threshold of 97%. As for the 16 rRNA, the logarithmic shape indicated that most of the considerable species richness of Hg (II) methylators occurring in the forest soils was accounted for in the combined dataset ( Fig. S1). Representative sequences from 22 families were found in the 34 analysed MeHg hotspots. Of all the hgcA sequences, 3.13% were not taxonomically assigned (unclassified), 0.28% were unclassified Euryarchaeota, and 7.28% could not be assigned beyond the rank of Bacteria (Unclassified Bacteria). The majority of the sequences annotated to the level of family clustered with Deltaproteobacteria, making up 85.4% of all the hgcA reads (     (Table 3). Unclassified Desulfuromonadales, Geobacteraceae, Ruminococcaceae, unclassified Desulfovibrionales, Desulfovibrionaceae, and unclassified Deltaproteobacteria seemed to contribute to differences in the composition of Hg (II) methylators in the studied soils (Fig. 3a). Among the measured geochemical parameters, the S% and the C/S seemed to have an impact on shaping the community composition of Hg (II) methylators (Fig. 3b). Moreover, Methanoregulaceae, Desulfovibrionaceae, Desulfuromonadaceae, Desulfarculaceae and Methanomassiliicoccaceae correlated positively with S% and negatively with C/S (Table S5). In the studied MeHg hotspots, S was strongly correlated with both C and N (Table S6), suggesting most of the measured sulphur in the hotspots is also likely presented in organic forms.
Phylogenetic analysis of hgcA genes. All the Proteobacteria families belonged to Deltaproteobacteria, a class with which most currently confirmed Hg (II) -methylating bacteria are affiliated 51,52 . When combined, the 20 most abundant OTUs accounted for 72% of the total reads. Noteworthy, phylogenetic analysis revealed that the most abundant Hg (II) -methylating OTUs ("OTU_0005", "OTU_0705", "OTU_0008", and "OTU_0012") in the studied forest soils were either taxonomically assigned as Geobacter sp. or phylogenetically related to Geobacter species (Fig. 4). Among the 20 most abundant OTUs, 17 were taxonomically annotated as Deltaproteobacteria. Among these 17 OTUs, 9 were taxonomically annotated as Geobacter and 8 were phylogenetically related to Geobacter species (Fig. 4). Summing the 9 OTUs taxonomically annotated as Geobacter with the 8 OTUs phylogenetically related to Geobacter species, the resulting 17 OTUs accounted for 62% of the total hgcA reads (Fig. 4). While the 5 th most abundant OTU was taxonomically denoted as Firmicutes (Ethanoligenens), the 6 th and 7 th most abundant OTUs could not be annotated beyond the bacterial domain.

Discussion
Community composition of Hg (II) methylators in boreal forest soils. Among the diverse microbial communities seen in the soil samples (Table 1), most of the previously identified Hg (II) methylating groups, e.g., Deltaproteobacteria, Chloroflexi, Firmicutes and Euryarchaeota could be detected (Table 3). Deltaproteobacteria have been considered a predominant Hg (II) methylating class in anaerobic soils 44,47,48 . In the present study, Deltaproteobacteria were also the predominant Hg (II) methylators at the hotspots with Geobacteraceae as the most represented family. When considering the OTUs that were taxonomically annotated, this family alone contributed over 30% of all hgcA reads. However, if we also account for the OTUs that were phylogenetically related to Geobacteraceae, this family might have contributed to more than 60% of the hgcA reads. The importance of Geobacteraceae could be seen at all the sampled sites and particularly in Strömsjöliden (Table 3). Iron reducing bacteria (FeRB) have previously been shown to be important for Hg (II) methylation in some environments 30,31,46,51 , and most Geobacter tested so far are particularly efficient at MeHg formation in the laboratory 31 . This suggests that the ability to methylate Hg (II) , a typical feature among the Geobacteraceae, is present in the studied soils and widely distributed in terrestrial and aquatic ecosystems. However, while previous studies have quantified the contribution of SRB (i.e molybdate inhibitor) and methanogens (i.e bromoethanesulfonate inhibitor) to MeHg While SRB are considered to be the principal Hg (II) methylators in aquatic systems [55][56][57][58][59] , not much information is available on Hg (II) methylators in soils. However, identified SRB in the hotspots only accounted for a minor portion of Hg (II) methylators (Table 3). However, it is nevertheless plausible that at least some of the hgcA sequences annotated as unclassified Deltaproteobacteria (Table 3) could be unknown Hg (II) methylating SRB or even Hg (II) methylating sulphate-reducing syntrophs, capable of syntrophic fermentation of simple organic acids in the absence of sulphate as the terminal electron acceptor 60,61 . Therefore, we cannot discard the possibility that also SRB contribute significantly to Hg (II) methylation in the studied systems. A previous study based on selective inhibitors and rate measurements indeed suggested SRB played an important role in MeHg formation in boreal forest soils 36 . Additionally it has been demonstrated that even when SRB belong to the 'rare biosphere' of peatlands, they contribute significantly to respiration processes 62 . Ruminococcaceae belongs to another newly confirmed representative of Hg (II) methylators, the Firmicutes 33 . Firmicutes contributed to Hg (II) methylating microbial communities at the water conservation areas of the Florida Everglades 48 but were not detected in boreal wetlands 44 . In the present study, Ruminococcaceae were prominent contributors to the hgcA pool in hotspots from Örebro and in all soils from Strömsjöliden (Table 3). They could thus play a role in shaping the composition of Hg (II) methylating community as further indicated by the negative correlation though weak between Ruminococcaceae and C/S, a primary geochemical factor shaping Hg (II) methylating communities in the hotspots (Table S5 and Fig. 3b). Not much research has been devoted to the possible relationship between organic S and Hg (II) methylating Ruminococcaceae. Considering the abundance of this group in forest soils, further efforts are needed to shed light on the metabolic or physiological pathways of Hg (II) methylating Ruminococcaceae.
Methanogens were early on suspected to be responsible for Hg (II) methylation 63 , but not until recently were they verified as a significant source of Hg (II) methylators in various environments 32,44 . In the hotspots in the studied soils, they were also detected, though they were not very abundant in the Hg (II) methylating microbial community. Chloroflexi has recently been identified as potential Hg (II) methylators in the water conservation areas, paddy soils and wetlands 44,48,64 . The hgcA data did not confirm any significant role of this group in MeHg production in boreal forest soils (Table 3), even though 16S rRNA data revealed non-Hg (II) methylating Chloroflexi (e.g. the class Anaerolineae) in soils from all three regions (Table 1).
Previous studies have mainly explored flooded environments such as paddy soils 47 , boreal wetlands 44 and the water areas of the Florida Everglades 48 . Hence our study provided important new information on the composition and diversity of Hg (II) methylating microbial communities in non flooded boreal forest soils and the boreal landscape, and in doing so identified Geobacteraceae as significant Hg (II) methylators in the terrestrial biome. The diversity of Hg (II) methylators described in this study need to be interpreted cautiously, as bias is inherent in methods employing PCR amplification of any variable target gene. The hgcA gene was only recently discovered and the optimization of the appropriate methods and, in particular the design of primers for the hgcA amplification, is still ongoing 65 . Additionally, DNA based methods only reveal the presence of organisms, while alternative approaches based on transcription data, proteomes or rate measurements are needed for verifying their activity. Our data nevertheless provide new insights about Hg (II) methylating microbial communities in boreal forest soils and can as such guide and serve as a resource for future research efforts in this field.

Interplay between bacterial communities and Hg (II) methylators. %MeHg has previously been used
as a proxy for methylation efficiency 66,67 , and high %MeHg has also in a few cases been shown to correlate positively with the abundance of Hg (II) methylators 21,68 . In the current study, sites with high %MeHg featured bacterial communities different from those observed at sites with low % MeHg (Fig. 2). Although, families known to contain Hg (II) methylators (Syntrophobacteraceae, Methanosarcinaceae, Methanoregulaceae, Desulfobulbaceae, Syntrophaceae, Desulfobacteraceae and Dehalococcoidaceae; 25) were found at sites with high %MeHg, there were also positive correlations between %MeHg and families that are not known to host Hg (II) methylators, such us unclassified Fibrobacterales, Methanothrix (formerly Methanosaeta), unclassified Ignavibacteriales, Spirochaetaceae, Holophagaceae and Anaerolineaceae (Table 2). This suggests that not only the Hg (II) methylators themselves, but also the supporting and interacting bacterial communities residing in the soil environment may influence MeHg formation across the studied regions. Anaerolineaceae, Spirochaetaceae and Holophagaceae are for example known to generate acetate by fermentation processes 69 . Fibrobacterales, have recently been suggested to have an important role in cellulose hydrolysis in anaerobic environments, including soils 70 . The Ignavibacteria class was recently described (Iino et al., 2010) and the physiology and metabolic capacities of this group is still poorly known, even if a distinctive feature of this group is the ability to grow on cellulose and its derivatives with the utilization of Fe(III) oxide as electron acceptor 71 . It may well be that these families, which correlated well with %MeHg (Table 2) and seem to be involved in the degradation of long chain OM compounds 72,73 , promoted MeHg production by providing appropriate substrates (e.g. acetate) for the Hg (II) methylators. Hg (II) methylators and non-Hg (II) methylating members of Desulfobulbaceae, known to oxidise organic substrates incompletely to acetate 74 , might also have provided the necessary substrate to Hg (II) methylators (Table 2). Based on our results, we propose an important role of also the non-Hg (II) methylating bacterial heterotrophs in sustaining the activity of the Hg (II) methylating microorganisms and thereby influencing MeHg formation in boreal forest soils. Moreover, the correlation between Methanothrix and %MeHg deserves special attention. It has been shown that Methanothrix can establish syntrophic cooperation with Anaerolineaceae 72 or Geobacteraceae 75 in methanogenic degradation of long chain carbon compounds (alkanes). As our results show that Geobacteraceae are major contributors to the Hg (II) methylating microbial community (Table 3), the high correlation found between Methanothrix and %MeHg could be the result of the interaction between the non-Hg (II) methylating Methanothrix and the Hg (II) methylating Geobacteraceae. In brief, we provide novel system-level information on putative trophic interactions between non-Hg (II) methylating and the Hg (II) methylating taxa. We further suggest that more in depth studies with metagenome-level sequencing and metabolic pathway reconstruction will be a logical next step to gain a more complete understanding of how Hg (II) methylating bacterial and archaeal species interact in soils.

Conclusions
A newly developed strategy that combine high-throughput hgcA amplicon sequencing with molecular barcoding revealed diverse clades of Hg (II) methylators in forest soils. This study confirms a predominant role of Deltaproteobacteria, and in particular Geobacteraceae, as key Hg (II) methylators in boreal forest soils. Firmicutes, and in particular Ruminococcaceae, were also abundant members of the Hg (II) methylating microbial community. Besides the identified Hg (II) methylators, we suggest that the non-Hg (II) -methylating bacterial community (e.g. Anaerolineaceae, Holophagaceae and Spirochaetaceae) might have contributed to the net MeHg formation SCIENtIFIC REPORTS | (2019) 9:518 | DOI:10.1038/s41598-018-37383-z (%MeHg) by processing OM and thereby providing low molecular mass OM compounds as a substrate to Hg (II) methylators (e.g acetate). By revealing linkages between Hg (II) methylators and non-Hg (II) methylators, our results call for further community-level work on the metabolic interactions in soil microbial communities to understand Hg (II) methylation. Such studies would need to go beyond the Hg (II) methylating microbial populations. Our findings provide a better understanding of Hg (II) methylating microbial communities in forest soils and the boreal landscape.

Materials and Methods
Site description. Soil samples were collected from 200 sites in October 2012 and were distributed across eight catchments in three boreal forest regions in Sweden (Table S1 and Table S4. The daily mean air temperatures during the 9 sampling days in September in 2012 varied between 7 and 12 °C in Örebro catchments and 4 and 11 °C in Balsjö and Strömsjöliden catchments. There were no major rain events during the sampling period and the temperature and precipitation was normal for the time of the year. Soil sampling. Soil samples were collected with a soil coring tube (Ø = 23 mm). In each catchment, around half of the samples (n = 12) were collected systematically along the topographic fall line of the hill slope, at set distances from the stream draining the area. These samples were collected from the upper 6 cm of the O horizons or the whole O horizons if these were less than 6 cm deep. The locations of the remaining sampling sites (n = 13) were chosen by actively looking for potential hot spots for MeHg formation, such as wet patches, driving tracks and stump holes. These targeted samples were also collected from various depths, e.g. depths where groundwater levels were most frequently fluctuating were of special interest for potential Hg (II) methylation.
Single-use plastic gloves were used and soil samples for chemical analyses were collected in plastic bags or acid washed Falcon tubes and stored on ice in a cooler during transport to the laboratory (within 8 hours). Soil samples for molecular analyses were collected following adequate aseptic sampling protocols. All sampling equipment was sterilized by washing in 70% ethanol in between samples. Samples were collected in sterilized plastic tubes and frozen in liquid nitrogen directly in the field, and then stored at −80 °C until further processing and analyses.
Chemical analyses. Soil samples were analysed for total Hg (THg), MeHg, water content, and mass percentage of carbon (C), nitrogen (N) and sulphur (S). Samples were freeze-dried and ground by hand in a mortar prior to analyses for THg, C%, N% and S%. Wet and dry weights were measured to estimate the water content. Total Hg was measured using a Perkin Elmer SMS100 total Hg analyser in accordance with US EPA method 7473. The method includes a thermal decomposition step, followed by amalgamation and atomic absorption spectrophotometric detection (working range 0.05-600 ng). Reproducibility and accuracy of measurements were checked by analyses of replicate samples and reference standards. Analyses of MeHg were done by using GC-ICPMS 77 on fresh samples immediately after thawing. C, N and S were analysed on dry soils packed tightly in tin capsules (Elemental Microanalysis, 6.4 mm) and subsequently measured by high temperature catalytic oxidation with a COTECH ECS 4010 elemental analyser calibrated with sulfanilamide standard (C 41.84%, N 16.27%, H 4.68%, O 18.58%, S 18.62%). Analytical precision was <±0.3% for C, ±1.5% for N and ±3.5% for S. Microbiological analyses. 16S rRNA gene. Microbial DNA was extracted from soil samples using the Power soil DNA isolation Kit (MoBio Laboratories Inc, CA, USA) and the quality of the extracted DNA was assessed by gel electrophoresis (1% agarose). Bacterial 16 S rRNA genes were amplified in two steps polymerase chain reaction (PCR) according to the protocol in Sinclair et al. (2015). Briefly, non-barcoded primers Bakt_341F and Bakt_805R (Table S7) were used for the 1 st PCR step of 20 cycles. The resulting PCR products were diluted 100 times before being used as template in a 2 nd PCR step of 10 cycles with similar primers carrying sample-specific 7-base DNA barcodes. All PCRs were conducted in 20 μL volume using 1.0 U Q5 high fidelity DNA polymerase (NEB, UK), 0.25 μM primers, 0.2 mM dNTP mix, and 0.4 μg bovine serum albumin. The thermal program consisted of an initial 95 °C denaturation step for 5 min, a cycling program of 95 °C for 40 seconds, 53 °C for 40 seconds, 72 °C for 60 seconds and a final elongation step at 72 °C for 7 minutes. Amplicons from the 2 nd PCR were purified using the Qiagen gel purification kit (Qiagen, Germany) and quantified using a fluorescence-based DNA quantitation kit (PicoGreen, Invitrogen). The final amplicons after two PCR steps were pooled in equal proportions to obtain a similar number of sequencing reads per sample. Amplicon sequencing was carried out following the protocol described in Sinclair et al. (2015) using the MiSeq instrument. Illumina sequencing was performed by the SNP/SEQ SciLifeLab facility hosted by Uppsala University using 300 bp chemistry. Chimera identification and OTU (Operational Taxonomic Unit) clustering by denoising was done using UNOISE (from USEARCH version 9, refs 69,70  , and a final extension at 72 °C for 7 min. Following this initial step, a 2 nd PCR was conducted to add sample-specific molecular barcodes. Reactions were carried out in 20 μL volumes using 1x Q5 reaction buffer, 0.2 mM dNTP mix, 0.1 μM barcode primers, purified 1 st PCR products and 1.0 U Q5 high fidelity DNA polymerase (NEB, UK) for an initial denaturation of 30 s at 98 °C followed by 18 cycles (10 s at 98 °C, 30 s 66 °C and 30 s at 72 °C), and a final extension at 72 °C for 2 min. The quality and size of the hgcA amplicons were assessed by gel electrophoresis and GelRed visualization on a 1% agarose gel (Invitrogen, USA) prior to purification by Agencourt AMPure XP (Beckman Coulter, USA) after both PCR steps. Quantifications of purified amplicons from the 2 nd stag PCR were performed using the PicoGreen kit (Invitrogen). Amplicons were sequenced using the same method as for the 16S rRNA gene. Forward read sequences were only used in data analysis due to long PCR product. Low quality sequences were filtered and trimmed using SICKLE 79 and adapter were removed by using CUTADAPT 80 . Subsequent processing of reads were performed by USEARCH and clustered at 60% identity cutoff using cd-hit-est 81 . HMMER 82 search was used for taxonomical annotation with manually curated database of Proteobacteria and sequences of Podar et al. (2015) ref. 49 . More details can be found in  ref. 46 .
Phylogenetic analysis. A phylogenetic analysis was performed for hgcA sequences representative for the OTUs observed for the 34 hotspots and existing hgcA entries in our curated database. The sequences were adequately curated and taxonomy homogenized using taxtastic (https://github.com/fhcrc/taxtastic) and the R-package taxize 83 . The obtained protein sequences were aligned with MUSCLE 84 (version 3.8.1551). The alignment was trimmed to the size of the amplicon, and a tree was generated using RAxML 85 (version 8.2.4) -with the PROTGAMMLG model and autoMR to choose the number of necessary bootstrap resamplings (n = 750). This tree and the corresponding alignment were used to generate a reference package for PPLACER 86 . The guppy tool of PPLACER was then used to classify the sequences with a likelihood threshold of 0.8.

Statistical analysis.
Family-level microbial community composition in the different samples were compared using non-metric multidimensional scaling (nMDS) based on Bray-Curtis similarities and using the software PRIMER 7 87 . Information on the common set of samples from community composition based on Bray-Curtis similarities and that from geochemical variables based on Euclidean distance was presented in one single ordination. A combined nMDS plot with bubble and vector plots of geochemical factors projected on the same ordination of community composition was constructed to reveal the relationships between community compositions and potentially explanatory geochemical variables 87,88 . Pearson's correlation coefficient (R) was assessed to reveal linear relationships between variables using a significance level of alpha < 0.05.