Metagenomic profiles of archaea and bacteria within thermal and geochemical gradients of the Guaymas Basin deep subsurface

Domain: earth_science
Published: 2023-11-27
Editorial outcome: major_revision

Abstract
Previous studies of microbial communities in subseafloor sediments reported that microbial abundance and diversity decrease with sediment depth and age, and microbes dominating at depth tend to be a subset of the local seafloor community. However, the existence of geographically widespread, subsurface-adapted specialists is also possible. Here, we use metagenomic and metatranscriptomic analyses of the hydrothermally heated, sediment layers of Guaymas Basin (Gulf of California, Mexico) to examine the distribution and activity patterns of bacteria and archaea along thermal, geochemical and cell count gradients. We find that the composition and distribution of metagenome-assembled genomes (MAGs), dominated by numerous lineages of Chloroflexota and Thermoproteota, correlate with biogeochemical parameters as long as temperatures remain moderate, but downcore increasing temperatures beyond ca. 45ºC override other factors. Consistently, MAG size and diversity decrease with increasing temperature, indicating a downcore winnowing of the subsurface biosphere. By contrast, specific archaeal MAGs within the Thermoproteota and Hadarchaeota increase in relative abundance and in recruitment of transcriptome reads towards deeper, hotter sediments, marking the transition towards a specialized deep, hot biosphere.

Introduction
The interplay between temperature stress and energy availability determines microbial survival in the subsurface biosphere, and delineates the extent and limits of life in the deep subsurface biosphere1,2. As microbial communities in cool, relatively shallow subsurface sediments transition into more deeply buried and increasingly warm and finally hot sediments, it should be possible to track how subsurface bacteria and archaea react to these gradually harsher regimes downcore on the levels of cellular activity and community change. While microbial abundance and diversity are generally expected to decline downcore3,4, it is also possible that particular subsurface-adapted microbial populations benefit from conditions that would eliminate others, and constitute a specialized deep, hot biosphere. Recent studies indicated active microbial populations in extremely deep and hot sediments, yet without sequence-based identification5,6. To learn more about bacterial and archaeal communities of the deep, hot biosphere from a genomic perspective, downcore trends of diversity and activity in increasingly hot sediments need to be examined, and microbial communities and their genomes have to be tracked downcore, as far as microbial biomass and DNA yield allow. Yet, investigating downcore changes in microbial abundance, community composition and activity in well-characterized geochemical and thermal gradients requires a suitable field site where extensive physical, chemical and microbial gradients can be sampled in adequate resolution by sediment coring and drilling.

An ideal natural laboratory for such a research task is Guaymas Basin, a hydrothermally-active ocean spreading center in the Gulf of California, covered by several hundred meters of sediment that host basaltic sill intrusions7 and strong geothermal heat flow8. Pyrolysis of buried organic carbon in these organic-rich sediments produces a complex milieu of petroleum hydrocarbons, including light hydrocarbons and methane, alkanes, and aromatic compounds, as well as carboxylic acids, and ammonia9,10. These compounds are transported via hydrothermal fluids through Guaymas Basin’s thick sediments, supporting diverse and active microbial communities11. Collectively, these communities not only perform chemosynthetic carbon fixation and heterotrophic organic matter remineralization, but they also assimilate fossil carbon into the benthic biosphere12. Yet, few studies to date have explored the microbiology of deep subsurface sediments in Guaymas Basin. Methanogens were enriched from sediments collected during Deep Sea Drilling Program Expedition 64 to Guaymas Basin13, and bacterial and archaeal communities in piston cores were surveyed using 16S rRNA amplicon sequencing14–16. Aside from these studies, the spatial extent, diversity and activity of the deep biosphere in Guaymas Basin have remained largely unknown.

International Ocean Discovery Program (IODP) Expedition 385 drilled into Guaymas Basin at eight locations that differ in their degree of hydrothermal influence and heat flow8, to survey their resulting characteristics17. Drilling sites followed broadly a northwest-to-southeast transect across the northern Guaymas axial trough (Fig. 1). Two neighboring sites (U1545 and U1546) on the northwestern end of Guaymas Basin18,19 essentially differ by the presence of a massive, thermally equilibrated sill between 350 to 430 meters below seafloor (mbsf) at Site U15467. Two drilling sites (U1547, U1548) targeted the hydrothermally active Ring vent area, approximately 28 km northwest of the spreading center15, where a shallow, recently emplaced and hot sill creates steep thermal gradients and drives hydrothermal circulation20. Drilling Site U154921 explores the periphery of an off-axis methane cold seep, Octopus Mound, located ~9.5 km northwest of the northern axial graben22.

These contrasting sites provide an opportunity for a comprehensive analysis of subsurface microbiota at different temperatures and depths. To assess the environmental distribution and genomic potential of microbes living in the deep biosphere of Guaymas Basin, we analyzed reconstructed metagenome-assembled genomes (MAGs) from depths ranging from 0.8 to 219.5 mbsf at these thermally and geochemically contrasting sites. We also provide evidence for the activity of specific bacterial and archaeal lineages by mRNA transcript mapping on bacterial and archaeal MAGs.

Results
Sampling sites and depths
Metagenomes were produced from sediment samples at drilling sites U1545B to U1549B that follow a northwest-to-southeast transect across the northwestern flanking region of Guaymas Basin (Fig. 1A) and include an off-axis hydrothermal system, the Ring vent site (Fig. 1B). The samples were selected to coordinate with depths used for separate ongoing analyses, and ranged from 1.7 m to 219.5 mbsf at Site U1545B, 0.8-16.3 mbsf at U1546B, 2.1-75.7 mbsf at U1547B, 9.1-69.4 mbsf at U1548B, and 16.5 mbsf at U1549B (Fig. 1; Table 1). For all samples, a wide range of geochemical parameters was analyzed shipboard (Supplementary Dataset 1). The sites represent distinctly different thermal gradients and cell densities; generally, sites with steeper downcore temperature gradients are characterized by more rapidly decreasing cell counts (Fig. 1C, D). U1545B is the reference site for IODP Expedition 385 because of the absence of seepage, hydrothermal influence, and massive sill intrusions18. Here, metagenome libraries extended down to 219.5 mbsf, at in-situ temperatures of 54.3°C. Cell count trends for sites U1545, U1546 and U1549 were similar, and showed a decrease over three orders of magnitude within 100 meters (Fig. 1C). At the hot Ring vent sites U1547B and U1548B8,20, comparable temperatures of 50–55°C were already reached near 70 mbsf (Fig. 1D), and cell counts decreased by four to five orders of magnitude within this depth range (Fig. 1C). To describe temperature-related trends in MAG recovery and diversity, we categorized our samples into three groups according to temperature; cool (2–20°C), warm (20–45°C) and hot (>45°C).

Subsurface Biogeochemical zonation
Most samples collected for metagenomes are from sediments within the sulfate-reducing zone where sulfate is still available at near-seawater concentrations (~28 mM) or becomes gradually depleted with depth (Table 1). At those same sediments hydrogen sulfide concentrations are gradually increasing towards multiple millimolar concentrations. Metagenome samples from site U1545B also include depths spanning the sulfate-methane transition zone (SMTZ) at ~64 mbsf where sulfate is consumed by microbial sulfate reduction, and methane begins to accumulate. At the SMTZ sulfate concentrations drop from 21.1 mM to 0.7 mM, sulfide reaches peak concentrations of 8.9 mM, and methane concentrations increase from picomolar to 1.5 mM (Table 1). High methane concentrations persist also in deeper samples from U1545B, and decrease only in the very deepest samples (>200 msbf). The deep subsurface methane reservoir at this and other sites results from long-term thermogenic and biological methane accumulation23. In contrast to site U1545B, samples from Ring vent sites U1547B and U1548B show gradual downcore sulfate consumption (from 27.9 to 18.8 mM) but not depletion, combined with hydrogen sulfide accumulation (max. 7.1 mM at 75.7 mbsf at U1547B); methane does not accumulate in these samples. Ammonia concentrations increase from <1 mM towards 3 to 5 mM downcore at most sites, and reach 9 to 25 mM below the SMTZ in U1545B. Dissolved inorganic carbon (DIC) and alkalinity concentrations are generally highest at Site U1545B where they peak in the SMTZ (~28 and 60 mM, respectively). Ammonia, DIC and alkalinity remain elevated not only in the upper sediment column but also in the deeper samples of Site U1545B, presumably due to cumulative bioremineralization of buried organic matter over time at this undisturbed site. In contrast, the Ring vent samples (sites U1547B and U1548B) generally have lower ammonia, alkalinity and DIC porewater concentrations, suggesting reduced remineralization of organic matter at these sites, most likely a consequence of hydrothermal activity due to recent volcanic sill emplacement15. Dissolved organic carbon (DOC) remained ~10 to 20 mg/L in most samples but increased towards 70 mg/L in the sulfate-methane transition zone of U1545B and remained between 20 and 50 mg/L in the deeper sediments of U1545B. This suggests DOC enrichment and decreased heterotrophic DOC consumption in deep methanogenic sediments of U1545B where energy-rich electron acceptors for heterotrophic carbon remineralization are not available. While total nitrogen and total organic carbon generally decrease with depth at all sites, the Ring vent sites have moderately elevated TOC/TN ratios (Table 1), likely reflecting the influence of nitrogen-depleted hydrothermal carbon sources14. Total petroleum hydrocarbon, saturated and polyaromatic hydrocarbon content remain each quite similar across a wide range of sediments and temperatures, before increasing considerably in hot sediments (>80°C) near deep sill intrusions (Supplementary Dataset 1, and Supplementary Figs. 1A, B).

MAG diversity, distribution, and evidence of activity
A total of 142 metagenome-assembled genomes (MAGs) were recovered from a co-assembly of all metagenomic samples (Supplementary Dataset 2). MAGs that matched those from negative controls were excluded from further analysis (Supplementary Dataset 3). For downstream analysis, we retained 89 bacterial and archaeal MAGs that had at least ≥50% bin completeness and ≤10% bin contamination24. Genome completeness ranged from ~50 to 97% (Supplementary Fig. 2; Supplementary Dataset 4).

Of these 89 MAGs, 26 MAGs were assigned to 6 archaeal phyla, and 63 MAGs were assigned to 13 bacterial phyla (Fig. 2); the phylogenetic spectrum includes lineages documented previously in 16S rRNA gene amplicon sequencing of shallow subsurface sediments14,15, and in metagenomic surveys of shallow hydrothermal sediments of Guaymas Basin25. In parallel to downcore decreasing cell numbers (Fig. 1), MAG diversity decreased downcore at all sites as temperatures increased (Fig. 2). In samples from cool (2–20°C) sediments from all sites, reads mapped to diverse bacterial and archaeal phyla, including the bacterial phyla Chloroflexota, Acidobacteriota, Desulfobacterota, WOR-3, Aerophobota, and Bipolaricaulota, and the archaeal phyla Thermoproteota, Thermoplasmatota, and Aenigmatarchaeota (Fig. 2). In samples with warm temperatures (20–45°C), reads were predominantly assigned to bacterial phyla Chloroflexota (mostly order-level group G1F9), Acidobacteriota, WOR-3 (order-level group UBA3073), Aerophobota and Bipolaricaulota, and to archaeal phyla Thermoproteota, Hadarchaeota, and Aenigmatarchaeota. At hot temperatures (45–60°C), bacterial reads mapped primarily to a single Chloroflexota MAG (class Dehalococcoidia), a single WOR-3 MAG and two Aerophobota MAGs (class Aerophobia). In contrast, several Archaeal MAGs show a marked preference for hot sediments, and mapped to the Thermoproteota (class Bathyarchaeia), and Hadarchaeota (class Hadarchaeia). Our recovered MAGs reflected metabolisms predicted for the deep biosphere including sulfur, nitrogen and methane cycling, hydrocarbon degradation, and carbon fixation26 (Supplementary Note, and Supplementary Figs. 3 and 4). Desulfobacterota MAGs linked to sulfate reduction contained the dsr operon (e.g., dsrB/J/K/D) that is essential for dissimilatory sulfate reduction27, and were recovered from shallow sediments with available sulfate. These MAGs also shared potential for iron reduction using extracellular electron transfer mechanisms, such as mtrA, mtoA, and eetB genes28,29. Marker genes of dissimilatory iron reduction28 (e.g., dmkA, dmkB, eetA, eetB, fmnA, fmnB, pplA, ndh2) with the potential for extracellular electron transfer (EET) were identified in 86/89 MAGs within all recovered phyla (Supplementary Note, and Supplementary Datasets 5, 6).

To determine any intra-phylum differences in metabolic activity, we mapped reads of the Guaymas Basin subsurface metatranscriptome30 to our recovered MAGs, for samples collected at the same sites (Fig. 3). Since the metagenome and metatranscriptome of the Guaymas Basin subsurface remain incompletely covered by sequence data, the absence of transcript read mapping to particular MAGs cannot be taken as evidence of microbial inactivity. Microbial activity of the deep biosphere is certainly constrained but not eliminated by substrate and energy limitation1. To avoid these ambiguities that are inherent in negative transcript mapping results, we focus on positive transcript mapping results that support the activity of specific MAGs in the subsurface. Actively transcribed genes are present for MAGs within all phyla discussed here, albeit at variable levels; some MAGs within individual phyla show no or much lower apparent activity than others (Fig. 3).

Most transcriptionally active bacterial and archaeal MAGs from warm and hot sediments belong to uncultured lineages previously detected in hydrothermal chimneys, sulfidic springs and seeps, and in Guaymas Basin surficial hydrothermal sediments (Fig. 3). Bacterial transcripts from warm sediments were affiliated with four MAGs (GMP_018, GMP_083, GMP_036, and GMP_057) of the Chloroflexota GIF19 lineage, a dominant group in carbonate hydrothermal chimneys31. Other transcripts from warm sediments mapped to MAG GMP_019 within the dehalogenating Dehalococcoides lineage, to MAGs GMP_007 and GMP_011 within the subsurface Aerophobota, and MAG GMP_58 within the Bipolaricaudota lineage UBA7950, found at the Lost City hydrothermal vents32. Archaeal transcripts in warm and hot sediment samples were mapped to MAG GMP_008 within the Hadarchaeota, MAG GMP_075 of the Aenigmatarchaeota QMZP01 lineage from a terrestrial sulfur spring33, MAG GMP_040l within the thermoproteotal brine pool lineage TCS6434, and to three Thermoproteota MAGs GMP_002, GMP_026, and GMP_039 within the B26-1 lineage from Guaymas Basin hydrothermal sediments35. Transcriptional activity of these MAGs suggests their inherent physiological adaptations to warm and reducing habitats are advantageous in the Guaymas Basin subsurface as well.

