First shotgun metagenomics study of Juan de Fuca deep-sea sediments reveals distinct microbial communities above, within, between, and below sulfate methane transition zones

Abstract

The marine deep subsurface is home to a vast microbial ecosystem, affecting biogeochemical cycles on a global scale. One of the better-studied deep biospheres is the Juan de Fuca (JdF) Ridge, where hydrothermal fluid introduces oxidants into the sediment from below, resulting in two sulfate methane transition zones (SMTZs). In this study, we present the first shotgun metagenomics study of unamplified DNA from sediment samples from different depths in this stratified environment. Bioinformatic analyses showed a shift from a heterotrophic, Chloroflexota-dominated community above the upper SMTZ to a chemolithoautotrophic Proteobacteria-dominated community below the secondary SMTZ. The reintroduction of sulfate likely enables respiration and boosts active cells that oxidize acetate, iron, and complex carbohydrates to degrade dead biomass in this low-abundance, low-diversity environment. In addition, analyses showed many proteins of unknown function as well as novel metagenome-assembled genomes (MAGs). The study provides new insights into microbial communities in this habitat, enabled by an improved DNA extraction protocol that allows a less biased view of taxonomic composition and metabolic activities, as well as uncovering novel taxa. Our approach presents the first successful attempt at unamplified shotgun sequencing samples from beyond 50 meters below the seafloor and opens new ways for capturing the true diversity and functional potential of deep-sea sediments.

Keywords: archaea; bacteria; deep biosphere; hydrothermal fluid; marine subsurface; metagenome assembled genomes; optimized DNA extraction.

Figure 1

(A) Scheme of the hydrogeological regime at the eastern flank of the JdF Ridge. Seawater enters the basement at the Grizzly Bare outcrop and circulates within the basement rock, where it undergoes chemical alteration. Sediment laying above the basement impedes fluid exchange with the ocean, but diffusion into sediment layers still occurs from above and below. At Baby Bare Seamount, the altered seawater from the basement is expelled again. The IODP drill site U 1301 and the depths from which the samples in this study were obtained are indicated (adapted after Engelen et al., 2008). (B) The chemical profile of pore water obtained during sampling by IODP in 2004 from cores U1301 C and D. The distribution of sulfate shows diffusion into the sediment both from above and below, while methane is only present when the sulfate concentration in the pore water is minimal, resulting in the formation of two sulfate methane transition zones (SMTZs) at ~30–50 m and at ~100–120 m (Fisher et al., 2005a) indicated as yellow bars. DIC: dissolved inorganic carbon. DOC: dissolved organic carbon. Data made available by IODP in 2008 (Ocean Drilling Data, 2021).

Figure 2

Alpha and beta diversities calculated from the abundance of contigs carrying 16S rRNA genes using ExpressBetaDiversity (Parks and Beiko, 2012). The 16S rRNA sequences used for diversity analysis were predicted from the contigs of the co-assembled metagenomes, and abundances were estimated by mapping the individual sample’s reads back onto the co-assembly. (A) Shannon diversity as a barplot. The diversity within the samples is shown while considering species abundance. The Shannon index weighs the abundance of highly abundant species higher. (B) ACE diversity as a barplot. The diversity within the samples is shown while considering species abundance. The ACE index places an additional focus on less abundant species. (C) Good’s coverage as a barplot. The estimated completeness of a sample is displayed as the negative decadic logarithm of the Good’s coverage, meaning that lower values indicate higher estimated completeness. (D) Jaccard diversity for every sample to every other sample as a heatmap. Displayed is the inverse of the similarity between every sample, meaning that two identical samples would have a dissimilarity of “0”.

Figure 3

(A) Information on the sample and key chemical components. The concentrations of methane, sulfate, and dissolved organic carbon (DOC) are shown as a simplified scheme; for a more detailed graph, see Figure 1. (B) Microbial abundance profiles of the four samples C1H2, C6H6, C11H1, and D2H3 based on rRNA, protein, and marker gene classification from IODP site U1301. Marker genes, 16S rRNA, and protein predictions were performed using the co-assembled metagenome, and their abundance within the original sample was determined using short-read mapping data. Classifications were obtained using the QIIME2 (Bolyen et al., 2019) classification outlined in the methods for 16S rRNA and by alignment against NCBI GenBank Release 235 for total protein predictions and marker gene predictions. Displayed are the classifications of identified sequences up to domain, phylum, and class, sequences that cannot be classified are displayed as part of the next higher level with no label. Sequences that cannot be attributed to a Domain are labeled as “unclassified.” CPR = Candidate Phyla Radiation. The Krona charts corresponding to the depicted figures can be found in Supplementary Data 2–4.

Figure 4

Heatmap of KO-Numbers collated into relevant functions associated with detected orthologues. OGs were predicted from the co-assembled metagenomes with the eggNOG Mapper (Huerta-Cepas et al., 2019) and their abundance was calculated by mapping the information of the corresponding encoding contigs back to the reads. Displayed is the relative abundance of the orthologues in each sample, normalized to the abundance of the universal housekeeping gene/function COG0468 (recA). In cases where closely related OGs were grouped into more general metabolic functions, the corresponding mean value is shown. Functional groups with a total relative abundance below 0.02 have been removed for better representability.

Figure 5

Metabolic profile of MAGs. The relative abundance of each MAG was determined based on short-read mapping information and is displayed as a heat map on the left side of the plot. The heat map shows the pathway completeness within the 42 selected metagenome-assembled genomes (MAGs) obtained from Juan de Fuca (JdF) samples that have more than 50% estimated completeness and less than 5% estimated contamination. To predict pathway completeness, OGs were predicted from the co-assembled metagenome using eggNOG Mapper (Huerta-Cepas et al., 2019) and checked to see whether they belonged to a MAG contig. OGs of MAGs were compared to the pathway modules outlined in Supplementary Table 3, and the number of orthologues found in a MAG is displayed as a percentage. Taxonomic classifications are based on GTDB-TK (Chaumeil et al., 2019) and the order in which MAGs are listed reflects the respective maximum likelihood-based phylogenetic relationships of the selected MAGs. Completeness and contamination estimates were obtained using MDMcleaner (Vollmers et al., 2022).