Microbial community structure in hadal sediments: high similarity along trench axes and strong changes along redox gradients

Hadal trench sediments are hotspots of biogeochemical activity in the deep sea, but the biogeochemical and ecological factors that shape benthic hadal microbial communities remain unknown. Here, we sampled ten hadal sites from two trench regions with a vertical resolution of down to 1 cm. We sequenced 16S rRNA gene amplicons using universal and archaea-specific primer sets and compared the results to biogeochemical parameters. Despite bathymetric and depositional heterogeneity we found a high similarity of microbial communities within each of the two trench axes, while composition at the phylum level varied strongly with sediment depth in conjunction with the redox stratification into oxic, nitrogenous, and ferruginous zones. As a result, communities of a given sediment horizon were more similar to each other across a distance of hundreds of kilometers within each trench, than to those of adjacent horizons from the same sites separated only by centimeters. Total organic carbon content statistically only explained a small part of the variation within and between trenches, and did not explain the community differences observed between the hadal and adjacent shallower sites. Anaerobic taxa increased in abundance at the top of the ferruginous zone, seeded by organisms deposited at the sediment surface and surviving burial through the upper redox zones. While an influence of other potential factors such as geographic isolation, hydrostatic pressure, and non-steady state depositional regimes could not be discerned, redox stratification and diagenesis appear to be the main selective forces that structure community composition in hadal sediments.

In the Atacama Trench, all sediment cores were obtained using a multicorer and had clear overlying water in the core liners, indicating that the sediment surfaces were only minimally disturbed by coring. In the Kermadec Trench, undisturbed cores were obtained from sites 4 and 6 using an autonomous benthic lander and a multicorer, respectively, while other samples from site 3, 4, 5 and 7 were obtained by sub-coring a 50x50 cm boxcorer. Porewater chemistry indicated losses of the upper ≤1, ≤1, and 5-7 cm at site 3, 4, and 5, respectively, whereas at site 7 the samples showed no sign of disturbance.

Sediment core processing
The sediment cores were sectioned in a 3 °C cold room immediately upon shipboard arrival, using sterilized utensils and two different sectioning schemes with coarser and finer spatial resolution. The coarser resolution scheme (referred to as CR) followed the guidelines of the ABYSS project with triplicate cores sectioned in 0-1 cm, 1-3 cm, 3-5 cm, 5-10 cm, 10-15 cm, 15-30 cm and 30 -bottom of the core intervals, placed into sterile plastic bags (VWR), homogenized and immediately frozen to -80 °C. The higher resolution (HR) scheme was included to resolve changes in community composition along the relatively steep redox gradients of the Atacama Trench and involved one core per site sectioned in 1 cm horizons down to 10 cm, 2.5 cm down to 30 cm, followed by 5 cm slices until the bottom of the cores. Subsamples from each horizon were taken using sterile, cut off syringes and stored in PCR graded Eppendorf Tubes at -80 °C.

Biogeochemical categorization
In parallel to the molecular work of this study, a multitude of studies from different fields was conducted on either the same sites or even same samples. Analyses included prokaryotic and viral counts [1], oxygen profiles [2], and porewater chemistry, the results of which serve as metadata for this study. The oxygen profiles and porewater analyses were used to assign the sediment horizons of this study to the redox categories oxic, nitrogenous, or ferruginous as defined by Canfield and Thamdrup (2009) and exemplified in Supplement Figure 1. Horizons were classified as oxic if their mid-depth was shallower than the oxygen penetration depth [2]. Similarly, deeper horizons were classified as nitrogenous if nitrate was present until the mid-depth or deeper (Supplement Table 1

Clean-up :
The PCR products were pooled and cleaned using 1X AMPure XP beads, and amplicon lengths were checked with the DNA High Sensitivity LabChip kit (Agilent Technologies, Santa Clara, CA, USA).
Subsequently the concentration of the purified PCR products was quantified with a Qubit fluorometer (Invitrogen, Carlsbad, CA, USA).

Amplicon library preparation
From each purified PCR product pool, one hundred ng were end-repaired, A-tailed and ligated to

Demultiplexing, ASV generation, taxonomic classification, phylogeny reconstruction and filtration of negative controls
Due to the aforementioned library preparation, each of the two paired-end read files contained a mixture of forward and reverse reads. We sorted, renamed, and primer clipped these via the combination of an in-house script and the cutadapt tool. Subsequently the identification and removal of reads without an associated forward or reverse counterpart was conducted using BBmap repair.
ASVs were exported to a FASTA file (seqinr R package version 3.4-5) and a de novo alignment calculated with MAFFT using default parameters [8,9]. Subsequently, a phylogenetic tree was constructed with FastTree, using the default model in double precision mode due to the low sequence dissimilarity [10].
Negative controls were taken and treated equally to the samples along the entire laboratory process.
This includes sampling procedure (containers for the sediment samples), DNA extractions and PCRs.
Ultimately, we pooled all negative control samples together and removed the thereby identified contaminating ASVs from the entire ASV table via prevalence (threshold = 0.5) with the decontam R package version 1.2 [11].

Microbiome analyses
Biostatistical analyses and visualizations of the microbiome data were conducted with the phyloseq and ampvis2 packages in RStudio [12,13]. Prior to calculating ordinations, ANOSIM, variation partitioning, and RDA, cumulative sum scaling normalization was applied using the metagenomeSeq R package [14]. For relative abundance, absolute abundance and core-microbiome analyses we rarefied our data to even sampling depth.
Ordination plots were produced using principle components analysis of Bray Curtis dissimilarities and distinctions between groups tested using the ANOSIM implementation in the vegan package [15]. Core-microbiome analyses were conducted in the ampvis2 environment, with a relative frequency cutoff (cut_f = 50) and relative abundance cutoff (cut_a = 0.05) to account for the large variability in sediments and fine phylogenetic resolution provided by ASVs [16].
We z-scored our metadata and Hellinger transformed the individual ASV counts for subsequent variation partitioning and redundancy analysis with the vegan package [17].

Supplementary Tables
Supplement Table 1: Sampling positions, maximum water depths, oxygen penetration depths (OPD) from Glud et al., (2021) and nitrate penetration depths (NPD) of each site in the Kermadec Trench and Atacama Trench. The * indicate samples were taken from a Boxcorer. The ** shows an estimated oxygen penetration based on an empirical relation between DOU and OPD while *** imply that the value is based on only one observation, the remaining profiles did not reach anoxia (see Glud