The influence of environmental factors on MAG composition
The relationship between environmental parameters (Supplementary Dataset 1) and the taxonomic composition of MAGs from cool (2–20°C), warm (20–45°C), and hot temperatures (45–62°C) was investigated using non-metric multidimensional scaling (nMDS) (Fig. 4) and Canonical Correspondence Analysis (CCA) (Supplementary Fig. 5). In both analyses, bacterial and archaeal MAGs clustered consistently by temperature, comparable to previous analyses of temperature-dependent microbial community composition in surficial sediments in Guaymas Basin22. In particular, MAGs from hot sediments aligned with temperature as the strongest influencing factor (Fig. 4). Total sulfide concentration (H2S) was aligned with temperature in the CCA plot (Supplementary Fig. 5). For samples from warm sediments, nMDS analyses revealed that methane, alkalinity, and dissolved inorganic carbon (DIC) concentrations exerted a significant effect (p < 0.05) on the MAG community (Fig. 4). For cool sediments, nMDS and CCA showed consistently that MAGs clustered in the direction of total organic carbon (TOC) and total nitrogen (TN) content, CO concentration, pH and salinity (Fig. 4, Supplementary Fig. 5). The influence of TN and TOC on MAG diversity in cool samples may reflect increased availability of labile sources of dissolved and particulate organic matter in near-surface sediments. The consistent impact of pH and salinity on MAG diversity in cool samples, in both CCA and nMDS analyses, reflects persistent downcore trends towards lower pH and slightly reduced porewater salinity (Supplementary Dataset 1).

To summarize, the environmental parameters that impact MAG composition change downcore, from surface-linked factors such as TN, TOC, pH and salinity that impact MAGs in cool sediments, to biogeochemical parameters reflecting terminal organic matter degradation, such as increasing DIC, alkalinity and methane concentrations, in deeper and warmer sediments. An organic substrate-depleted, DIC- and methane-enriched deep subsurface environment may select for specific phyla or taxa with autotrophic capabilities (e.g., Hadarchaeota). For MAGs from deep and hot samples, carbon or nitrogen substrates, or other chemical factors, become secondary to the impact of temperature itself.

Metagenomic features with wide subsurface distribution
In addition to genes for core metabolic processes (e.g., glycolysis, biosynthesis of nucleotides and amino acids), Guaymas Basin MAGs contain widespread genomic features that extend across multiple bacterial and archaeal phyla. Some of these widespread genomic features have obvious adaptive value and are thus retained for survival, while others challenge assumptions on subsurface adaptations and evolutionary constraints under subsurface conditions.

Among genes that confer survival advantages, two-component systems (TCSs) can induce metabolic shifts, and are used extensively by bacteria and some archaea to respond and adapt to environmental changes36. Generally, archaea acquire TCS genes through horizontal gene transfer from bacteria37. In Guaymas Basin, TCS genes occur in the majority of bacterial MAGs but not in archaeal MAGs (Supplementary Datasets 5, 6), and they may help cells to adapt to long term burial. For example, the KinABCDE-Spo0FA system is present in almost all our bacterial MAGs and plays a role in sporulation by shifting cellular metabolism from active growth to dormancy/sporulation38. Likewise, the RegB/RegA redox-signaling mechanism involved in carbon fixation, hydrogen oxidation and anaerobic respiration39 is present in the majority of the bacterial MAGs.

Among widely distributed genes, we found a variety of transporters and efflux pumps associated with microbial defense, and biosynthetic gene clusters involved in synthesis of diverse secondary metabolites (Supplementary Datasets 5, 6). While all bacterial and archaeal MAGs encoded transporters26, efflux pumps (found in 67% of all MAGs) included a large proportion of multidrug resistance pumps, detected in 44% of all MAGs. In 25% of all MAGs, biosynthetic gene clusters were involved in the biosynthesis of diverse secondary metabolites (Supplementary Dataset 7). Archaeal biosynthetic gene clusters were primarily annotated as polyketide synthases, ribosomally synthesized and post-translationally modified peptides. Additionally, archaeal genes encoded the synthesis of terpenes (e.g., geranylgeranyl diphosphate synthase; Supplementary Dataset 6) that can be part of their lipid membranes, or function as pigments, antimicrobial agents and (in plants) as thermoprotectants40.

Genes involved in chemotaxis (cheA/B/R/W/Y) and motility (flgB/C/E/G/H/I and fliE/F/G/) were present in 56% of all bacterial and archaeal MAGs (Supplementary Datasets 5, 6). These findings suggest potential for cell-cell interaction, cell movement and competition for resources in the Guaymas Basin subsurface microbial community – a surprising result given deep biosphere microorganisms are trapped in tight pore spaces that limit movement and interaction41. While cell motility genes are gradually depleted downcore in the marine subsurface42, they do not disappear. Cell motility and secondary metabolite biosynthesis genes were present and expressed in marine subsurface MAGs from Peru Margin and Canterbury Basin sediments at depths down to 345 mbsf43.

Our results can be interpreted as evidence that evolution in the deep biosphere proceeds at extremely slow rates. Cells deep below the sediment surface must use available energy to maintain their cellular integrity over possibly geological timeframes while greatly attenuating cell division and genome replication1,41, unless some physical disturbance or fluid flow returns them to the sediment surface. Under subsurface conditions, attenuated gene loss slows down the impact of selection that gradually shapes the subsurface biosphere3. However, the adaptive value for genes of motility and competition may be reduced with depth but is unlikely to expire entirely, since pore space constraints do not preclude slow microbial movement (millimeters over months), as demonstrated experimentally by gradual recolonization of deep subsurface sediments44.

Characteristics and distribution of dominant bacterial and archaeal groups
The Guaymas Basin subsurface yields predominantly MAGs affiliated with specific phylum- and order-level lineages that show distinct mesophilic and thermophilic preferences. This suggests lineages appearing at specific depths and temperature ranges respond to environmental factors, which in turn shape their occurrence patterns. Our central working hypothesis is that the Guaymas Basin subsurface community is not a random assemblage, but reveals phylogenetic and functional structure that can be tracked downcore. Our account of this structured community focuses on dominant bacterial and archaeal phyla (Chloroflexota, Thermoproteota, Hadarchaeota); an extended overview on further bacterial and archaeal MAGs is provided in the Supplementary Note.

Dominant subsurface bacteria
Of 63 bacterial MAGs found in the Guaymas Basin subsurface, 23 are members of the phylum Chloroflexota, one of the dominant phyla in marine sediments with metabolically diverse fermentative and dehalogenating lineages45 (Supplementary Note). Within the Guaymas subsurface, Chloroflexota MAGs comprise 12 order-level lineages, and account for a significant fraction of recruited metagenomic reads per sample (up to 8.3%) (Fig. 2, Supplementary Fig. 6). At site U1545B, Chloroflexota MAGs were widespread within cool samples (2–20°C) and persisted occasionally into deep and warm sediments; at Ring vent site U1547B they were ubiquitous in cool samples but also widely found in warm sediments (20–45°C) (Fig. 2, Supplementary Fig. 6). MAGs that occur in warm sediments are affiliated with the subsurface and hydrothermal GIF9 group31,46, the VGOG01 lineage from the sulfidic, warm water column of tropical Lake Tanganyika47, and the dehalogenating Dehalococcoidales lineage. In hot sediments above 45°C, Chloroflexota MAGs appear only in traces (Supplementary Fig. 6). Thus, the Guaymas Basin subsurface Chloroflexota generally prefer cool or moderately warm habitats, and avoid temperatures above ca. 40°C.

Metagenomes were assembled and annotated for all Chloroflexota in our datasets to gather additional information about their metabolic potential (Supplementary Fig. 3, and Supplementary Dataset 8). Using the KEGG framework for functional annotation, within the general category “central carbohydrate metabolism” we find core genes that can participate in the TCA cycle, glycolysis, gluconeogenesis, and the pentose phosphate pathway (Supplementary Fig. 3). This category includes one specific module (K0378) that encodes an aldehyde ferrodoxin oxidoreductase (AOR), a tungsten-containing enzyme identified in mesophilic bacteria that can reduce aromatic compounds48. Within the category of “other carbohydrate metabolism” we find genes affiliated with galactonate/galactose degradation that can be linked to biosynthesis of alkaloids (e.g., terpenoid alkaloids). We also detect modules (assigned as “photorespiration”) that are involved in the glycine cleavage system and shared between different amino acid biosynthetic pathways (Supplementary Note). In the general category “metabolic capacity” we detected genes assigned to oxygenic and anoxygenic photosynthesis, nonetheless, these are genes (e.g., pyruvate phosphate dikinase and citryl-CoA lyase) involved in carbon fixation. An expanded KEGG module analysis of the whole community metagenome (Supplementary Dataset 9) reveals many of the same genes, including those within the categories “central carbohydrate metabolism” and “other carbohydrate metabolism”, reflecting the dominance of Chloroflexota among the recovered MAGs Hadarchaeota (Supplementary Fig. 4).

Dominant subsurface archaea
Although the archaea contributed only 26 MAGs compared to 63 bacterial MAGs to our total, and represent fewer phylum-level lineages, they exhibit greater thermal range (Fig. 2). MAGS of two dominant archaeal phyla—the Thermoproteota (11 MAGs) and the Hadarchaeota (4 MAGs)—prefer warm and hot subsurface sediments, and are introduced here in greater detail; additional archaeal lineages are discussed in the Supplementary Note.

Archaeal MAGs were dominated by the Thermoproteota, an archaeal phylum consisting of four major lineages, the Thaumarchaeota, Aigarchaeota, Korarchaeota and Bathyarchaeia49. All 11 Thermoproteota MAGs belonged to the uncultured class Bathyarchaeia; these were detected at all examined sites but primarily at the Ring vent sites U1547B and U1548B (Fig. 2 and Supplementary Fig. 6). Order-level identification of bathyarchaeial MAGs reveals linkages to subsurface, seep and hydrothermal sediment habitats. Five MAGs assigned to the order-level lineages TCS64 and 40CM-2-53-6 were recovered primarily between 0.8-15 mbsf sediments at sites U1545B and U1547B with cool temperatures ranging from 2.8–17.4°C. These bathyarchaeial orders have been reported also from deep sea brine pool samples34 and from soil samples50. The order-level lineage B26-1, previously found in Guaymas Basin hydrothermal sediments25,35, included three MAGs from warmer sediments (19–40°C) below 63.8 mbsf at U1545B, and from warm to hot sediments (24–47°C) between 19.3 and 65.8 mbsf at Ring vent site U1547B. The order-level lineage RBG-16-48-13, recovered previously from terrestrial subsurface cores51, was represented by a MAG detected at site U1548 at 20–45°C (Fig. 2). Two bathyarchaeial MAGs could not be classified at the order level, but one of these MAGs was abundant at temperatures between 39.5–47°C at U1547B (Supplementary Fig. 6). The detection of bathyarchaeial MAGs over a wide temperature spectrum, and the link of bathyarchaeial orders to specific temperature regimes, suggests distinct thermal preferences among different lineages of Bathyarchaeia52. The ubiquitous presence of Bathyarchaeia in anaerobic sediments53, including hydrothermal sediments35, can be attributed to their capacity to metabolize multiple organic substrates, e.g., polysaccharides, urea, acetate, detrital proteins, and aromatics compounds such as benzoate and lignin54, potential substrates in the hydrocarbon-rich Guaymas Basin subsurface. Based on MAG gene content, Bathyarchaiea can potentially utilize formaldehyde and shuttle it into carbon fixation via the Wood-Ljungdahl pathway (Supplementary Note). Lineage-specific thermophilic adaptations among the Bathyarchaeia include reverse DNA gyrase that facilitates DNA supercoiling under extreme temperatures54.

Hadarchaeota thrive in subsurface sediments by a combination of heterotrophic traits (fermentation of carbohydrates) with autotrophic energy generation, specifically the oxidation of carbon monoxide and hydrogen55. Hadarchaeota were previously recovered from surficial hydrothermal sediments in Guaymas Basin25. Consistently, the 4 hadarchaeotal MAGs (GMP_008, GMP_020, GMP_027, GMP_034) did not recruit any reads from cool samples but only from warm and hot samples, indicating a preference for elevated temperatures (Fig. 2). In contrast to changing thermal preferences for MAGs from different bathyarchaeial orders, the Hadarchaeota, originally detected in hot and deep terrestrial aquifers56, consistently prefer elevated temperatures in deep sediments of the Guaymas Basin subsurface (Fig. 2).

Hadarchaeotal genomic features
Abundant and highly expressed hadarchaeotal MAGs were examined for characteristic features in their genomes and transcriptomes. One Hadarchaeota MAG, Guaymas_P_008, recruited ~19% of all metagenomic reads at 74.3 msbf (in-situ temperature 51°C) at Ring vent site U1547B (Figs. 2, 3). This MAG contained genes for carbohydrate hydrolysis (α-RHA, β-galactosidase) and nucleoside uptake and degradation (nucleoside transporters, purine nucleosidases) that suggest purine/pyrimidine synthesis from nucleosides. This MAG also contained carbon monoxide oxidation genes (coxM, coxS) that were absent in the other Hadarchaeota MAGs that encoded genes for fermentation (porA, ack, acdA) and aromatics degradation (ubiX) (Supplementary Datasets 6, 7). The ability to utilize a wider range of carbohydrates may support higher temperature tolerance, as reported for thermally-adapted Bathyarchaeia genomes52. The potential for hydrocarbon utilization in Hadarachaeota and other phyla (Supplementary Note) might contribute to reduced hydrocarbon concentrations at intermediate sediment depths and temperatures (Supplementary Fig. 1). One of our Hadaracheota MAGs (P_034) contained homologs to mcrC and mcrG that regulate the expression and assembly of the alkyl/methyl coenzyme M reductase operon57, the essential methane and alkane-activating genes in archaeal methanogens, methane oxidizers and short-chain alkane oxidizers58. Finally, we note the presence of KaiC histidine in Hadarchaeota, a circadian clock protein that regulates cell division and allow prokaryotes to adapt to changes in environmental conditions59, and the gene for programmed cell death (protein 5) that is linked to anti-virus defense and triggers dormancy under hostile conditions60.

