Genetic diversity in terrestrial subsurface ecosystems impacted by geological degassing

Abstract

Earth's mantle releases 38.7 ± 2.9 Tg/yr CO₂ along with other reduced and oxidized gases to the atmosphere shaping microbial metabolism at volcanic sites across the globe, yet little is known about its impact on microbial life under non-thermal conditions. Here, we perform comparative metagenomics coupled to geochemical measurements of deep subsurface fluids from a cold-water geyser driven by mantle degassing. Key organisms belonging to uncultivated Candidatus Altiarchaeum show a global biogeographic pattern and site-specific adaptations shaped by gene loss and inter-kingdom horizontal gene transfer. Comparison of the geyser community to 16 other publicly available deep subsurface sites demonstrate a conservation of chemolithoautotrophic metabolism across sites. In silico replication measures suggest a linear relationship of bacterial replication with ecosystems depth with the exception of impacted sites, which show near surface characteristics. Our results suggest that subsurface ecosystems affected by geological degassing are hotspots for microbial life in the deep biosphere.

Fig. 1. Metagenomic and microscopic characterization of the community in subsurface fluids discharged by Geyser Andernach.

A RpS3-based phylogenetic diversity of the organisms in the Geyser Andernach. Centroid rpS3 sequences (after clustering at 99% similarity using cdhit) were used for the calculation of the phylogenetic tree using IQTree. The colors of the different branches signify different phyla. Matching recovered draft genomes in each sample (A–C for samples GA_E1-1, GA_E1-2, and GA_E2-1, respectively), i.e., genomes binned from these samples, are provided as green boxes (otherwise left white). The presence of marker genes based on a marker gene search using HMMs on these genomes for specific chemolithoautotrophic pathways is shown as green boxes (otherwise left white). C signifies carbon fixation with (1) CBB, (2) rTCA, and (3) WL, C₁ for C₁-metabolism with (4) carbon monoxide oxidation, (5) formaldehyde oxidation, and (6) methanol oxidation, O for oxygen metabolism with (7) cytochrome c bd, (8) cytochrome c bo, (9) cytochrome c caa₃, and (10) cytochrome cbb₃, H for hydrogen metabolism with (12) FeFe-Hydrogenases type A, (13) NiFe-Hydrogenases type 3b, (14) NiFe-Hydrogenases type 3c, (15) Nife-Hydrogenases, (16) NiFe-Hydrogenases type 4 and (17) NiFe-Hydrogenases type 1, N for nitrogen metabolism with (18) Nitrate reduction, (19) Nitric oxide reduction, (20) nitrite reduction and (21) nitrous oxide reduction, S for sulfur metabolism with (22) sulfide oxidation, (23) sulfite reduction with dsr, (24) sulfite reduction with asr, (25) sulfur oxidation with dsr, (26) sulfur oxidation with sor, (27) sulfur oxidation with sdo, (28) sulfate reduction via APS with sat and (29) Thiosulfate disproportionation. Olive bars show the average iRep value of the respective bacterial population, brown bars show the maximal growth rate of the representative genome as estimated by growth red, and blue bars show the average log10-scaled coverage. B Morphologies of microorganisms as determined via DAPI staining and fluorescence microscopy (scale bars = 5 µm) are shown. The morphologies were documented in two sampling campaigns (June 2016 and February 2018 with three and two samples in technical duplicates, respectively).

Fig. 2. In situ bacterial replication rates across subsurface ecosystems ordered by ecosystem depth.

The figure depicts a beeswarm plot of iRep values of genomes (x-axis) across ecosystems (y-axis) with genomes colored according to their predicted metabolic potential and the black dot representing the median iRep value (individual iRep values in Supplementary Data 2). C represents carbon, N₂ nitrogen, H₂ hydrogen, O₂ oxygen, and S sulfur. Colored squares depict the sample type. Samples impacted by geological degassing and a sediment sample along with the respective aquifer sample are plotted separately. The top y-axis shows the sampling depth of the different ecosystems (Supplementary Table 5). In total, 895 genomes were used for this analysis with ≧ 70% completeness and ≦ 10% contamination based on 51 bacterial and 38 archaeal single-copy genes. The order of samples is given in Supplementary Table 5. p-Values are derived from two-sided student’s t-tests. The exact p-values from top to bottom are p < 2.2 × 10¹⁶ (minimal value in R) and p = 0.0003934, respectively.

Fig. 3. Chemolithoautotrophic metabolic potential across ecosystems.

The heatmap shows the read-normalized abundance of chemolithoautotrophic pathways, Z-score scaled for the respective metabolisms. Colored squares on the right depict the sample type. If multiple biological replicates of samples were available, up to three were depicted. Sample order is according to Supplementary Table 5. Supplementary Fig. 5 and Supplementary Fig. 6 display the Z-scaled number of hits (Supplementary Fig. 5) or normalized abundance (Supplementary Fig. 6) of the individual genes aggregated into their pathways in this figure.

Fig. 4. Geographical distribution and chemolithoautotrophic potential of Ca. Altiarchaeota.

A Global map with locations from which Ca. Altiarchaeota genomes were recovered. B Metabolic potential of Ca. Altiarchaeota genomes. Genomes belonging to the Alti-1 clade are highlighted in dark gray, Alti-2 genomes in beige. If multiple genomes from a specific site were available, they were all used to identify the metabolic potential. The bar chart shows averaged growthpred-predicted minimal generation times across all genotypes recovered from a specific genome, with error bars denoting the averaged standard deviations (growthpred returns both an average minimal generation time and a standard deviation for this value). In addition, the mean minimal generation time for each genome is indicated by black dots. The circled numbers below the heatmap depict the genes identified as markers and stand for (1) codhC, (2) codhD, (3) rubisco form III, (4) fae, (5) fmtf, (6) mtmc, (7) NiFe-Hydrogenase group 4, (8) NiFe-Hydrogenase group 3b, (9) NiFe-Hydrogenase group 1, (10) FeFe-Hydrogenase, (11) hdh, (12) ars. C Phylogeny of Alti-1 genotypes based on 30 universal ribosomal proteins (5136 aa positions, IQTree JTTDCMut+F + G4) and using the Alti-2 genome IMC4 as the outgroup. Branch supports correspond to ultrafast bootstraps (1000 replicates), the SH-aLRT test (1000 replicates), and the approximate Bayes test, respectively (a tree with outgroup in Supplementary Fig. 8). Details on Altiarchaeales genomes in Supplementary Table 3.

Fig. 5. Metabolic capacities of Ca. Altiarchaeum pangenome.

Previously identified genes in Ca. Altiarchaeum hamiconexum IMS was used as the basis to query the other genomes of known Altiarchaea clade members (see Fig. 4 for all members used in this analysis). To expand the predictable metabolic capacity of the genomes, METABOLIC was used to annotate genes, which mainly resulted in peptidases and glycosylhydrolases. If multiple genomes copies per site were available, they were all used to query for the respective genes. All gene functions are listed in Supplementary Data 3.