Genome size trends in the deep biosphere
Comparisons of estimated genome sizes for all MAGs that recruited at least 0.1% of metagenomic reads from cool, warm, and hot sediments revealed a difference in average genome size. The most abundant genomes in cool sediments were on average significantly (two-sided partially overlapping samples t-test, adjusted p < 0.05) larger (~32%) than those recovered from hot sediments (Fig. 5A, B). The estimated genome size of MAGs recovered from our shallow (2–15 mbsf) samples was also ~22% larger on average than those detected in deeper (>60 mbsf) and warmer (>30–40°C) sediments. Linear regression analysis demonstrated a general reduction in average genome size in our samples as both temperature and depth increased (Fig. 5C, D). Elevated in situ temperatures are thought to select for smaller genome sizes via genome streamlining61, for example increased gene loss after duplication; the effects of genomic streamlining are pervasive and result in the elimination of hundreds of genes all over the genome62. Reduced genome size lowers the metabolic cost required for microbial DNA replication, as suggested for hadal microorganisms in the Challenger Deep at Mariana Trench63. Microbes with smaller genomes would gain a relative survival advantage and gradually dominate the microbial community in the subsurface, as metabolically more demanding microbial community members with large genomes die off. Such a mechanism would contribute to the selection of subsurface-adapted microbial communities that has been documented already within the top few meters below seafloor3, and it would explain the small size of microbial cells in deep subsurface sediments, near 0.5 micrometer64.

Temperature impact on MAG recovery
The environmental stresses that increasingly exclude microbial lineages, reduce genome size and reduce overall microbial population size (and thus, quantity of recovered DNA) are reflected in decreased recovery of MAGs in warmer and deeper samples from all sites (Fig. 6A, B). Plotted against depth, MAG recovery declines more quickly for the two hotter Ring vent sites U1547B and U1548B than for the cooler sites U1545B, U1546B and U1549B (Fig. 6B). When plotted against temperature, declining MAG recovery for the hot Ring vent sites and the cooler sites converge towards a shared minimum between ca. 50 and 60°C (Fig. 6A). These comparisons show that the decline of MAG recovery with depth is locally modified, behaves differently at different sites, and does not follow a uniform depth-related decay rate. In contrast, the influence of increasing temperature is pervasive, reduces microbial diversity at all sites, and occludes the emergence of MAGs representing new microbial lineages beyond approximately 50–60°C.

Conclusions and outlook
While improved DNA and RNA recovery could potentially compensate for declining downcore cell density, and extend the recovery of new bacterial and archaeal MAGs towards deeper and hotter sediments, the observed trend towards increasingly limited microbial diversity in the subsurface stands in marked contrast to the numerous bacterial and archaeal lineages that thrive in surficial hydrothermal sediments of Guaymas Basin, where fluidized sediments are permeated by pulsating, extremely hot (>80°C) and highly reducing fluids25. We ascribe the difference to contrasting energy supply65, and suggest that relatively moderate temperatures in IODP boreholes have a disproportionally greater impact on the energy-limited microbial deep biosphere, whereas surficial microbial communities that are well-supplied with energy-rich circulating hydrothermal fluids can tolerate high temperatures. The latter conditions select for thermophilic and hyperthermophilic, frequently chemolithoautotrophic bacteria and archaea66,67. We suggest that this difference ultimately results in distinct microbial communities in surficial hydrothermal sites, and in subsurface sediments where dominant bacteria and archaea (Chloroflexota, Thermoproteota, Acidobacteriota, Desulfobacterota) resemble the largely heterotrophic and mesophilic microbiota of non-hydrothermal benthic sediments68,69.

Yet, we note that specific archaea, in particular the Hadarchaeota, show a preference for deep, hot sediments of Guaymas Basin. These archaea extend consistently into the deep sediment column, not only by MAG detection but also in 16S rRNA gene surveys30, and appear to represent deep subsurface thermophiles that are sustained by substrates and energy sources of deep, hot sediments. Observations of microbial cells and activity in extremely deep, hot subsurface environments5,70, could indicate such thermophile communities that have adapted to deep subsurface conditions. Since candidate archaea for deep, hot biosphere communities were consistently detected in the hydrothermally influenced Ring vent sites where a hot volcanic sill is emplaced into organic-rich marine sediments, we extrapolate further that the mineralogically and morphologically complex basalt interface20 could provide microbial substrates and energy sources71, calling for further studies.

Methods
Sample collection
Sediment cores were collected during IODP Expedition 385 using the drilling vessel JOIDES Resolution. Holes at each site were first advanced using advanced piston coring (APC), then half-length APC, and then extended core barrel (XCB) coring as necessary. Temperature measurements used the advanced piston corer temperature (APCT-3) and Sediment Temperature 2 (SET2) tools8. Downhole logging conducted after coring used the triple combination and Formation MicroScanner sonic logging tool strings. After bringing core sections onto the core receiving platform of the D/V JOIDES Resolution, whole round samples for microbiology were retrieved within ~30 minutes using ethanol-cleaned spatulas. Samples for biogeochemical measurements were obtained and processed shipboard17. Whole round samples for DNA-based studies were capped with ethanol-sterilized endcaps, transferred to the microbiology laboratory, and stored briefly at 4°C in heat-sealed tri-foil gas-tight laminated bags flushed with nitrogen until processing. Masks, gloves and laboratory coats were worn during sample handling in the laboratory where core samples were transferred from their gas-tight bags onto sterilized foil on the bench surface inside a Table KOACH T500-F system, which creates an ISO Class I clean air environment (Koken Ltd., Japan). In addition, the bench surface was targeted with a fanless ionizer (Winstat BF2MA, Shishido Electrostatic Co., Ltd., Japan). Within this clean space, the exterior 2 cm of the extruded core section were removed using a sterilized ceramic knife. The core interior was transferred to sterile 50-mL Falcon tubes, labeled, and immediately frozen at –80°C for post cruise analyses. For RNA-based studies, sampling occurred immediately after core retrieval on the core receiving platform by sub-coring with a sterile, cut-off 50 cc syringe into the center of each freshly cut core section. These sub-cores were immediately frozen in liquid nitrogen and stored at –80°C.

DNA extraction and sequencing
DNA was extracted from selected core samples using a FastDNA SPIN Kit for Soil (MP Biomedicals). Up to 5 grams of sediment were processed following a modified manufacturer’s protocol72. Briefly, each sediment sample was homogenized twice (vs. once that the manufacturer suggests) in Lysing Matrix E tubes for 40 seconds at speed 5.5 m/s, using the MP biomedicals benchtop homogenizer equipped with 2 ml tube adaptors. Between the two homogenization rounds the samples were placed on ice for 2 minutes. After the second homogenization the samples were centrifuged at 14,000 x g for 5 minutes. For each sample, the supernatant and the top layer of the pellet was transferred to a clean 2 ml tube where proteins were precipitated by the addition of the protein precipitation solution (PPS) provided in the extraction kit. The rest of the extraction protocol followed the manufacturer’s recommendations. When parallel extractions were performed, the extracts were pooled and concentrated using EMD 3kDa Amicon Ultra-0.5 ml Centrifugal Filters (Millipore Sigma). A control extraction, in which no sediment was added, was included to account for any laboratory contaminants (Supplementary Materials). All libraries for metagenome sequencing (n = 29; 26 samples and 3 controls; Supplementary Data 2) were prepared from genomic DNA extracts that were submitted at the University of Delaware DNA Sequencing & Genotyping Center. Thirteen libraries were sequenced with NovaSeq S4 PE150 (Illumina) at the University of California, Davis Genome Center, and thirteen libraries were sequenced with NextSeq 550 (Illumina) at the University of Delaware DNA Sequencing & Genotyping Center. Metagenome sequence reads were deposited to the National Center for Biotechnology Information Sequence Read Archive under access numbers SRR23614663-23614677 and SRR22580794-SRR22580807 (Bioproject PRJNA909197).

Metagenomic co-assembly, binning, dereplication and taxonomic assignment
Metagenomic reads originating from adjacent regions (such as adjacent depths targeted in this study) are likely to overlap in the sequence space, increasing the mean coverage and extent of reconstruction of MAGs when using a co-assembly approach. Before assembly, reads were trimmed for quality and adapters removed using Trimmomatic v0.3973 (parameters: leading: 20; trailing: 20; sliding window: 0-24; min length 50). The quality of reads was verified with FastQC v0.11.974. For MAG reconstructions, we used the trimmed reads of metagenomic datasets from all 29 Guaymas samples sequenced in this study (Supplementary Data 2). The 26 metagenomes were co-assembled into contigs using MEGAHIT 1.2.975 with default parameters. For determining non-redundant MAGs, assembled contigs were binned using three different binners, MetaBAT2 2.12.18376, MaxBin 2.2.777, as well as CONCOCT 1.1.078. Output bins from all three binning algorithms were refined and dereplicated using DASTool 1.1.679. DasTool determines a unique MAG through a single-copy gene (SCG) scoring strategy coupled to an iterative bin de-replication procedure that produces the highest-scoring set of non-redundant bins (in terms of SGC completeness/contamination) from input bins generated by different binners. Completeness, size, and contamination levels of the reconstructed genomes were estimated using CheckM2 1.0.026. Only MAGs that were at least 50% complete and contained less than 10% contamination were used for downstream analyses (Supplementary Data 4). The taxonomic placement of the MAGs was performed with GTDB-Tk 2.1.080.

To account for seawater and laboratory contamination (Supplementary Note), control samples (Supplementary Data 2) identified MAGs of lab/control contaminants, including Patescibacteria (Paceibacteria, Microgenomatia), Actinobacteriota (Actinomycetia, Humimicrobia), Gammaproteobacteria (Pseudomonadales, Burkholderiales), and Firmicutes (Staphylococcales); these were removed from downstream analyses (Supplementary Data 3).

Calculation of MAG relative abundances
Metagenomic reads from 26 samples were mapped to each MAG using the CoverM 0.6.1 (https://github.com/wwood/CoverM) command line tool with the BWA 2.0 aligner81. The CoverM tool automatically concatenated all the MAGs into a single file, and metagenomic reads were recruited to MAG contigs, setting the parameter --min-read-percent-identity to 95 and --min-read-aligned-percent to 50. The “Relative Abundance” CoverM method on the “genome” setting was used to calculate the percent of total metagenomic reads per sample that mapped to each of the 89 MAGs. A custom R script was utilized to concatenate all coverM output files into a single file in a matrix format (with each sample representing a column and each row representing total percent of DNA-Seq reads per sample that mapped to a MAG) and was used for heatmap plotting.

Gene annotation, and prediction of KEGG metabolic module presence/absence using MetaPathPredict
Genes were called for all MAGs using Prodigal 2.6.382 and then annotated using Prokka 1.14.683, KofamScan 1.3.084, and METABOLIC 4.085 using default settings. KEGG modules for bacterial MAGs were reconstructed using gene annotations from the KofamScan 1.3.0 command line tool, and the presence or absence of incomplete modules in the genomes was predicted using MetaPathPredict 1.0.086 with default settings. MetaPathPredict cannot yet be applied to archaeal MAGs. Briefly, Prodigal was used to call genes, and KofamScan was used to annotate them. Gene annotations were generated for predicted genes from bacterial MAGs, and were used as input to MetaPathPredict, which generated predictions for the presence or absence of KEGG modules based on the gene annotations of all bacterial MAGs.

CCA and nMDS analyses of metagenomic abundance datasets and associated environmental parameters
The abundances of metagenomic reads mapped to MAGs were normalized using the “transcripts per million” normalization method87 with the read mapping "counts" output from coverM (https://github.com/wwood/CoverM). The abundance data were analyzed using canonical correlation analysis (CCA) as well as non-metric multidimensional scaling (nMDS) and were fitted with the environmental parameters in Supplementary Data 1 using R88. The cca and metaMDS functions were used for CCA and nMDS analyses, respectively, as well as the envfit function from the vegan 2.6-4 package89. The results were plotted using ggplot2 3.3.690 with sample shapes corresponding to temperature regime.

Estimated genome size analysis
The estimated genome size of all 89 MAGs was calculated by dividing the MAG assembly size (total base pair length of the MAG) by the fractional CheckM2 completeness of the MAG (the default CheckM2 completeness output divided by 100; a number between 0 and 1). Difference in genome size distributions for MAGs that recruited at least 0.1% of metagenomic reads from samples across temperature (cool [2–20°C], warm [20–45°C], hot [45–62°C]) and depth (shallow [2–15 mbsf], intermediate [15–60 mbsf], deep [>60 mbsf]) regimes was assessed using the two-sided partially overlapping samples t-test91, and resulting p-values were adjusted for multiple comparisons via Benjamini-Hochberg correction. The average estimated genome size of MAGs that recruited at least 0.1% of reads from metagenomic samples (n = 26) was fitted using linear regression against temperature and depth measurements affiliated with the samples.

MAG recovery at sample sites versus temperature and depth
The number of non-redundant MAGs that recruited at least 0.1% of reads from metagenomic samples (n = 26) was fitted using linear regression against temperature and depth measurements affiliated with the samples. Temperature values were interpolated for each sample using linear regression of the local thermal gradient (°C/m) multiplied by depth (mbsf), plus the y-axis intercept: U1545B, T = 0.225 x depth + 4.899; U1546B, T = 0.221 x depth + 2.627; U1547B, T = 0.511 x depth + 13.01; U1548B, T = 0.804 x depth + 6.499; U1549A/B, T = 0.194 x depth + 3.532. Reads from the 13 metatranscriptome samples collected from sites that metagenomic samples were taken from were mapped to each MAG using the CoverM 0.6.1 (https://github.com/wwood/CoverM) command line tool with the BWA 2.0 aligner81. The CoverM tool automatically concatenated all the MAGs into a single file, and metatranscriptome reads were recruited to MAG contigs, setting the parameter --min-read-percent-identity to 95 and --min-read-aligned-percent to 50. A custom R script was utilized to concatenate all coverM output files into a single file in a matrix format, with each sample representing a column and each row representing total percent of RNA-Seq reads per sample that mapped to a MAG. The output was used in this study for heatmap plotting to examine evidence for activity of the taxa for which were covered MAGs. Metatranscriptome reads were deposited to the National Center for Biotechnology Information Sequence Read Archive under accession numbers SRR22580929-SRR22580947 (Bioproject PRJNA909197).

Scanning of MAGs for secondary metabolite biosynthetic gene clusters
All 89 MAGs were individually scanned for secondary metabolic biosynthetic gene clusters using antiSMASH 6.092 with default parameters. Resulting gene cluster prediction results (in GenBank format) were parsed and their gene content was analyzed. Clusters with a total length less than 5kb were discarded from downstream analysis to minimize the inclusion of fragmented biosynthetic clusters in the analysis.

RNA extraction, library preparation, sequencing, and mapping of RNA reads to the MAGs
RNA was extracted from 19 sediment samples from sites U1545B-U1552B and a blank sample (control) using the RNeasy PowerSoil Total RNA Kit (Qiagen) following the manufacturer’s protocol with modifications which are discussed below. RNA samples were prepared from samples spanning the depths 0.8 to 101.9 mbsf. All samples, including a blank control, were first washed twice with absolute ethanol (200 proof; purity ≥ 99.5%), and sterile DEPC water (once) to reduce hydrocarbons and other inhibitory elements that otherwise resulted in low RNA yield. In brief, 13-15 grams of frozen sediments were transferred into UV-sterilized 50 ml Falcon tubes (RNAase/DNase free) using clean, autoclaved and ethanol-washed metallic spatulas. Each sample transferred into the 50 ml Falcon tube received an equal volume of absolute ethanol and was shaken manually for 2 min followed by 30 seconds of vortexing at full speed to create a slurry. Samples were spun in an Eppendorf centrifuge (5810R) for 2 minutes at 2000 x g. The supernatant was decanted and after the second wash with absolute ethanol, an equal volume of DEPC water was added into each sample and samples were spun for 2 minutes at 2000 x g. The supernatant was decanted, and each sediment sample was immediately divided into three 15 mL Falcon tubes containing beads provided in the PowerSoil Total RNA Isolation Kit (Qiagen). The RNA extraction protocol was followed as suggested by the manufacturer, with the modification that the RNA extracted from the three aliquots was pooled into one RNA collection column. All steps were performed in a UV-sterilized clean hood equipped with HEPA filters. Surfaces inside the hood and pipettes were thoroughly cleaned with RNase AWAYTM (Thermo Scientific™) before every RNA extraction and in between extraction steps.

Trace DNA contaminants were removed from RNA extracts using TURBO DNase (Thermo Fisher Scientific) and the manufacturer’s protocol. Removal of DNA was confirmed by negative PCR reactions using the bacterial primers BACT1369F/PROK1541R (F: 5’CGGTGAA-TACGTTCYCGG3’, R: 5’AAGGAGGTGATCCRGCCGCA3’) targeting the 16S rRNA gene93. Each 25 μl PCR reaction was prepared using 0.5 U μl–1 GoTaq® G2 Flexi DNA Polymerase (Promega), 1X Colorless GoTaq® Flexi Buffer, 2.5 mM MgCl2, (Promega) 0.4 mM dNTP Mix (Promega), 4 μM of each primer (final concentrations), and DEPC water. PCR reactions used an Eppendorf Mastercycler Pro S apoprotect (Model 6321) thermocycler with following conditions: 94°C for 5 min, followed by 35 cycles of 94°C (30 s), 55°C (30 s), and 72°C (45 s). The PCR products were run in 2% agarose gels (Low-EEO/Multi-Purpose/Molecular Biology Grade Fisher BioReagents™) to confirm absence of DNA amplification. Amplified cDNAs from the DNA-free RNA extracts were prepared using the Ovation RNA-Seq System V2 (Tecan) following manufacturer’s suggestions. All steps through cDNA preparation were completed the same day to avoid freeze/thaw cycles. cDNAs were submitted to the Georgia Genomics and Bioinformatics Core for sequencing using NextSeq 500 PE 150 High Output (Illumina). The cDNA library generated from our control did not contain detectable DNA. It was nonetheless submitted for sequencing, but it failed to generate any sequences that met the minimum length criterion of 300-400 base pairs.

Reads from the 13 metatranscriptome samples collected from sites that metagenomic samples were taken from were mapped to each MAG using the CoverM 0.6.1 (https://github.com/wwood/CoverM) command line tool with the BWA 2.0 aligner81. The CoverM tool automatically concatenated all the MAGs into a single file, and metatranscriptome reads were recruited to MAG contigs, setting the parameter --min-read-percent-identity to 95 and --min-read-aligned-percent to 50. A custom R script was utilized to concatenate all coverM output files into a single file in a matrix format, with each sample representing a column and each row representing total percent of RNA-Seq reads per sample that mapped to a MAG. The output was used in this study for heatmap plotting to examine evidence for activity of the taxa for which were covered MAGs. Metatranscriptome reads were deposited to the National Center for Biotechnology Information Sequence Read Archive under accession numbers SRR22580929-SRR22580947 (Bioproject PRJNA909197).

Cell counts
The sediment sampling for cell counts occurred immediately after core retrieval on the core receiving platform by sub-coring with a sterile, tip-cut 2.5 cc syringe from the center of each freshly cut core section. Approximately 2 cm3 sub-cores were immediately put into tubes containing fixation solution consisting of 8 mL of 3xPBS (Gibco™ PBS, pH 7.4, Fischer) and 5% (v/v) neutralized formalin (Thermo Scientific™ Shandon™ Formal-Fixx™ Neutral Buffered Formalin). If necessary, the mixture was stored at 4°C.

Fixed cells were separated from the slurry using ultrasonication and density gradient centrifugation94. For cell detachment, a 1 mL aliquot of the formalin-fixed sediment slurry was amended with 1.4 mL of 2.5% NaCl, 300 μL of pure methanol, and 300 μL of detergent mix95, 100 mM ethylenediamine tetraacetic acid [EDTA], 100 mM sodium pyrophosphate, 1% [v/v] Tween-80). The mixture was thoroughly shaken for 60 min (Shake Master, Bio Medical Science, Japan), and subsequently sonicated at 160 W for 30 s for 10 cycles (Bioruptor UCD-250HSA; Cosmo Bio, Japan). The detached cells were recovered by centrifugation based on the density difference of microbial cells and sediment particles, which allows collection of microbial cells in a low-density layer. The sample was transferred onto a set of four density layers composed of 30% Nycodenz (1.15 g cm–3), 50% Nycodenz (1.25 g cm–3), 80% Nycodenz (1.42 g cm–3), and 67% sodium polytungstate (2.08 g cm–3). Cells and sediment particles were separated by centrifugation at 10,000 × g for 1 h at 25°C. The light density layer was collected using a 20G needle syringe. The heavy fraction, including precipitated sediment particles, was resuspended with 5 mL of 2.5% NaCl, and centrifuged at 5000 × g for 15 min at 25°C. The supernatant was combined with the previously recovered light density fraction. With the remaining sediment pellet, the density separation was repeated. The sediment was resuspended using 2.1 mL of 2.5% NaCl, 300 μL of methanol, and 300 μL of detergent mix and shaken at 500 rpm for 60 min at 25°C, before the slurry sample was transferred into a fresh centrifugation tube where it was layered onto another density gradient and separated by centrifugation just as before. The light density layer was collected using a 20G needle syringe, and combined with the previously collected light density fraction and supernatant to form a single suspension for cell counting.

For cell enumeration, a 50%-aliquot of the collected cell suspension was passed through a 0.22-μm polycarbonate membrane filter. Cells on the membrane filter were treated with SYBR Green I nucleic acid staining solution (1/40 of the stock concentration of SYBR Green I diluted in Tris-EDTA [TE] buffer). The number of SYBR Green I–stained cells were enumerated either by direct microscopic counts70 or image-based discriminative counts96. For image-based discriminative counting, the Count Nuclei function of the MetaMorph software (Molecular Devices) was used to detect and enumerate microbial cells.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Geller-McGrath et al. used metagenomic and metatranscriptomic approaches to study the environmental distribution, function and activity patterns of bacteria and archaea along thermal, geochemical, and cell count gradients in five different drilling sites of the Guaymas Basin (Gulf of California, Mexico). The authors found that the composition and distribution of MAGs dominated by Chloroflexota and Thermoproteota were determined by biogeochemical parameters in sediment layers with moderate temperature (lower than 45ºC),but become influenced by temperature as it exceeds 45ºC. Overall, MAG genome sizes and diversity decreased with increasing temperature. Additionally, the study revealed that specific microbial populations of Thermoproteota and Hadarchaeota showed an increasing trend in relative abundance and recruitment of reads from the transcriptome with increasing depth temperature.\n\nThe strength of the manuscript is the systematic anlaysis of the biogeochemical parameters along with the genomes and metranscriptomes of the microorganisms in the deep subsurface of hydrothermal regions. However, there are several key shortcomings greatly limited the quality of the MS, including: (1)lacking of enough details on methodologies, especially detailed information about the omics, which lowered the possibilities to evaluate the quality of the manuscript. For example, The concentration, total amount and integrity of the DNA and RNA for omics sequencing? How did the metagenome and metatranscriptome libraries been constructed? What are the sequencing depths? (2) The scientific question of the manuscript is not quite clear (or not specific enough). The major body of the manuscript is mainly descriptive information, and the results and discussion focused too much on scattered details and lack integration and indepth analysis. (3) Additionally, the significance and impact of the findings were not clearly discussed and illustrated, which need to be further elaborated.\n\nSome detailed comments:\nIntroduction: The introduction on biodiversity, activity and functions of microbial communities of the Guaymas Basin, are not sufficient to allow readers to understand the significance and necessarity of this study; and the scientific question and objectives of the study were not clearly stated.\nLine 78-79, this sentence is not clear. what exactly do you mean by \"extent and characteristics\" ?\nLine 92-93, this seems to be the scientific question the MS wanted to address, however, there were too many questions within this sentence, e.g., environmental distribution, genomic potential, adaptation, stress tolerrance. Although each of the above questions was mentioned in the manuscript, none of them were clearly and solidly illustrated. I would suggest the authors to be more focus on one or two points and make in-depth analysis and discussion.\nLine 129, what do you mean by \"incomplete downcore sulfate depletion\"?\nLine 142-143, the authors measured petroleum hydrocarbon, why not invoving them in later biostatistical analysis?\nLine 148-149, why did negative controls have sequenced reads and matched with MAGs?\nLine 194-213, As one of the major conclusion in the manuscript, the influence of envieornmental factors on MAG composition is largely descriptive and lacking of sufficient biostatistical supporting. For example, Line 203-206, the authors claimed differnt environmental parameters that affected the MAG composition in different temperature zone. However, the results were mainly based on the descriptions on nMDS plot, which was lacking of direct statistical support. Also, the authors only described the parameters with positive correlation with the MAG composition and ignored those with negative correlations. In addition, the validity of using nMDS analysis to study the influence of environmental parameters need to be evaluated. Why not using CCA, RDA etc. ?\nLine 215-246, the \"Metagenomic subsurface adaptation\". (1)What are environmental features of studied sites that are special and distinct from other marine habitats? How does the annotated functions provide adaptation advantage for the MAGs, associated with the special environmental conditions of the studied sites? The linkage between the genomic features and specific environmental condition of the study sites is weak. (2) The identified \"adaptation features\" seemed to be commonly present in other marine habitats such as marine sediments, so what are the special adaptation mechanisms for the thermally heated subsurface microorganisms?\nLine 239-240, how common did the related functions in the MAGs (percentage)?\nLine 273-276, this sentence reads quite strange, please modify.\nLine 358 and other locations, what do you mean by \"unique MAGs\"? A definition is needed.\nLine 357-367, the relationship between MAGs recovery and temperature may need to be revisit by considering the sequencing depth of the metagenomes. As the variation on sequencing depth would greatly impact the MAG recovery rate and reads recruitment. For example, if the metagenomes from the hot regions were generally with less sequencing depth (smaller data size), the observed trends of decreasing MAG recovery with temperature would be questionable.\nLine 424-425, More details in library construction are needed. Were they whole genome amplified before sequencing? What were the concentration and total amount of DNA extracted from the samples? What is the sequencing depth? How many Gb of data for each metagenome?\nLine 445-449, Why were there contaminants ? and how could they been sequenced and with so many reads?\nLine 475, not clear. How could abundance of metagenomic reads been normalized as \"transcripts\" per million?\nLine 536-537, Please provide the concentration and integrity of the RNAs extracted. Also, provide details on cDNA libraries construction. Were the RNA samples with enough concentration and total amount for transcriptome sequencing?\nFigure 4, Description is not clear. Are those data based on composition of MAGs or entire metagenomes from each site? Also, the resolution of the plot was low.\nThe methodologies for biogeochemical parameter measurement seemed to be missing?", "reviewer_2": "In the manuscript “Metagenomic Profiles of Archaea and Bacteria within Thermal and Geochemical Gradients of the Guaymas Basin Deep Subsurface” Geller-McGrath and Mara et al present a genomic resolution of microbes and their transcripts along a temperature gradient in the Guyamas Basin. The authors combine geochemical measurements with genome-resolved metagenomics and transcriptomics. While I think this paper is of merit and it definitely has the potential to be appear in Nature Communications, I do have some reservations regarding the presentation of the results:\n1. The manuscript is very technical. While the results are outstanding and interesting per se, I’m missing an interpretation that makes these findings important, e.g. ecosystem functioning. Even the title is very descriptive. So, the story, the meaning of the results need to be developed.\n2. The headers of the results/discussion section should present a major finding rather than be technical and – in the worst case – state a method or simply a phylum name like “Chloroflexota”, which doesn’t tell the reader anything.\n3. Although I do like the correlation of MAG recovery/size to temperature or depth, I’d love to see how the functions in the ecosystem are maintained in the metagenomes (not only binned genomes) and metatranscriptomes (mapped to entire (assembled) metagenomes).\n4. How do the results of the current study compare to “Rapid metabolism fosters microbial survival in the deep, hot subseafloor biosphere” by Beulig et al. Nature Communications volume 13, Article number: 312 (2022)? I’m missing this comparison because it seems that the current study contradicts the findings of Beulig et al.\nMinor comments:\nL 36 and elsewhere in the manuscript: It should be “By contrast, xxx” or “In contrast to XY, xxx” but not “In contrast, xxx”.\nL 95: consider rephrasing “mRNA evidence”.\nL199-203: Have the authors considered using a bioENV to figure out the best combination of environmental variables that described the community composition?\nL336: Have the authors considered genome completeness of the MAGs? Taking the completeness into account when it comes genome sizes might be important but a metric other than CheckM needs to be used, for which markers are equally distributed across the genome or generally takes the unequal distribution into account.\nL369: How can DNA/RNA recovery overcome cell density? I do not understand the ecological context here.\nL431: Have the authors considered combined assemblies of metagenomes and metatranscriptomes to promote recruitment of reads? I’m not saying this is a necessity, but it could improve the mapping frequency. What is the rationale behind co-assembling all the 26 metagenomes rather than individually assembling the metagenomes?\nFigure 1: In-situ should be italics (in figure and in legend)\nFigure 6: MAG recovery can be affected by many biological features (e.g., strain heterogeneity). I think a different metric for diversity should be used, which is based on marker genes like ribosomal protein S3 to identify if temperature and depth have an impact on the recovery of genomes / diversity per se.", "reviewer_3": "In “Metagenomic profiles of archaea and bacteria within thermal and geochemical gradients of the Guaymas Basin deep subsurface,” the authors examine trends in microbial community composition and metabolic potential by examining several MAGs recovered from cores at various points in Guaymas Basin.\nOverall, the manuscript is well-written and the methods are sound. My primary piece of feedback is more of a question: why focus only on MAGs for this study rather than also include results from the full metagenomes? As the authors note, there are limitations and perhaps biases to MAG recovery that result from challenges like poor assembly, poor read recruitment, etc, and it is possible that these limitations may introduce biases with depth, skewing the results. It seems to me that including some results from an analysis of the full metagenome (i.e. a broader survey of community composition using something like mmseqs2, and/or a broader survey of total metabolic potential focusing on relative abundance of key metabolic genes independent of MAGs) would provide broader context for how the whole community changes with depth independent of the MAGs, and would ameliorate potential biases in MAG recovery that may occur with depth.\nIn addition to this, the authors put a lot of information in the supplement, but some parts of the Results/Discussion feel overly broad and/or incomplete. I would suggest adding some details back into the Results/Discussion, especially as it pertains to specific metabolic potential of specific taxa/bins and how that changes with depth. As a further example, Supplementary Figure 2 appears to simply be a copy of Figure 2, with a bit more taxonomic detail. It seems to me that these details should be included in the main Figure 2 and there doesn’t need to be an additional figure in the supplement.\nFurther comments are listed below.\nIntroduction\nLine 48: maybe clarify that this would be downcore in cool/moderate temps, right? These are subsurface sediments away from hydrothermal activity?\nMethods:\nLine 435 What read trimming technique was used?\nThis is just a suggestion, but if the authors are trying to reconstruct metabolic modules and the pipelines they use can’t do archaea, perhaps try KEGGdecoder, which is more agnostic.\nFor the most part, though, the methods make sense and it seems that appropriate controls were run and appropriate precautions were taken regarding contamination.\nResults and Discussion:\nLine 116: remove “at”: “cell counts decreased by four to five orders of magnitude”\nFor the section on Subsurface biogeochemical zonation: a lot of chemistry data is reported in table 1. This is a bit hard to visualize in context of samples—would it be possible to show this as a figure somehow?\nIt is unclear to me what units are in the x-axis of supplementary figure 1.\nFor the section on MAG diversity, distribution and evidence of activity:\nThe authors comment on relative abundances of MAGs in different depths, which is interesting, but why are specific MAGs within a specific phylum more abundant than others? More description of differences between specific MAGs (sort of like a mini pangenome analysis) would help contextualize these results. For example, why are some Chloroflexi at high abundance and not others? Is there information that might contextualize this based on the metabolic pathway completeness analyses?\nAlong the lines of the above comment, a figure showing bin metabolic pathway completeness heatmap to compare with figures 2 and 3 would be helpful. (KEGGDecoder makes heatmaps like these, but I am sure other pipelines do too.)\nLines 171-178: The section on metabolisms, and the mention of CRISPRs, needs to be contextualized, as it is seemingly thrown in there with no further discussion. What can we derive from this?\nFigure 2: For relative abundance, did the authors normalize according to the length of the MAG? The values are reported as what percent of reads per sample mapped to a MAG, so it’s normalized for the number of reads per sample, but in addition to that it seems like MAGs that are longer would have more reads that recruit to them. I checked the Materials and Methods and this isn’t clear—is this normalization calculated as part of CoverM?\nFor the section on metatranscriptome mapping (lines 179-194):\nAs above, this section would benefit from a bit more detail and context. Why are specific MAGs more active than others within specific phyla? Were certain metabolisms more active than others as depth and/or temperature increased? Were these metabolisms specific to certain phyla? How did activity compare with overall DNA abundance?\nIn examining the supplementary material, it seems that a lot of this information is in the supplementary material, but without this context in the main manuscript, the results are overly vague. Please consider perhaps condensing much of the info in the supplementary materials and putting it in the main manuscript.\nAdd “that” to line 219: “widespread genomic features that extend”\nLines 241-248: this is an interesting result; the authors should have the data indicating whether genes for chemotaxis were actually expressed in the MAGs, no? It seems that would be worth reporting here.\nSupplementary Figure 2 seems like it’s exactly figure 2, but with more detail regarding which taxa each row correspond to. It seems like that should be included in figure 2, and supplementary figure 2 could be dispensed with.\nConsider reorganizing the text: most of the phylum-specific MAG description is focused on abundance (DNA), not RNA. One suggestion would be to add the metatranscriptomic reults later, and/or add metatranscriptomic results to phylum-specific descriptions.\nAs I mentioned above, I also think this section would be helped by a metabolic pathway completion map like those produced by KEGGdecoder. Would help visualize differences between these MAGs, visualize metabolic potential.\nRika Anderson"}, "author_response": "REVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nGeller-McGrath et al. used metagenomic and metatranscriptomic approaches to study the environmental distribution, function and activity patterns of bacteria and archaea along thermal, geochemical, and cell count gradients in five different drilling sites of the Guaymas Basin (Gulf of California, Mexico). The authors found that the composition and distribution of MAGs dominated by Chloroflexota and Thermoproteota were determined by biogeochemical parameters in sediment layers with moderate temperature (lower than 45ºC), but become influenced by temperature as it exceeds 45ºC. Overall, MAG genome sizes and diversity decreased with increasing temperature. Additionally, the study revealed that specific microbial populations of Thermoproteota and Hadarchaeota showed an increasing trend in relative abundance and recruitment of reads from the transcriptome with increasing depth temperature.\nResponse: This are indeed the major points of our manuscript, well summarized.\nThe strength of the manuscript is the systematic analysis of the biogeochemical parameters along with the genomes and metranscriptomes of the microorganisms in the deep subsurface of hydrothermal regions. However, there are several key shortcomings greatly limited the quality of the MS, including: (1)lacking of enough details on methodologies, especially detailed information about the omics, which lowered the possibilities to evaluate the quality of the manuscript. For example, The concentration, total amount and integrity of the DNA and RNA for omics sequencing? How did the metagenome and metatranscriptome libraries been constructed? What are the sequencing depths? (2) The scientific question of the manuscript is not quite clear (or not specific enough). The major body of the manuscript is mainly descriptive information, and the results and discussion focused too much on scattered details and lack integration and indepth analysis. (3) Additionally, the significance and impact of the findings were not clearly discussed and illustrated, which need to be further elaborated.\nResponse: Thank you for identifying three major areas of improvement. In response to 1), we are now reporting DNA and RNA yield, and corresponding libraries (before and after QC) in a new table (Supplementary Table 2). We are reporting methodologies in greater depth (lines 474-475; 502-523; 621-625). 2) The scientific rationale for metagenomic exploration of the deep biosphere was explained more explicitly in the revised introduction (new lines 51 to 56), and the manuscript text is revised extensively to reveal the overall scientific perspective underlying the manifold findings that we are reporting (as recommended by other reviewers as well). 3) We have developed the Conclusion section further (now “Conclusion and Outlook”) and added new research perspectives for deep, hot biosphere studies at the interface of sediment/hot basalt intrusions. These chemically complex interfaces are underlying the deep Ringvent sediments where we have consistently found high-temperature archaeal communities, and might constitute a new deep subsurface ecological niche (lines 451-461).\nWith regard to text, the headings and introductions to Results & Discussion sections have been revised, to emphasize the nature of the most significant findings and their context: “Metagenomic features with wide subsurface distribution”, line 245 [formerly “metagenomic adaptations”]; “Characteristics and distribution of dominant bacterial and archaeal groups”, line 292, with subsections “Dominant subsurface bacteria” (line 302), “Dominant subsurface archaea” (line 335), and “Hadarchaeotal genomic features” (line 379). Finally, new headers were included for the last sections “Genome size trends in the deep biosphere” (line 401), and “Temperature impact on MAG recovery” (line 422).\nSome detailed comments:\nIntroduction: The introduction on biodiversity, activity and functions of microbial communities of the Guaymas Basin, are not sufficient to allow readers to understand the significance and necessarity of this study; and the scientific question and objectives of the study were not clearly stated.\nResponse:We have now stated our research rationale explicitly, first by introducing the recent detection of microbial cells and activity in the deep, hot biosphere (Beulig et al. 2022; Heuer et al. 2020), and the lack of further identification that makes it impossible to specify the biological diversity and potential of the different kinds of bacteria and archaea in the deep, hot biosphere. Therefore we note that bacterial and archaeal members of the deep, hot biosphere and their activity need to be identified by sequence-based approaches (new text lines 51 to 56). A new opening line for the next paragraph (line 61) establishes the link to Guaymas Basin as a suitable field site for this research task. The scientific goals are now stated more specifically and accurately in lines 98 to 104.\nLine 78-79, this sentence is not clear. what exactly do you mean by \"extent and characteristics\" ?\nResponse: We have replaced these words by the more specific terms “the spatial extent, diversity and activity” (lines 74-75).\nLine 92-93, this seems to be the scientific question the MS wanted to address, however, there were too many questions within this sentence, e.g., environmental distribution, genomic potential, adaptation, stress tolerance. Although each of the above questions was mentioned in the manuscript, none of them were clearly and solidly illustrated. I would suggest the authors to be more focus on one or two points and make in-depth analysis and discussion.\nResponse: We have re-focused this section on environmental distribution and genomic potential of subsurface microbiota; adaptation and stress tolerance are special aspects of genomic adaptation that will be highlighted at suitable occasions in the manuscript text. We have also inserted an introductory line for this paragraph “These contrasting sites provide an opportunity for a comprehensive analysis of subsurface microbiota” (line 89), to link the drilling sites (described briefly in the previous paragraph) to the scientific questions that are introduced here.\nLine 129, what do you mean by \"incomplete downcore sulfate depletion\"?\nResponse: We have changed this term to “show gradual downcore sulfate consumption (from 27.9 to 18.8 mM)”; the sulfate concentrations in parentheses make it clear that this is not depletion to zero sulfate. (line 131).\nLine 142-143, the authors measured petroleum hydrocarbon, why not involving them in later biostatistical analysis?\nResponse:The problem is that these data have numerous gaps; the high cost of these analyses prevented analyzing a fully matching suite of samples. Thus, we are reporting these data in. Supplementary Table 1 and in Supplemental Figures 1A and 1B, but due to extensive data gaps we cannot include them in the correlational analysis.\nLine 148-149, why did negative controls have sequenced reads and matched with MAGs?\nResponse:\nWhile laboratory extraction (laboratory and kit contamination potential) controls may not produce DNA in quantifiable concentrations (as for our two blank controls), they nearly always produce sequence data, and these same data can show up in the data sets for samples prepared with those kits, particularly when there is very little biomass in the original samples (the “kitome”, Salter et al. 2014). Therefore, taxa detected in negative controls must be removed from sample data. We have three negative controls that produced sequence data, two blank controls (only lab reagents for datasets produced by the Teske and Edgcomb labs using the same kits) and a drilling fluid control (which is surface seawater plus drilling lubricant). The reads recovered for each of the three controls before and after QC are reported in Supplementary Table 2. Originally, contaminants were removed based on taxonomic identity in contamination controls and commonly retrieved suspect taxa in the literature (Salter et al. 2014). As requested, to confirm accurate contamination removal, we mapped control reads to our metagenome assembly and removed all contigs that received mapping with minimum 98% identity over minimum 75% of the read length. Contaminant MAGs were identified in this way (Supplementary Table 3), and they were consistent with published contamination profiles (Salter et al. 2014). This did not remove any of the MAGs discussed in this paper (see lines 524-528, and Supplementary Text paragraph “Accounting for Seawater and Laboratory Contamination”).\nLine 194-213, As one of the major conclusions in the manuscript, the influence of environmental factors on MAG composition is largely descriptive and lacking of sufficient biostatistical supporting. For example, Line 203-206, the authors claimed different environmental parameters that affected the MAG composition in different temperature zone. However, the results were mainly based on the descriptions of nMDS plot, which was lacking of direct statistical support. Also, the authors only described the parameters with positive correlation with the mAG composition and ignored those with negative correlations. In addition, the validity of using nMDS analysis to study the influence of environmental parameters need to be evaluated. Why not using CCA, RDA, etc. ?\nResponse:To comment on the assumption that the nMDS plot lacks direct statistical support, it is important to understand that the plot shows the vectors of environmental and biogeochemical parameters that correlate with MAGs after passing a statistical threshold filter. The environmental factors plotted and shown in the nMDS plot have a p value of less than 0.05 based on permutation tests using the environmental factors and MAG abundance data sets.\nHowever, we agree that running a CCA would also be informative since the CCA approach allows to take into consideration non-linear regression between different parameters, which nMDS does not do. CCA results underscored most of the earlier nMDS results, since (in addition to the central importance of temperature) central geochemical parameters (among others, total sulfide, TOC, TN, pH, salinity) are shown to impact MAG composition for distinct thermal regimes. We are now commenting on both analyses (nMDS and CCA) in the main manuscript (see extended text in lines 217-243), and we include the CCA plot as supplementary Figure 5.\nLine 215-246, the \"Metagenomic subsurface adaptation\". (1) What are environmental features of studied sites that are special and distinct from other marine habitats? How does the annotated functions provide adaptation advantage for the MAGs, associated with the special environmental conditions of the studied sites? The linkage between the genomic features and specific environmental condition of the study sites is weak. (2) The identified \"adaptation features\" seemed to be commonly present in other marine habitats such as marine sediments, so what are the special adaptation mechanisms for the thermally heated subsurface microorganisms?\nResponse: We are discussing widespread genomic features that were found in GB subsurface microbes, and appear to be maintained in this environment; their persistent occurrence indicates that they are useful and relevant (genes involved in the regulation of sporulation, for example). We do not claim that these features are unique for the Guaymas Basin subsurface, and we are reframing this paragraph to clarify expectations (see lines 245-250). The new title “Metagenomic features with wide subsurface distribution” more accurately describes the results that are presented here. Most importantly, we imply that a sizeable portion of these widely distributed and preserved genes does in fact have adaptive value, and confers an advantage to microbes that contain and maintain them. However, other widely distributed genes require more subtle interpretations and we weigh the pros and cons of adaptive explanations while providing perspective on the (slow) pace of subsurface evolution (see especially lines 281-290).\nLine 239-240, how common did the related functions in the MAGs (percentage)?\nResponse note: We have now included percentage numbers for these functions, and integrated them into the extensively revised text of this section (lines 261-273).\nLine 273-276, this sentence reads quite strange, please modify.\nResponse: The entire section and the preceding paragraph were rewritten, the sentence is now subdivided as reads as follows (line 341-344): “Archaeal MAGs were dominated by the Thermoproteota, an archaeal phylum consisting of four major lineages, the Thaumarchaeota, Aigarchaeota, Korarchaeota and Bathyarchaeia (Oren & Garrity, 2021). All 11 Thermoproteota MAGs from the Guaymas subsurface belonged to the class Bathyarchaeia; these MAGs were detected at all examined sites...”\nLine 358 and other locations, what do you mean by \"unique MAGs\"? A definition is needed.\nResponse: “Unique” means that MAGs were dereplicated before further analysis (see new methods text, lines 513-523). After dereplication, the 89 resulting MAGs were distinct from each other, and none of these MAGs were redundant. After this clarification, there is no need for the term “unique” and it was omitted.\nLine 357-367, the relationship between MAGs recovery and temperature may need to be revisited by considering the sequencing depth of the metagenomes. As the variation on sequencing depth would greatly impact the MAG recovery rate and reads recruitment. For example, if the metagenomes from the hot regions were generally with less sequencing depth (smaller data size), the observed trends of decreasing MAG recovery with temperature would be questionable.\nResponse: We now present as Supplementary Table 2 the statistics on sequencing depth for all samples. This shows the number of raw reads obtained per sample, the number of reads post QC, the number of assembled contigs per sample, and the number or contigs retained after removing possible contaminants. We ran ANOVA on the fastQ file sizes/metagenome sequencing depths, and we found no significant differences in metagenome sequencing depth between any of the temperature groups. We don’t think that data interpretation is significantly impacted by differences in data recovery from different depths/temperature regimes.\nLine 424-425, More details in library construction are needed. Were they whole genome amplified before sequencing? What were the concentration and total amount of DNA extracted from the samples? What is the sequencing depth? How many Gb of data for each metagenome?\nResponse: We are adding more methods detail (lines 621-625) in section “RNA extraction, library preparation, sequencing, and mapping of RNA reads to the MAGs”: “Amplified cDNAs from the DNA-free RNA extracts were prepared using the Ovation RNA-Seq System V2 (Tecan) following manufacturer’s suggestions. All steps through cDNA preparation were completed the same day to avoid freeze/thaw cycles.” cDNAs were submitted to the Georgia Genomics and Bioinformatics Core for library preparation and sequencing using NextSeq 500 PE 150 High Output (Illumina).”\nNo whole genome amplification was performed at any point in our study, since the consequences of this amplification step for metagenome composition are just too unpredictable. DNA concentrations are now tabulated in new supplementary table 2 that is also reporting other basic parameters relating to mRNA recovery, and sequencing depth in reads per library.\nLine 445-449, Why were there contaminants ? and how could they have been sequenced and with so many reads?\nResponse: We assume that this comment refers to seawater, drilling and laboratory contamination, not the contamination that is methodologically inherent in the metagenomic reconstruction of MAGs from highly complex natural microbial communities. Deep-sea drilling employs a mixture of seawater and lubricant mud during drilling operations, which carries surface seawater-derived microbial contaminants into the sediment samples. Sample processing on the ship includes steps to remove the surfaces of recovered core samples prior to collection for different analyses. This minimizes the chances of contamination by drilling fluids. Contamination monitoring (using polyfluorinated chemical tracers) was run throughout the drilling operations (see Lever et al. 2006 Geomicrobiology Journal 23:517-530 for details) but it is never perfect. Thus, the mixed microbial community in seawater and drilling lubricant (which can be substantial!) needs to be accounted for by a drilling fluid control.\nThe second type of contaminants are introduced through handling, DNA extraction kits and lab chemicals (the “kitome”, see Salter et al. 2014. BMC Biol. 12:87). This contaminant community is commonly less diverse than the seawater + drilling fluid control, but needs to be accounted for as well, based on blank extractions (only lab kit and chemicals, no sample) and the relevant literature. To be sure, we have run two blank extraction controls, since DNA was extracted in two labs (using identical kits); at USC (Amend lab) and at WHOI (Edgcomb lab).\nThis background on contamination control data handling is now included in a text paragraph and references in Supplementary Materials (“Accounting for Seawater and Laboratory Contamination”), and referenced in the main text (lines 524-528).\nLine 475, not clear. How could abundance of metagenomic reads been normalized as \"transcripts\" per million?\nResponse: In the context of normalizing metagenome reads, the term “transcript” is an unfortunate misnaming, but is currently established in metagenomic analysis and in the relevant bioinformatics literature. “Transcript per million” or TPM is a unit term used in the bioinformatic analysis of both metatranscriptomes and metagenomes; it is unrelated to the process of transcription. We have revised this sentence and the entire paragraph to better describe TPM as a bioinformatics scale parameter (line 553-562). See also our response to a related question from Reviewer 3 about the normalization procedure inherent in CoverM.\nLine 536-537, Please provide the concentration and integrity of the RNAs extracted. Also, provide details on cDNA libraries construction. Were the RNA samples with enough concentration and total amount for transcriptome sequencing?\nResponse: We have added information on RNA and DNA quantity extracted from each sample, and sequencing outcomes for these extracts, to Supplementary Table 2. Regarding the RNA samples, since these samples are extremely difficult to extract RNA from, we elected not to spend any of our RNA for RNA integrity assessments, but to prioritize removal of DNA, and confirmation of removal of all DNA by PCR after DNAse treatment. RNA integrity is likely suboptimal due to the additional steps we introduced in our extraction protocols that included ethanol washes of the sediments. However, numerous other extraction methods and variants that we tried failed to yield RNA in any workable amounts, and the methods used here were the only ones that yielded anything.\nFigure 4, Description is not clear. Are those data based on composition of MAGs or entire metagenomes from each site? Also, the resolution of the plot was low.\nResponse: The data in Figure 4 are based on the composition of MAGs. The figure legend (Figure 4) was revised for clarity and brevity as follows:\n“The nMDS plot depicts the correlation of Guaymas Basin MAG occurrence with in-situ environmental parameters. The directions of the arrows indicate a positive or negative correlation among the environmental parameters with the ordination axes, for statistically significant (p < 0.05) environmental parameters (plot stress: 0.106). Arrow length reflects the strength of correlation between the environmental parameter and MAG occurrence, with longer lines indicating stronger correlations. The samples are color-coded by site, and their temperature regimes are indicated by shape (circles for 0-20°C; triangles for 20-45°C and squares for 45-60°C).”\nWe are providing new figures (nMDS and CCA versions) in the revision, at high resolution.\nWe have modified the figure legend in the CCA version (Supplementary Figure 5).\nThe methodologies for biogeochemical parameter measurement seemed to be missing?\nResponse: The biogeochemical methods are outlined in the Methods chapter of the IODP 385 Expedition report, now included in the reference list as Teske et al. 2021a (lines 474-475).\nReviewer #2 (Remarks to the Author):\nIn the manuscript “Metagenomic Profiles of Archaea and Bacteria within Thermal and Geochemical Gradients of the Guaymas Basin Deep Subsurface” Geller-McGrath and Mara et al present a genomic resolution of microbes and their transcripts along a temperature gradient in the Guaymas Basin. The authors combine geochemical measurements with genome-resolved metagenomics and transcriptomics. While I think this paper is of merit and it definitely has the potential to be appear in Nature Communications, I do have some reservations regarding the presentation of the results:\n1. The manuscript is very technical. While the results are outstanding and interesting per se, I’m missing an interpretation that makes these findings important, e.g. ecosystem functioning. Even the title is very descriptive. So, the story, the meaning of the results needs to be developed.\nResponse:Thank you for pointing out the need to translate results per se (the technical aspect) into meaningful statements on ecosystem function. Reviewer 1 made a similar comment that the significance was not well elaborated. To emphasize the relevance of our results for ecosystem functioning, we have amended our manuscript as follows:\nWe elaborate the paragraph on subsurface biogeochemical zonation and add context on the observed trends in major biogeochemical parameters, as well as site-specific observations (lines 128-162); Table 1 was checked and expanded to include TN values. The revised paragraph sets the stage for MAG occurrence patterns and MAG-based ecosystem inferences in the Guaymas subsurface.\nWe are adding an expanded transcript section, to improve our presentation of MAG expression within the subsurface ecosystem (line numbers 200-215). Interestingly, MAGs that recruit transcriptomes at warm and hot temperatures often belong to order-level lineages from warm sulfur springs, sulfidic warm water, hydrothermal chimneys and sediments.\nWe re-emphasize our observation that at moderate and warm temperatures the microbial community reacts to chemistry (meaning, it is subject to ecological interactions), whereas in hot sediments chemical factors fade away and MAGs relate to temperature as primary factor (lines 217-243).\n2. The headers of the results/discussion section should present a major finding rather than be technical and – in the worst case – state a method or simply a phylum name like “Chloroflexota”, which doesn’t tell the reader anything.\nResponse:The phylum names are replaced with more general and informative headings, “dominant subsurface bacteria” or “dominant subsurface archaea”, etc. These sections have now more informative introductions that emphasize the ecophysiological relevance of these dominant bacterial and archaeal groups. The introductory lines (now 292-301) are extensively rewritten, followed by suitably modified text in lines 302ff. The section on subsurface Chloroflexota was fleshed out with a new paragraph on their metabolic potential (lines 317-334). Interestingly, MAGs of the dominant bacterial and archaeal phyla in warm and hot sediments are related to order-level lineages that were previously obtained from warm, sulfidic, or hydrothermal environments (lines 200-217), not cold benthic sediments.\nWe have also replaced the technical headings of sections on MAG genome size and thermal MAG diversity limits with more informative headings (“Genome size trends in the deep biosphere” [line 401] and “Temperature impact on MAG recovery” [line 422]).\n3. Although I do like the correlation of MAG recovery/size to temperature or depth, I’d love to see how the functions in the ecosystem are maintained in the metagenomes (not only binned genomes) and metatranscriptomes (mapped to entire (assembled) metagenomes).\nResponse:A full presentation of functional and metabolic analyses of these metagenomes is outside of the scope of this paper, especially since a graduate student in another lab who has joined us on the Guaymas drilling expedition is working on such an analysis. However, we have followed the reviewer’s suggestion and have performed a full metagenome analysis, and we now include a summary figure of KEGG categories for the interested reader as Supplementary Figure S4 (with matching Supplementary Table S9). Many functional categories are well maintained, even in the deepest samples for which we could obtain metagenomes. We note the similarity of this overall KEGG category figure to the Chloroflexota version (Supplementary Figure S3), a possible consequence of high Chloroflexota representation in the subsurface microbiome (lines 331-334).\n4. How do the results of the current study compare to “Rapid metabolism fosters microbial survival in the deep, hot subseafloor biosphere” by Beulig et al. Nature Communications volume 13, Article number: 312 (2022)? I’m missing this comparison because it seems that the current study contradicts the findings of Beulig et al.\nResponse:The metagenome profiles in this study (dominated by Hadarachaeota at the high temperature end) approach in-situ temperatures of approx. 60C, whereas Beulig et al. 2022, and the related study by Heuer et al. 2020 report microbial activities (without sequence data) at extreme temperatures up to 120C. The difference is easily understood by the obligatory DNA demands for metagenomic analysis; generally at least 0.5 ng good-quality DNA is required for metagenomic library construction (while avoiding whole genome amplification due to its uncontrolled biases). The very small cell numbers that were observed by Heuer and Beulig in extremely hot subsurface sediments (occasional peaks of ca. 100 cells/cm3 sediment) are not sufficient for metagenomic detection. Without belaboring the point, we are referring to this issue (lines 53 to 56). We are currently finalizing a separate methods manuscript explicitly on this topic.\nMinor comments:\nL 36 and elsewhere in the manuscript: It should be “By contrast, xxx” or “In contrast to XY, xxx” but not “In contrast, xxx”.\nResponse: Corrected as suggested (line 36).\nL 95: consider rephrasing “mRNA evidence”.\nResponse: what about this rephrased sentence “...We also provide evidence for the activity of specific bacterial and archaeal lineages by mRNA transcript mapping on bacterial and archaeal MAGs.” (lines 93-95)\nL199-203: Have the authors considered using a bioENV to figure out the best combination of environmental variables that described the community composition?\nResponse: We considered bioENV but did not elect this type of analysis that is based on Euclidean distances. Instead, we used different multivariate analyses (direct and indirect) as per Reviewer 1’s request. We believe that multivariate analyses such as Canonical Correspondence Analysis (CCA) that account for compositional variances are appropriate.\nL336: Have the authors considered genome completeness of the MAGs? Taking the completeness into account when it comes genome sizes might be important but a metric other than CheckM needs to be used, for which markers are equally distributed across the genome or generally takes the unequal distribution into account.\nResponse: We used CheckM2 to estimate completeness, as indicated in the methods. We agree that CheckM has its limitations because the marker genes it looks for are not equally distributed across the genome. CheckM2, by comparison to CheckM, uses machine learning models trained on genome properties including genome length, number of coding sequences, individual amino acid counts, as well as annotation of predicted proteins using KEGG (in the context of KEGG modules and pathways) (Chklovski et al. 2023). A literature search for alternate methods confirmed that CheckM2 is currently the gold standard for genome completeness assessments, acknowledged even by authors who are introducing alternate bioinformatics tools (often fine-tuned for specialized microbial groups).\nL369: How can DNA/RNA recovery overcome cell density? I do not understand the ecological context here.\nResponse: We wanted to say that improved DNA/RNA recovery can compensate the limitations of low cell numbers (and thus, DNA/RNA content). “While improved DNA and RNA recovery could potentially compensate for declining downcore cell density, ...”(line 434)\nL431: Have the authors considered combined assemblies of metagenomes and metatranscriptomes to promote recruitment of reads? I’m not saying this is a necessity, but it could improve the mapping frequency. What is the rationale behind co-assembling all the 26 metagenomes rather than individually assembling the metagenomes?\nResponse:This is an interesting question; there is virtually no limit to the possibilities of different types of analyses one could do with such data. We did not consider a co-assembly of transcriptomes and metagenomes in this study, but this could be attempted separately, as we believe adding additional analyses would overburden our current paper that is already “criticized” for being too complex.\nRegarding the metagenomes, 16S rRNA gene profiles (bacterial data not shown in this paper, however archaeal 16S data are presented in Mara et al. 2023 in press, The ISME J) indicate a gradient of community composition with depth (as opposed to distinctly separate communities in every sample), thus we suspected that metagenomes would similarly exhibit a gradient of overlap in composition. Our results support this. Given the complexity of sediment communities we felt that a co-assembly would improve the yield of MAGs. Higher metagenome sequencing depth obtained by co-assembly can yield a more robust assembly that captures more of the diversity in the system, and it facilitates the comparison across samples by giving us one reference assembly to use. It can also improve our ability to recover genomes from metagenomes through differential coverage.\nFigure 1: In-situ should be italics (in figure and in legend)\nResponse: done as suggested.\nFigure 6: MAG recovery can be affected by many biological features (e.g., strain heterogeneity). I think a different metric for diversity should be used, which is based on marker genes like ribosomal protein S3 to identify if temperature and depth have an impact on the recovery of genomes / diversity per se.\nResponse: We agree that MAG diversity should never be interpreted as portraying the complete picture of phylogenetic diversity, and we also never claim that our MAGs represent total diversity. The objective of analysis of assembled MAGs in this study is to learn more about the genomic potential of the taxa whose MAGs are recovered. A separate paper on marker gene profiling of the bacterial and archaeal community is currently being drafted in our group.\nReviewer #3 (Remarks to the Author):\nIn “Metagenomic profiles of archaea and bacteria within thermal and geochemical gradients of the Guaymas Basin deep subsurface,” the authors examine trends in microbial community composition and metabolic potential by examining several MAGs recovered from cores at various points in Guaymas Basin.\nOverall, the manuscript is well-written and the methods are sound. My primary piece of feedback is more of a question: why focus only on MAGs for this study rather than also include results from the full metagenomes? As the authors note, there are limitations and perhaps biases to MAG recovery that result from challenges like poor assembly, poor read recruitment, etc, and it is possible that these limitations may introduce biases with depth, skewing the results.\nIt seems to me that including some results from an analysis of the full metagenome (i.e. a broader survey of community composition using something like mmseqs2, and/or a broader survey of total metabolic potential focusing on relative abundance of key metabolic genes independent of MAGs) would provide broader context for how the whole community changes with depth independent of the MAGs, and would ameliorate potential biases in MAG recovery that may occur with depth.\nResponse: Thank you for your positive and thoughtful evaluation of our efforts, and your suggestions to develop this study further. Examining the total metabolic potential of the subsurface community in greater depth is certainly a possible research direction, but as indicated above in a response to another reviewer, we have not taken this path for two reasons: 1) a total metagenomic survey is already being done by a graduate student who sailed on IODP expedition 385 (although on a somewhat different sample set), and we have deliberately coordinated our manuscripts to avoid overlap. This division of labor goes back to the Guaymas drilling cruise in 2019 when this was agreed upon. 2) In this paper, we really want to emphasize the phylogenomic dimension, and examine how microorganisms and their genomic capabilities change downcore; necessarily this means working with MAGs. On their own, core metabolic genes are affected too much by lateral gene transfer.\nNonetheless, also in response to reviewer 2, we have added a KEGG category re-analysis of the metagenomes for all our samples, and we present the resulting overview figure of KEGG functional gene categories for the interested reader (Chloroflexota metagenome in Supplemental Figure 3, and complete microbial community metagenome in Supplemental Figure 4), to provide an idea of the overall functional potential and its downcore changes within our sample set (lines 317-334). Tables of the KEGG abundance results is also added to the supplementals (Supplementary Tables 8 and 9).\nIn addition to this, the authors put a lot of information in the supplement, but some parts of the Results/Discussion feel overly broad and/or incomplete. I would suggest adding some details back into the Results/Discussion, especially as it pertains to specific metabolic potential of specific taxa/bins and how that changes with depth. As a further example, Supplementary Figure 2 appears to simply be a copy of Figure 2, with a bit more taxonomic detail. It seems to me that these details should be included in the main Figure 2 and there doesn’t need to be an additional figure in the supplement.\nResponse: In writing this manuscript, we have explored numerous possibilities, starting with minimal supplements and a very lengthy manuscript, then turning around to a rigorously pruned manuscript, and finally crafting a compromise where significant text sections were re-integrated into the manuscript (the submitted version). It is certainly possible to re-insert further supplemental materials into the manuscript, yet that another reviewer was unhappy about the relative abundance of descriptive material on various microbial lineages. To solve this dilemma within the context of a manuscript that was growing considerably during revision (ca. 1000 words), we have constructed more informative linkages between the main manuscript and the supplements (beginning with lines 178-188). In the end, we hope to have crafted a compromise that is acceptable for readers who love microbial detail, and for those who tend to skip this over.\nAbout Figure 2, we are inserting the fully labeled version into the manuscript, and remove the supplemental version.\nFurther comments are listed below.\nIntroduction\nLine 48: maybe clarify that this would be downcore in cool/moderate temps, right? These are subsurface sediments away from hydrothermal activity?\nResponse: We have modified the sentence as follows: “As microbial communities in cool, relatively shallow subsurface sediments transition into more deeply buried and increasingly warm and finally hot sediments, it should be possible to track how subsurface bacteria and archaea react to these gradually harsher regimes downcore on the levels of cellular activity and community change.” (lines 44-48)\nYes, the sediments are not hydrothermal sediments in the usual sense (poorly consolidated sediments that are permeated by pulsating hydrothermal fluids); they are consolidated subsurface sediments, heated by the strong geothermal heatflow that is pervasive in Guaymas Basin (Neumann et al. 2023, in ref. list).\nMethods:\nLine 435 What read trimming technique was used?\nResponse: thank you for pointing out this methods gap. We have added the following sentence and references to the methods: “Before assembly, reads were trimmed for quality and adapters removed using Trimmomatic v0.39 (Bolger et al. 2014) (parameters: leading:20; trailing:20; sliding window: 0-24; min length 50). The quality of reads was verified with FastQC v0.11.9 (Andrews et al. 2012).” (lines 505-508).\nReferences: Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014 Aug 1;30(15):2114–20.\nAndrews S, Krueger F, Segonds-Pichon A, Biggins L, Krueger C, Wingett S. 2012. FastQC: a quality control tool for high throughout sequence data. Babraham Institute, UK.\n:http://www.bioinformatics.babraham.ac.uk/projects/fastqc.\nThis is just a suggestion, but if the authors are trying to reconstruct metabolic modules and the pipelines they use can’t do archaea, perhaps try KEGGdecoder, which is more agnostic.\nFor the most part, though, the methods make sense and it seems that appropriate controls were run and appropriate precautions were taken regarding contamination.\nResponse: Thank you for the thoughtful suggestion, it would be a natural follow-up to our metagenomic analyses. However, if we reconstructed metabolic modules and emphasized them in our paper, we would infringe on the research (and the manuscript) of a graduate student who sailed with us, and potentially cause problems at a very vulnerable career stage.\nHowever, the implicit suggestion to at least sketch out a physiological downcore profile makes good sense, and we include an overview of functional gene potential within all samples based on KEGG categories as new Supplementary Figure 4, and a corresponding overview of the dominant bacterial phylum Chloroflexota as Supplementary Figure 3. The KEGG abundance data underlying these supplementary heatmaps are also provided as Supplementary tables 8 and 9.\nResults and Discussion:\nLine 116: remove “at”: “cell counts decreased by four to five orders of magnitude”\nResponse: good catch, corrected.\nFor the section on Subsurface biogeochemical zonation: a lot of chemistry data is reported in table 1. This is a bit hard to visualize in context of samples—would it be possible to show this as a figure somehow?\nResponse: Conceptually, such a figure would combine data from different sampling sites and different gradients, in very rough resolution of tens of meters. Individual data points taken out of site-specific context (see figures in IODP site chapters) would be hard to follow. Each datapoint from a chemical profile fits ultimately into the context of its own site and chemical species, as plotted in the site chapters (which we reference).\nThe contextualization of these environmental data is also presented in the nMDS plot (Figure 4) and in the CCA plot (Supplementary Figure 5). In fact, the purpose of these plots is to extract from the table the factors that are driving MAG distribution.\nYet, to help to track the chemistry in the context of the samples, we have rearranged Table 1 to show each site in downcore sequence.\nIt is unclear to me what units are in the x-axis of supplementary figure 1.\nResponse: Good catch, this is temperature in Celsius (the C is cut off by mistake). We are providing a corrected figure, now with panels 1A and 1B for total petroleum hydrocarbons and saturated hydrocarbons, respectively.\nFor the section on MAG diversity, distribution and evidence of activity:\nThe authors comment on relative abundances of MAGs in different depths, which is interesting, but why are specific MAGs within a specific phylum more abundant than others? More description of differences between specific MAGs (sort of like a mini pangenome analysis) would help contextualize these results. For example, why are some Chloroflexi at high abundance and not others? Is there information that might contextualize this based on the metabolic pathway completeness analyses?\nResponse: This is an excellent comment, and this kind of comparison of genome content for presumably abundant taxa vs rarer taxa would frame a potential follow-up study. For this study, we wanted to dig into this a bit for the Chloroflexota, which appear to be one of the most successful and diverse groups present in our MAGs and in the marine sedimentary subsurface in general (Fincker et al. 2020). We have created a new KEGG category figure that shows KEGG category abundances for all Chloroflexota in a heatmap format (Supplementary Figure 3, and matching abundance table (Supplementary table 8) (text lines 317-334). The figure shows that the Chloroflexota are remarkably consistent in their downcore metagenomic profiles across different samples, and are gradually reduced only at ca. 40C and higher. If the editor wishes to approve an additional main text figure, it can be moved into the main text.\nAlong the lines of the above comment, a figure showing bin metabolic pathway completeness heatmap to compare with figures 2 and 3 would be helpful. (KEGGDecoder makes heatmaps like these, but I am sure other pipelines do too.)\nResponse: We are providing the underlying KEGG table for Chloroflexi MAGs (in excel format) as Supplemental Table 8, which matches the heatmap figure (Supplementary Figure 3). This large table provides access to the abundance of relevant Chloroflexi pathways and gene categories from the Guaymas subsurface, a useful service for the (considerable) group of Chloroflexi researchers who might want to dig unto these data further.\nLines 171-178: The section on metabolisms, and the mention of CRISPRs, needs to be contextualized, as it is seemingly thrown in there with no further discussion. What can we derive from this?\nResponse: This information is included and discussed in more detail in the supplementary information, which is now referenced explicitly in expanded lines of the manuscript text (lines 179-188). Mentions of CRISPR detection are now limited to the supplement, where more context is available.\nFigure 2: For relative abundance, did the authors normalize according to the length of the MAG? The values are reported as what percent of reads per sample mapped to a MAG, so it’s normalized for the number of reads per sample, but in addition to that it seems like MAGs that are longer would have more reads that recruit to them. I checked the Materials and Methods and this isn’t clear—is this normalization calculated as part of CoverM?\nResponse: This question (and a related comment from Reviewer 1 about normalization units) prompted us to check the inner workings of normalization in CoverM with the developer. It turned out that CoverM's \"relative abundance\" metric (that we used when mapping metagenomic and metatranscriptomic reads to our MAGs) accounts for genome size by using the mean coverage of the genomes (as opposed to the total coverage). For example, if we had two MAGs that were equally abundant, but one MAG was 1Mbp in length and the other was 2Mbp, double the reads would map to the MAG that was 2Mpb in length. To account for this, CoverM compares the proportion of mean coverages of the input set of MAGs in its calculation of relative abundance, which would be the same for these two genomes in this example. Its calculation therefore not only accounts for total metagenomic/metatranscriptomic library size, but genome length as well.\nHowever, to test for potential effects of different normalization approaches, we also normalized reads mapped to kilobase of genome, per million reads (the rkpm normalization option in CoverM, see table of normalization options in https://github.com/wwood/CoverM\\#calculation-methods.\nUsing this alternate normalization, Figures 2 and 3 looked practically indistinguishable for the original versions; thus we keep the current (original) versions, after improving them with taxonomy annotation.\nFor the section on metatranscriptome mapping (lines 179-194) (now starting at line 200):\nAs above, this section would benefit from a bit more detail and context. Why are specific MAGs more active than others within specific phyla? Were certain metabolisms more active than others as depth and/or temperature increased? Were these metabolisms specific to certain phyla? How did activity compare with overall DNA abundance?\nResponse: Prompted by this comment, we re-examined Figure 3, which shows the transcriptionally active MAGs. Interestingly, it turned out that most bacterial and archaeal MAGs that were actively transcribed at warm and hot temperatures belonged to lineages that were previously detected in warm sulfidic springs, in hydrothermal carbonate chimneys, and even in hydrothermal sediments of Guaymas Basin. We recognize that we have definitely missed some opportunities here, and are now providing phylogenetic and environmental context for this section (lines 200-215). For easy sample identification, the samples of Figure 3 are taxonomically annotated (on the y-axis).\nIn examining the supplementary material, it seems that a lot of this information is in the supplementary material, but without this context in the main manuscript, the results are overly vague. Please consider perhaps condensing much of the info in the supplementary materials and putting it in the main manuscript.\nResponse: We have added more specific material to the main manuscript: a sharper scientific rationale in the introduction; revised Figures 2 and 3 fully annotated with MAG taxonomy; revised Figure 4 (and matching Supplementary Figure 5) with updated geochemistry correlations; new discussion of transcriptomics in warm and hot samples; and numerous text emendations and revisions throughout the manuscript. With all these additions, we are now running over the 5000-word limit. However, we have inserted more supplement references (text, figure and table links) while eliminating unnecessary duplications in the main text (see lines 179-188).\nAdd “that” to line 219: “widespread genomic features that extend”\nResponse: Good, done\nLines 241-248: this is an interesting result; the authors should have the data indicating whether genes for chemotaxis were actually expressed in the MAGs, no? It seems that would be worth reporting here.\nResponse: Yes, the chemotaxis genes and motility genes are expressed primarily between 0.8 and 36.8 mbsf, and this is reported in a separate paper in press (Mara et al. ISME J). We have inserted a new paragraph where we highlight the slow pace of evolution in the subsurface, and the potential of microorganisms to retain chemotaxis and motility genes in the subsurface (Lines 272-290).\nSupplementary Figure 2 seems like it’s exactly figure 2, but with more detail regarding which taxa each row correspond to. It seems like that should be included in figure 2, and supplementary figure 2 could be dispensed with.\nResponse: We are reorganizing our figures; Figure 2 is now fully labeled with taxonomy information on the y-axis, and then the supplementary version is no longer needed. Figure 3 is also fully labeled for easier reference.\nConsider reorganizing the text: most of the phylum-specific MAG description is focused on abundance (DNA), not RNA. One suggestion would be to add the metatranscriptomic results later, and/or add metatranscriptomic results to phylum-specific descriptions.\nResponse: In this manuscript, metatranscriptomic results are essentially used to validate DNA and MAG results by transcript mapping, something that was not possible in our separate overview manuscript on subsurface metatranscriptomics in Guaymas Basin (Mara et al. ISME J., in press). We make this connection more explicit so that readers will not miss it (lines 189-191): “To determine any intra-phylum differences in metabolic activity, we mapped reads of the Guaymas Basin subsurface metatranscriptome (Mara et al., in press) to our recovered MAGs, for samples collected at the same sites (Figure 3).” To avoid overlap with our own paper (Mara et al., in press) we limit the transcriptomics results to validate MAG activity, and for this purpose they are introduced within the section “MAG diversity, distribution, and evidence of activity” (lines 200-215).\nAs I mentioned above, I also think this section would be helped by a metabolic pathway completion map like those produced by KEGGdecoder. Would help visualize differences between these MAGs, visualize metabolic potential.\nResponse: We would like to keep the phylogenomic gradient (MAG diversity changes downcore) in the center of our manuscript, and avoid overlap with the ongoing metabolism-centered study of our graduate student colleague. However, for overview purposes we are now including a KEGG category analysis of the metagenome samples as Supplemental Figure 4, and its Chloroflexota-focused equivalent as Suppl. Figure 3.\nRika Anderson\nReviewer #1 (Remarks to the Author):\nThe authors made great effort to modify the manuscript and my concerns on methodologies are largely elimnated. However, some shortcomings still remains (1)The major body of the manuscript is mainly descriptive information, and the results and discussion focused too much on scattered details and lack sufficient integration and indepth analysis. (2) the significance and impact of the findings were not clearly discussed and illustrated.\nThe metagenomic and metagenomic data sets from such extreme habitat on Earth are undoubtably with great values, and the discovery of temperature regulated distribution of microbial communities and preference of different taxa in different temperature are very interesting. However, the manuscript still need to be substatially re-organized, making good summarization, integration and discussion of all geochemical, distribution, as well as functional features. The hypotheisis ' that the Guaymas Basin subsurface community is not a random assemblage, but reveals phylogenetic and functional structure that can be tracked downcore. ' is a good and attractive point. Unfortunately, current version of MS focused mainly on phylogenetic structure while the patterns of functional characteristics along the temperature gradient is not clearly shown. I suggest the authors to enhance the analysis on changes of functional composition at different temperature zone, revealing the possible survival strategy and ecological/biogeochemical effects of microbes at different temperature zone, and if possible, find out the genomic features for microbes to adapt to high temperature. Correspondingly, at least one figure showed be included in the maintext to show the functional features of the microorganisms. Figures showing similar kind of data can be integrated, for example, Figure 2 and 3 can be integrated into one.\nPlease note that \"enhance the analysis of functional characteristics\" not means putting more details of functional annotations in the main text, but to make in-depth and integrative anlaysis on a particular question such as \"changes of functional characteristics along the temperature gradient and their implications\".\nWould it be better to show the changes of MAG recovery (i.e., MAG diversity) and MAGs composition along the temperature gradient first, then shown the changes of genomic features( Genome size, and functional characteristics composition) to reveal the functional and metabolic basis for observed species patterns? (This is my personal feeling about the logic, just for communication)\nSome other comments:\nLine 52, why would \"sequence-based identification\" be so important?\nLine 55-57, why is it important to \"learn more about bacterial and archaeal communities of the deep, hot biosphere from a genomic perspective\"?\nThe authors provided details of DNA/RNA extractions in table S2. I found that many samples were with low DNA/RNA concentration (especially for the negative controls), how can they be successfuly sequenced without amplification (low quatities of nucleic acid usually lead to failure in library construction)? What kit was used for metagenomic library construction? Was it the same kit for library construction for the other samples?\nLine 245-399, this part is very long but I feel it can be organized and analyzed in a better way. A lot of content here are repeating the distribution of different species.\nReviewer #2 (Remarks to the Author):\nThe authors have - though many times referring to other studies that are in preparation / accepted elsewhere - all my questions. The only remaining question that I have relates to Figure 6. The authors perform a correlation analysis with unique MAGs. However, there might be other explanation than temperature and depth for a decline / increase in recovered MAGs. Can the authors please also test alternate explanations and correlate the number of unique MAGs with sequencing depth (no of sequenced bps per sample) or DNA yield (from DNA extractions per sample)? Lastly a FDR-correction of the p-values retrieved from these correlations would be appropriate.\nI congratulate the authors on this great study.\nReviewer #3 (Remarks to the Author):\nThe authors have sufficiently and thoughtfully addressed my comments, and I feel that the paper is now suitable for publication.\nReviewer responses\nReviewer #1 (Remarks to the Author):\nThe authors made great effort to modify the manuscript and my concerns on methodologies are largely elimnated. However, some shortcomings still remains (1)The major body of the manuscript is mainly descriptive information, and the results and discussion focused too much on scattered details and lack sufficient integration and indepth analysis. (2) the significance and impact of the findings were not clearly discussed and illustrated.\nResponse: We respectfully disagree on both points. The manuscript has been carefully calibrated to strike an acceptable balance between emphasizing the overall scientific message and adding illustrative detail, first during repeated rewrites among the author team, and then during the review process. Colleagues, editors and reviewers clearly understand the general significance and the scientific implications of our manuscript.\nThe metagenomic and metagenomic data sets from such extreme habitat on Earth are undoubtably with great values, and the discovery of temperature regulated distribution of microbial communities and preference of different taxa in different temperature are very interesting. However, the manuscript still need to be substatially re-organized, making good summarization, integration and discussion of all geochemical, distribution, as well as functional features. The hypothesis ' that the Guaymas Basin subsurface community is not a random assemblage, but reveals phylogenetic and functional structure that can be tracked downcore. ' is a good and attractive point. Unfortunately, current version of MS focused mainly on phylogenetic structure while the patterns of functional characteristics along the temperature gradient is not clearly shown. I suggest the authors to enhance the analysis on changes of functional composition at different temperature zone, revealing the possible survival strategy and ecological/biogeochemical effects of microbes at different temperature zone, and if possible, find out the genomic features for microbes to adapt to high temperature. Correspondingly, at least one figure showed be included in the maintext to show the functional features of the microorganisms. Figures showing similar kind of data can be integrated, for example, Figure 2 and 3 can be integrated into one.\nResponse: The fundamental problem in a paper of high complexity is to choose a particular perspective that maximizes insight and understanding, and we have chosen phylogeny in first place and function in second place. The reason is that phylogenetic composition changes downcore more significantly and more visibly than functional potential. Supplementary figures S3 and 4, the KEGG heat maps of functional gene categories, show most functional categories initially persist and then decline in a similar manner as depth and temperature increase, reminiscent of similar observations in surficial hydrothermal sediments (Su et al., 2023. AEM, 10.1128/aem.00018-23). So, instead of arguing from a functional perspective first and foremost, we anchor our manuscript in MAG phylogenetic diversity trends, and then link functional inference to noteworthy MAGs. Of course, the interested reader will find plenty of functional information in the text and tables, especially in the supplements. This is a goldmine that is open for further exploration.\nFigures 2 and 3 show results for metagenomic and transcriptomic analyses, respectively. Each of these large figures is already near the limit of complexity and readability (for example, regarding taxon labels on the y-axis). Apart from the unsolved question how to “combine” these figures, a combination would produce an overly large figure that does not fit on a page and (if reduction in scale is attempted) becomes impossible to decipher.\nPlease note that \"enhance the analysis of functional characteristics\" not means putting more details of functional annotations in the main text, but to make in-depth and integrative analysis on a particular question such as \"changes of functional characteristics along the temperature gradient and their implications\".\nResponse: This is again a matter of perspective. In our strategy, the “integrative analysis” focuses on MAG diversity first, and from there it is possible to infer functional diversity. This strategy has the advantage of greater specificity; it is more accurate to argue from the organisms (represented by MAGs, see Figures 2 and 3) towards their potential ecophysiological function, than to start with function and trying to split up widely shared and distributed functionalities (see Supplementary Figures 3 and 4) into numerous and diverse MAGs that are participating.\nWould it be better to show the changes of MAG recovery (i.e., MAG diversity) and MAGs composition along the temperature gradient first, then shown the changes of genomic features( Genome size, and functional characteristics composition) to reveal the functional and metabolic basis for observed species patterns? (This is my personal feeling about the logic, just for communication)\nResponse: We have placed changes in MAG recovery last, in order to highlight the limits of MAG detection at 60C, and to showcase again the crucial role of temperature that delimits the deep subsurface biosphere. Ending this study with a focus on temperature limits also facilitates comparison with other deep biosphere studies (which provide temperature information but usually do not contain information on MAG genome size trends).\nSome other comments:\nLine 52, why would \"sequence-based identification\" be so important?\nResponse: Sequence-based identification provides a phylogenetic and physiological identity to microbial populations that are otherwise only known from cell counts.\nLine 55-57, why is it important to \"learn more about bacterial and archaeal communities of the deep, hot biosphere from a genomic perspective\"?\nResponse: Since these deep subsurface communities are resistant to cultivation and lab-based study of isolates, genomic analysis is (for now) the best way to study them. The next generations of deep subsurface microbiologists will hopefully have more cultures at their disposal.\nThe authors provided details of DNA/RNA extractions in table S2. I found that many samples were with low DNA/RNA concentration (especially for the negative controls), how can they be successfully sequenced without amplification (low quantities of nucleic acid usually lead to failure in library construction)? What kit was used for metagenomic library construction? Was it the same kit for library construction for the other samples?\nResponse: Our sequencing libraries were prepared at the Sequencing center of the University of Delaware, which specializes in successful library construction from low-DNA samples. We are discussing DNA yields and methodological requirements in an upcoming paper that has been submitted to the Results volume of the IODP 385 Expedition reports. The low RNA concentrations were something that we expected because the Guaymas subsurface sediments were hydrocarbon-rich. We constructed the cDNA libraries from total RNA using the Ovation RNA-Seq System V2 kit (NuGEN) which uses picograms of template RNA. As in the case of DNA, the cDNA libraries were sent to a sequencing facility (Georgia Genomics and Bioinformatics Core facility) which also specializes in challenging, low biomass samples.\nLine 245-399, this part is very long but I feel it can be organized and analyzed in a better way. A lot of content here are repeating the distribution of different species.\nResponse: In the re-revised manuscript, text in these lines is subdivided into contrasting sections that are thematically distinct, while progressing logically from shared genomic features of the subsurface community towards more specific content on particular bacteria and archaea: “Metagenomic features with wide subsurface distribution” (lines 236 to 279), “Characteristics and distribution of dominant bacterial and archaeal groups” (line 281 to 289), “Dominant subsurface bacteria” (lines 290 to 321), “Dominant subsurface archaea” (lines 322-363), and “Hadarchaoetal genomic features” (lines 363 to 382). Please note that these sections do not discuss “species” but MAGs, and major taxonomic groups.\nReviewer #2 (Remarks to the Author):\nThe authors have - though many times referring to other studies that are in preparation / accepted elsewhere - all my questions. The only remaining question that I have relates to Figure 6. The authors perform a correlation analysis with unique MAGs. However, there might be other explanation than temperature and depth for a decline / increase in recovered MAGs. Can the authors please also test alternate explanations and correlate the number of unique MAGs with sequencing depth (no of sequenced bps per sample) or DNA yield (from DNA extractions per sample)? Lastly a FDR-correction of the p-values retrieved from these correlations would be appropriate.\nI congratulate the authors on this great study.\nResponse: Alternative interpretations were tested using the information in Supplementary Table 2, which contains the sequencing depth for all samples (number of raw reads obtained per sample, the number of reads post QC, the number of assembled contigs per sample, and the number of contigs retained after removing possible contaminants). We ran ANOVA on the fastQ file sizes/metagenome sequencing depths, and found no significant differences in metagenome sequencing depth between any of the three temperature groups. Thus, we don’t think that data interpretation is significantly impacted by differences in sequence recovery from different depths/temperature regimes.\nHowever, DNA yield is another matter, since declining DNA recovery (a proxy for declining microbial biomass and cell numbers) will intrinsically limit the diversity of MAGs that can be recovered. At some point, DNA “extinction” will stop the recovery of MAGs. This issue is now acknowledged explicitly by rewriting the introductory sentence in the paragraph “Temperature impact on MAG recovery”:\n“The environmental stresses that increasingly exclude microbial lineages, reduce genome size and reduce overall microbial population size (and thus, quantity of recovered DNA) are reflected in decreased recovery of MAGs in warmer and deeper samples from all sites (Figure 6A, B).” (lines 404-407)\nIn a recently submitted data report to the Expedition 385 IODP volume, we are revisiting this issue more extensively (some reviewers might think obsessively) by examining declining DNA yield trends with depth and temperature for multiple sample sets. DNA yields (and parallel cell counts) drop dramatically at hot sites (Ringvent) and much more gradually at more temperate sites (all others). We do not rule out that improvements in DNA extraction could at some point yield MAGs from sediments with higher temperatures. As long as cells exist, there should be DNA present; the question is how to recover it in sufficient amounts and quality.\nThank you for your constructive and eminently actionable comments, they have contributed substantially to the revision and improvement of our study.\nReviewer #3 (Remarks to the Author):\nThe authors have sufficiently and thoughtfully addressed my comments, and I feel that the paper is now suitable for publication.\nResponse: Thank you for your constructive and specific comments, they have contributed substantially to the revision and improvement of our study."}]

Decision letter