1 Introduction

Prokaryotes, which include the two domains bacteria and archaea, play important roles in the transformation of organic matter and the cycling of biogeochemical elements in marine sediment. Bacterial communities in marine sediments typically consist of a number of ubiquitous phyla, including Proteobacteria, Chloroflexi, Planctomycetes and Acidobacteria (Inagaki et al., 2006; Nunoura et al., 2012; Liu et al., 2015; Mahmoudi et al., 2015). Based on phylogenetic surveys, among these bacteria, Proteobacteria are predominant in marine sediments, and associated with biogeochemical cycles such as the carbon, nitrogen, or sulfur cycles (McCaig et al., 1999; Lenk et al., 2011; Varon-Lopez et al., 2014). The phylum Chloroflexi was reported to be dominant in deep subsurface sediments, as it could occupy up to 70% of the bacterial 16S rRNA sequences in the Mediterranean sediment (Coolen et al., 2002). Different subgroups of Chloroflexi are suggested to have ecological functions, such as CO₂ fixation, aromatic and fatty acid oxidization, acetate production and dimethyl sulfoxide utilization (Hug et al., 2013; Wasmund et al., 2014). Archaea are commonly found in marine sediments, with Thaumarchaeota often dominating the microbial community in marine sediments (Inagakiet al., 2003, 2006; Nunoura et al., 2012; Liu et al., 2015). Thaumarchaeota include almost all aerobic ammonia-oxidizers, suggesting that they might play important roles in biogeochemical nitrogen cycling (Könneke et al., 2005; Brochier-Armanet et al., 2008). Moreover, some cases, for example, Bathyarchaeia within Crenarchaeota, can constitute up to 71% of the microorganisms in volcanic ash layers in the Okhotsk Sea, as revealed by 16S rRNA gene sequence analysis (Coolen et al., 2002; Inagaki et al., 2003).

The South China Sea (SCS) is one of the largest marginal seas around the Western Pacific Ocean, covering a total area of 3.5 million km² with an average depth of over 2000 m (Xie et al., 2003). Previous studies have focused on the diversity and the functional roles of microbial communities in various areas of the SCS (Zhang et al., 2012, 2014; Jiao et al., 2015; Liu et al., 2015; Gong et al., 2017). For instance, in the Shenhu area, microbial communities were influenced by hydrate and geochemical variable, and could be differentiated from other hydrate-related sediments around the world (Jiao et al., 2015). And in the Taixinan Basin, the dominant bacterial/archaeal groups observed were mainly involved in organic matter decomposition, sulfur-oxidation and methane generation (Gong et al., 2017). In the cold seep of the northern South China Sea, the microbial structures revealed a coupled reaction of sulphate reduction and methane oxidation (Zhang et al., 2012). All mentioned above showed that the microbes are highly diverse and play a vital role in the biogeochemical cycles in the sediments of the SCS. These results were based on an analysis of the 16S rRNA gene and functional genes. However, amplicon sequencing cannot overcome PCR in obtaining accurate qualitative and quantitative biodiversity data. In contrast, high-throughput metagenomic sequencing that obtains whole genome sequences from environmental samples can avoid PCR bias. This method can simultaneously demonstrate both the taxonomic and functional diversity of a microbial community and become a powerful tool to provide an all-inclusive picture of the functional potential of an ecosystem.

In the present study, we conducted genotyping assays via two approaches (sequencing of 16S rRNA gene amplicons and metagenomic analysis) on the Illumina HiSeq 2500 platform to analyze the microbial community structures and functional profiles in sediments from the PRMB in the northeastern SCS. Our primary objectives were i) to reveal the vertical distribution of the microbial communities and their potential environmental drivers and ii) to investigate the functional genes and microbes involved in nitrogen and sulfur biogeochemical cycles.

2 Materials and Methods 2.1 Sample Collection

Two sediment cores were obtained in the site S8 in the PRMB of the SCS during a research cruise in May 2016. The location was 19.88˚N, 115.15˚E (Fig. 1). One sediment core (S8.1, length: 40 cm) was collected by using a box corer, and the other sediment core (S8.2, length: 200 cm) was collected by using a gravity corer. Water depth, temperature, salinity and oxygen saturation of the bottom waters were measured using a conductivity-temperature-depth (CTD) system (Table 1). After collection, the 1-, 20- and 40-cm layers of sediment core S8.1 were sampled for microbial analysis, and the 1-, 50-, 100-, 150- and 200-cm layers of sediment core S8.2 were sampled for microbial analysis. It is worth mentioning that certain shallow marine sediments break away when using the gravity corer. Thus, we used sediment core S8.1 to make up in the analysis for the loss of shallow sediments in core S8.2.

2.2 Environmental Factor Analysis

Sediment pH was measured after mixing the sediments with deionized water free of CO₂ at a ratio (sediment/water) of 1:5. NH₄⁺, NO₂⁻, and NO₃⁻ were extracted from sediments with 2 mol L⁻¹ KCl as previously described (Hou et al., 2013) and quantified with a QuAAtro nutrient auto analyzer (Seal Analytical Ltd., UK). The total organic carbon (%TOC) was measured after removing inorganic carbon by digesting with HCl and using an elemental analyzer (PE2400II, UK).

2.3 DNA Extraction

Total genomic DNA was extracted from approximately 0.25 g of wet weight sediment using a PowerSoil DNA Isolation Kit (Mo Bio, USA) according to the manufacturer's protocol.

2.4 Quantitative PCR

Quantitative real-time PCR was performed to quantify the abundance of the total bacteria and archaea in the eight core samples as previously described (He et al., 2015) with minor modification. The bacterial 16S rRNA gene copies were quantified using the universal primer pair 338F/806R (Peiffer et al., 2013), and the archaeal 16S rRNA gene copies were quantified using the primer pair U519F/806R (Porat et al., 2010). A 20-µL reaction for the amplification of the 16S rRNA gene contained 10 μL of FastStart Universal SYBR Green Master (Rox) (Roche Diagnostics, Germany), 0.3 μL of each primer, 0.2 μg μL⁻¹ bovine serum albumin (BSA) and 2.0 μL of sediment DNA. The procedure employed for the PCR was as follows: an initial denaturation step at 95℃ for 10 min and 40 cycles of 95℃ for 15 s and 58℃ for 2 min. The assays were performed on an ABI PRISM® 7500 Sequence Detection System (Applied Biosystems, USA) with ABI PRISM 7500 software (version 1.3.1). The standard curves for each assay were obtained by using serial 10-fold dilutions of plasmids carrying the target gene fragments. Only standard curves with linear relationships and R² > 0.99 were used. Amplification efficiencies were 0.90 and 0.99 for bacteria and archaea, respectively. All qPCRs were performed in triplicate for each sample.

2.5 16S rRNA Gene Pyrosequencing and Data Analysis

Total DNA from the eight core samples was used as templates for amplification using barcoded primers (515F/806R for bacteria and Arch519F/Arch915R for archaea) (Coolen et al., 2004; Hori et al., 2014) targeting the 16S ribosomal RNA (rRNA) genes. Sequencing was performed on the Illumina HiSeq 2500 PE250 platform (Total Genomics Solution Institute, Shenzhen, China). The raw data were processed using the QIIME pipeline (Caporaso et al., 2010). Raw reads that had a quality score higher than 20 over a 50-bp window size and a minimum length of 50 bp were retained (Bokulich et al., 2013). The pair-end reads were joined with at least a 10-bp overlap using FLASH (Magoč and Salzberg, 2011). After quality control, UPARSE was employed to cluster all of the clean reads into operational taxonomic units (OTUs) at a 97% dissimilarity level with default parameters. The most common sequences from each OTU were selected as the representative sequences. Then, taxonomy was assigned against the SILVA 132 database via QIIME with a confidence level of 80%.

To fairly compare all the samples at the same sequencing depth, sequences were randomly reduced in number to the smallest read numbers in each sample. Then, diversity estimators (Chao1, Shannon and Simpson) and Good's coverage estimator in each library were calculated using the QIIME software package. For beta diversity, the cluster dendrogram based on Bray-Curtis dissimilarity at the OTU level was constructed with the R package vegan to describe the relationship of bacterial and archaeal communities among samples. Redundancy analysis (RDA) with Monte Carlo permutation tests (999 unrestricted permutations, P < 0.05) was conducted using the R package to reveal the correlation between community composition and environmental factors. Spearman's correlation analysis for microbial abundance, alpha diversity, percentage composition of taxa and environmental factors was performed by using SPSS statistics software (v17.0 for Windows).

2.6 Metagenomic Analysis

The metagenomic DNA isolated from 1- and 40-cm samples from S8.1 was sent to the Total Genomics Solution Institute (Shenzhen, China) for metagenome library construction and Illumina high-throughput sequencing. Paired-end (PE) raw reads were generated after sequencing, and clean data were extracted from the raw reads following the removal of adaptor fragments and low-quality reads. After quality control, clean reads were then assembled into scaffolds by using IDBA software (Peng et al., 2010). Then, open reading frames (ORFs) were predicted using MGA (version 2.10), and a BlastP search (v2.2.30) was conducted on the ORF sequences, comparing these sequences to those in the National Center for Biotechnology Information (NCBI) non-redundant protein ('nr') and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases with an E-value cut-off of 1e⁻⁵ to obtain functional and taxonomy annotation.

2.7 Nucleotide Sequence Accession Numbers

The sequence data generated in the present study have been deposited into the NCBI SRA database under the accession numbers SRP102424 (16S rRNA gene) and SRP102424 (metagenomic datasets).

3 Results 3.1 Environmental Factors at Sampling Stations

Detailed environmental factors of the sediment depth profile were shown in Table 2. In terms of inorganic nitrogen contents, higher NH₄⁺ concentrations were observed in the shallow layers (S8.1:1–40 cm and S8.2–1 cm) than in deep layers (S8.2: 50–200 cm). There was no obvious trend for NO₃⁻ and NO₂⁻ the sediment depth. Sediment pH varied from 7.62–8.37, increasing with sediment depth. The TOC content showed a clear decreasing trend with increasing sediment depth, with a range of 0.49%– 0.86%.

3.2 Quantitative PCR

The results of the real-time PCR experiment showed that the abundance of bacteria and archaea ranged from 1.29 × 10⁷ copies g⁻¹ to 6.93 × 10⁸ copies g⁻¹ (wet weight) and 2.04 × 10⁶ copies g⁻¹ to 8.20 × 10⁷ copies g⁻¹ (wet weight), respectively, in the eight core samples (Fig. 2). The archaeal gene copy numbers were equivalent to approximately 9.18%–15.79% of the bacterial gene copy numbers at corresponding depths. Meanwhile, the abundance of the bacterial and archaeal 16S rRNA genes decreased with increasing depth, decreasing by one order of magnitude at depths greater than 50 cm.

3.3 Bacterial and Archaeal 16S rDNA Analysis

Illumina-based analysis of the 16S rDNA amplicons recovered a total of 222033 high-quality bacterial sequences with an average length of 256 bp from the eight core samples (Table 3). When OTUs were grouped at a 97% similarity level, there were a total of 3813 OTUs, ranging in number from 886 to 1918, in the eight core samples (Table 3). Good's coverage estimator revealed that 96.59%–98.83% of the bacterial OTUs were present in each sediment sample, which indicated that a majority of the bacterial OTUs had been captured. The Chao1 index varied between 1083 and 2567, and the Shannon index varied between 7.29 and 8.97. These indices showed that both species diversity and richness were decreased in the deep sediments (50–200 cm) of sample S8.2.

A total of 53 bacterial phyla were identified from the eight core samples. The communities were dominated by Chloroflexi, Proteobacteria, Planctomycetes and Acidobacteria (Fig. 3). In this study, we identified Chloroflexi (31.27% of the bacterial sequences) and Proteobacteria (22.12% of the bacterial sequences) as the predominant bacterial groups. Dehalococcoidia (DEH) was the most abundant class within Chloroflexi. The relative abundance of this class increased with sediment depth, making it the most dominant class in the deep layers (50-200 cm) of sample S8.2. Within Proteobacteria, the predominant group was Deltaproteobacteria (11.10%), followed by Gammaproteobacteria (6.99%) and Alphaproteobacteria (4.02%). These two groups were dominant in the shallow sediment layers (1-40 cm) of S8.1 and S8.2, and their relative abundance decreased with sediment depth. Sequences belonging to Planctomycetes were detected at all sediment depths, and this phylum mainly dominated the shallow sediments.

A total of 206612 high-quality archaeal sequences with an average read length of 385 bp were obtained from the eight core samples, and each sample contained 14, 988 to 31842 reads (Table 3). Each archaeal library was composed of 971 to 1315 OTUs at 97% sequence similarity. High Good's coverage revealed that ≥ 96% of the archaeal OTUs were well captured in each sediment sample, indicating that sufficient sequencing depth was achieved for archaeal sequences. The Chao1 index varied between 1236 and 1893, and the Shannon index varied between 6.13 and 7.54. The Chaol and Shannon indices showed that the archaeal richness and diversity were lower than the richness and diversity of the bacterial communities.

Nine different archaeal phyla were detected across the sediment samples. The predominant phyla in the sediment samples were Crenarchaeota (38.90%) and Thaumarchaeota (21.89%), followed by Nanoarchaeota (12.87%), Asgardaeota (11.33%) and Euryarchaeota (7.91%) (Fig. 4). Within Crenarchaeota, a majority of the sequences were assigned to the class Bathyarchaeia (38.50%). The relative abundance of this class increased with sediment depth, making Bathyarchaeia the most dominant class in the deep sediment layers (50–200 cm) of sample S8.2. Nitrososphaeria was the most abundant class within Thaumarchaeota, which dominated the shallow sediment layers (1–40 cm) of S8.1 and S8.2 and even accounted for 75.63% of the surface sediment (1 cm) of S8.1. The archaeal sequences associated with Woesearchaeota had the highest relative abundances in the shallow sediment layers (1–40 cm) of samples S8.1 and S8.2.

3.4 Microbial Community Comparison and Environmental Factors Explaining Community Variations

A cluster analysis of the sediment samples showed similar trends for bacterial and archaeal communities; the shallow sediment layers (1–40 cm) of the S8.1 and S8.2 samples showed strong clustering, and the deep sediment layers (50–200 cm) of the S8.2 sample also showed strong clustering (Fig. 5). The results also revealed that sediment depth was one of the factors contributing to the variability in community structure. Notably, the bacterial community in the surface sediment of S8.2 was similar to that of the 20-cm layer of S8.1, while the archaeal community in the surface sediment of S8.2 was similar to that of the 20-cm and 40-cm sediment layers of S8.1. This similarity is because certain shallow sediment layers in S8.2 were lost when collection was performed using a gravity corer.

Correlations between the environmental parameters and the bacterial and archaeal communities were analyzed by redundancy analysis (RDA) (Fig. 6). Bacterial and archaeal communities in different sediment samples were structured by distinct environmental factors. For bacteria, TOC (r² = 0.98, P = 0.001), sediment depth (r² = 0.78, P = 0.022) and pH (r² = 0.72, P = 0.03) were the most significant factors in explaining the total variation (Fig. 6a). For archaea, TOC (r² = 0.93, P = 0.004), sediment depth (r² = 0.72, P = 0.044) and NH₄⁺ (r² = 0.68, P = 0.048) were significantly influential factors (Fig. 6b).

Spearman's analysis showed that TOC, NO₂⁻ and NH₄⁺ were positively correlated with the abundances of both bacteria and archaea (P < 0.05), while the pH and sediment depth showed the opposite trend (P < 0.05) (Fig. 7). For the bacterial communities, sediment depth and pH were negatively correlated with alpha diversity indices (P < 0.05), including Shannon, Simpson and Chao1 indices. However, TOC, NO₂⁻ and NH₄⁺ were positively correlated with the Shannon index of the bacterial community (P < 0.05) (Fig. 7). However, no environmental factor was found to have a significant correlation with alpha diversity indices of the archaeal community (Fig. 7).

3.5 Metabolic Profiles of Microbial Communities Based on Metagenomics Analysis

In this study, more than 20 Gb of sequencing data were generated for the sediment samples. After trimming, there were 10125.84 and 11710.60 Mb of sequencing data that remained for the 1-cm and 40-cm samples, respectively. By assembly and gene prediction, 656194 (808188 ORFs; average length, 377.03 bp) and 766739 (944967 ORFs; average length, 375.23 bp) scaffolds were obtained for the 1-cm and 40-cm samples of S8.1, respectively (Table 4).

3.5.1 Genes involved in nitrogen metabolism

The detected genes were involved in five nitrogen metabolism pathways: nitrate reduction, dissimilatory nitrate reduction, denitrification, nitrogen fixation, and nitrification pathways (Fig. 8a). The assimilatory nitrate reduction pathway includes the nasA, nasB, narB, nR, and nirA genes. Genes including nasA, narB and nirA were detected in the two sediment samples but had lower relative abundance. The narB and nasA gene product catalyzes the first step in this pathway, specifically the conversion of nitrate to nitrite, and the nirA gene product catalyzes the nitrite to ammonium reaction (Fig. 9). Deltaproteobacteria and Planctomycetes contributed most to narB in the studied sediments (Fig. 9a). Members that contributed to the nasA gene were Gammaproteobacteria, Actinobacteria and Verrucomicrobia in the two samples. Nitrospina of Nitrospinae and Nitrospira of Nitrospirae contributed a large proportion to nirA expression in both sediment samples (Fig. 9b). The dissimilatory nitrate reduction pathway includes the nirB, nirD, nrfA, and nrfH genes. All these genes were detected in the two sediment samples but had lower abundance (Fig. 8a). Proteobacteria, including Deltaproteobacteria, Gammaproteobacteria and Alphaproteobacteria, contributed mostly to the nirB, nirD, nrfA and nrfH genes (Fig. 9a). However, Planctomycetes, Verrucomicrobia and Actinobacteria also contributed a large proportion to nirB, nirD, nrfA, and nrfH (Fig. 9a). Biological denitrification is a key process in the marine nitrogen cycle in which NO₃⁻ is converted to gaseous products, including NO, N₂O and N₂. This pathway includes the narGHI, napAB, nirKS, norBC, and nosZ genes. In the denitrification pathway, all genes were present in the two sediment samples and were higher in abundance (Fig. 8a). Proteobacteria, including Deltaproteobacteria, Gammaproteobacteria, Alphaproteobacteria and Betaproteobacteria, were the most important contributors of all genes involved in this pathway, whereas other members, such as Planctomycetes, Chlorobi, and Bacteroidetes, were also shown to play important roles in denitrification (Fig. 9a). Importantly, archaea Thaumarchaeota contributed significantly to the presence of nirK in S8.1 and S8.2 (Fig. 9a). The nitrogen fixation pathway was not active in the sediment samples, and only the nifH gene was detected in the two samples (Fig. 8a). Methanocella and Methanosarcina were the most important contributors to this pathway (Fig. 9b). The nitrification pathway includes the microbial oxidation of ammonia to nitrite and subsequently to nitrate. Ammonia oxidation comprises two steps: the first step, the oxidation of ammonia to hydroxylamine (NH₂OH), is catalyzed by ammonia monooxygenase (Amo), and hydroxylamine oxidoreductase (Hao) then further oxidizes NH₂OH to NO (Caranto and Lancaster, 2017). The genes amoA, amoB, and amoC were mainly expressed by Thaumarchaeota, Crenarchaeota and Betaproteobacteria in two sediment samples, whereas hao was mainly expressed by Candidatus Scalindua within the Planctomycetes phylum (Fig. 9b). The nxrA and nxrB genes encode products that catalyze the conversion of nitrite to nitrate in the nitrification pathway. Chlorobi and Candidatus Acetothermia were the most important contributors of nxrA and nxrB, respectively, whereas Gammaproteobacteria also contributed a large proportion to nxrA and nxrB genes (Fig. 8a). However, the anaerobic ammonium oxidation (anammox) pathway, which seems to be of ecological importance in marine environments (Oshiki et al., 2016), was not detected in the two sediment samples. These results indicated that denitrification and nitrification pathways were active in the microbial communities of our studied sediments and played important roles in the nitrogen cycles. In addition, the relative abundances of most genes involved in nitrogen metabolism other than narB, nirB and nifK were higher in the 1-cm layer than in the 40-cm layer, implying that the nitrogen metabolism might more active in the surface layer.

3.5.2 Genes involved in sulfur metabolism

According to the KEGG database, the sulfur metabolism includes assimilatory sulfate reduction, dissimilatory sulfate reduction and oxidation and SOX (sulfur oxidation) systems Assimilatory sulfate reduction is a pathway used by a wide range of organisms to convert inorganic sulfate to sulfide, which is further incorporated into the biosynthesis of S-containing amino acids (Brunold, 1993). In this pathway, sulfate is first activated by the catalysis of sat and cysND gene products, forming adenosine 5'-phosphosulfate (APS). Then, APS is reduced to sulfite by the cysCH gene products and further reduced to sulfide by cysJI and sir gene products. All of the above genes were detected in the two sediment samples (Fig. 8b). Proteobacteria contributed the most to all these functional genes (Fig. 10a). Moreover, Candidatus Scalindua within Planctomycetes also contributed a large proportion to cys-CDHN genes, and Thaumarchaeota contributed a large proportion to sat, cysCH and sir. (Fig. 10). The dissimilatory sulfate reduction and oxidation pathways include sat, arpAB and dsrAB. Key genes (dsr and apr) potentially enabling reverse sulfite dissimilatory reduction and potentially involved in the oxidation of sulfide to sulfite were identified in both sediment samples. Deltaproteobacteria and Gammaproteobacteria of Proteobacteria contributed mostly to dsrAB, while Alphaproteobacteria, Deltaproteobacteria and Planctomycetes contributed mostly to aprAB (Fig. 10a). The Sox (sulfur oxidation) system is a well-known sulfur oxidation pathway that is catalyzed by the Sox enzyme system. The model of the Sox enzyme system comprises the four periplasmic complexes SoxXA, SoxYZ, SoxB and SoxCD, which catalyze thiosulfate oxidation (Meyer et al., 2007a). Proteobacteria, including Alphaproteobacteria, Betaproteobacteria, and Gammaproteobacteria, contributed most to the genes encoding the Sox enzyme system. Moreover, candidate division NC10 contributed significantly to the presence of soxY, whereas Acidithiobacillia and Aquificae also contributed significantly to the presence of soxA and soxB in the sediment samples (Fig. 10a).

4 Discussion 4.1 Bacterial and Archaeal Abundance in the Sediments

Bacteria were more abundant than archaea in the various marine environments (Schippers et al., 2006; Schippers et al., 2012; Liu et al., 2015), as confirmed by the results of the present study. The trends observed in this study for microbial abundance, which decreased with increasing depth, were consistent with the global distribution pattern of benthic microbial abundance (Kallmeyer et al., 2012). In this study, the bacterial and archaeal 16S rRNA gene abundances were positively correlated with NH₄⁺ (P < 0.01), NO₂⁻ and TOC (P < 0.05) but negatively correlated with pH (P < 0.01). Organic matter greatly constrains microbial abundance because this abundance highly depends on organic matter as a carbon source and main electron donor (Zhou et al., 2002; Drenovsky et al., 2004).

4.2 Shift in Microbial Community Structure with Depth and Potential Environmental Drivers

Our study revealed that the microbial community structure shifted significantly with increasing depth in the studied sediment. Bacterial communities were significantly influenced by sediment depth, TOC and pH, while archaeal communities were significantly influenced by sediment depth, TOC and NH₄⁺. For the bacterial community, Deltaproteobacteria and Gammaproteobacteria were the most abundant classes of Proteobacteria, and their relative abundance generally decreased with depth. Groups within Deltaproteobacteria are well known as sulfate-reducing bacteria (SRBs) and can incompletely oxidize products to acetate or completely oxidized products such as carbon dioxide. As the availability of fresh organic matter decreases with sediment depth, SRBs might be comparably less competitive than other metabolic guilds in the anaerobic food chain and thus decrease with depth (Jochum et al., 2017). However, predominant bacterial groups Chloroflexi were believed to have an anaerobic and heterotrophic lifestyle that depends on the recalcitrant organic matter buried in deep sediments (Oni et al., 2015), leading to a higher relative abundance in the deeper sediments. Gammaproteobacteria in this study were mainly comprised of Woeseiaceae and Halomonadaceae. The Woeseiaceae were mostly represented by the type genus Woeseia, and Halomonadaceae was mostly represented by the type genus Halomonas in this study. Woeseia showed a significantly positive relationship with TOC (r = 0.778, P < 0.05), which agreed with its chemoheterotrophic lifestyle using organic molecules as energy and carbon sources (Muβ mann et al., 2017). Halomonas preferred relatively oligotrophic environments (Keshtacher-Liebso et al., 1995), as evidenced by the significantly negative relationship with TOC in our study (r = −0.731, P < 0.05). In addition, pH was regarded as another environmental parameter that determines the bacterial community according to RDA. Krause et al. (2012) found that moderate changes in pH can cause shifts in bacterial communities. pH may affect the bacterial community through direct physiological mechanisms and may also reflect the indirect influences of other unmeasured factors (Yannarell et al., 2005; Krause et al., 2012).

The archaeal community also showed a clear stratified distribution, which was confirmed by RDA. The relative abundances of Nitrososphaeria and Woesearchaeota were enriched in the shallow sediments. The class Nitrososphaeria, which is mainly involved with aerobic oxidization of ammonium, was enriched in the surface layer and soon decreased with increasing depth. Similarly, Woesearchaeota (formerly the Deep-sea Hydrothermal Vent Euryarchaeotic Group 6, DHVEG 6) prefer to grow in relatively oxidative environments (Nunoura et al., 2012) and were indicated to play a role in organic mineralization (Lin et al., 2014). Thus, oxygen or redox state might have been the crucial factor contributing to the differences between shallow and deep layer communities. Compared to bacteria, archaea have major importance in the prokaryotic community structure in deep layer, where less-reactive organic matter is buried. It has been proposed that archaea can be united by a universal ecological ability to cope with energy stress (Gold, 1992). In the present study, Bathyarchaeia and Thermoplasmata were enriched in the deep layers and increased accordingly with depth. The phylum Bathyarchaeota, formerly named the Miscellaneous Crenarchaeotal Group, is widespread and abundant in marine subseafloor sediments and is composed of important components of anaerobic microbial communities (Kubo et al., 2012). Studies have shown that these microbes can anaerobically utilize detrital proteins (Lloyd et al., 2013), aromatic compounds (Meng et al., 2014), lignin (Yu et al., 2017) or other buried organic matter as growth substrates for biosynthesis and energy production. Microorganisms from the archaeal class Thermoplasmata were found to be linked to methanogenic activities (Iion et al., 2013; Poulsen et al., 2013). The increased abundance of these groups suggested that they might play an important role in the anaerobic degradation of organic carbon in the studied sediments.

4.3 Metabolic Profile 4.3.1 Nitrogen metabolism

The ocean's nitrogen cycle is driven by complex microbial transformations, including nitrogen fixation, assimilation, and loss. However, genes of nitrogen fixation and assimilatory nitrate reduction had lower abundance (Fig. 8a), suggesting that these two pathways might not be active in the studied sediment. Gammaproteobacteria, Deltaproteobacteria, Alphaproteobacteria and Betaproteobacteria of the Proteobacteria were found to play the most important role in the nitrogen cycle. These microbial groups contributed most to functional genes involved in nitrification, denitrification, assimilation and dissimilatory nitrate reduction. In this study, Candidatus Scalindua within the phylum Planctomycetes was the most abundant bacterial genus revealed by 16S rRNA gene sequencing. Candidatus Scalindua is often described as a marine anammox bacterial genus because the 16S rRNAs of anammox bacteria retrieved from marine ecosystems are exclusively affiliated with this bacterial group (Schmid et al., 2007; Woebken et al., 2008; Dang et al., 2013). However, according to the metagenomic analysis in our study, no genes affiliated with the anammox pathway were detected, implying that anammox might not be active in the studied sediments. In contrast, metagenomic analysis in our study showed that Candidatus Scalindua contributed most to the genes hao, nxrA, narG, nirS, nirB, nrfH and norC, which are involved in nitrification, dissimilatory nitrate reduction and denitrification. This observation is similar to those reported from metagenome analysis for 'Candidatus Scalindua' species, in which functional genes such as hao, narG and narS were discovered in these microbes (van de Vossenberg et al., 2013; Oshiki et al., 2017; Rambo et al., 2019). In summary, it was confirmed that the most abundant genus, Candidatus Scalindua, might play a more important role in nitrification, dissimilatory nitrate reduction and denitrification in nitrogen metabolism than anammox. In the nitrification pathway, Crenarchaeota contributed a large proportion of the amoA gene, and Thaumarchaeota contributed most to the amoB and amoC genes in the two sediment samples (Fig. 9a), implying that archaea are important players in the nitrification pathway in our study.

4.3.2 Sulfur metabolism

In our study, metagenomic analysis revealed that Proteobacteria contributed most to the functional genes involved in sulfur metabolism, suggesting that these microbes possibly played an important role in sulfur reduction and oxidation. This observation was in line with the previous observations that Proteobacteria dominated the sulfate-reducing prokaryote (SRP) and sulfur-oxidizing prokaryote (SOP) communities analyzed by functional gene (e.g., dsrB, aprA and soxB) sequencing (Meyer et al., 2007b; Jiang et al., 2009; Zhang et al., 2017). However, Thaumarchaeota and Candidatus Scalindua (Planctomycetes phylum, Candidatus Brocadiales order) were also found to contribute a large proportion of functional genes involved in assimilatory sulfate reduction (Fig. 10b). A previous report found that two Thaumarchaeota genomes recovered from metagenomic data in the Amazon River can encode a complete assimilatory sulfate reduction pathway (Pinto et al., 2020). Another previous report predicted Candidatus Brocadiales within Planctomycetes to be an important and unrecognized versatile player in C, N and S cycling in anoxic groundwater (Starke et al., 2017). It is noteworthy that we identified Thaumarchaeota and Candidatus Scalindua as potential sulfate reducers in the deep-sea sediments of the SCS for the first time, which expands our knowledge of the taxonomic distribution of these microbes in deep-sea systems.

5 Conclusions

In summary, the microbial communities and the potential functions of these communities in sediments from the PRMB in the SCS were analyzed. Proteobacteria, Nitrososphaeria and Woesearchaeota were dominant in the shallow sediments, while Chloroflexi and Bathyarchaeota were dominant in the deep sediments (50-200 cm). According to metagenomics analysis, we confirmed that Proteobacteria contributed most to almost all pathways in nitrogen and sulfur metabolism. Moreover, Thaumarchaeota and Candidatus Scalindua also play a vital role in some pathways of nitrogen and sulfur metabolism.

Acknowledgements

We thank Mr. Hong Qiu for analyzing the TOC content. We are also grateful to all staff on the Dongfanghong 2 for assistance with the collection of samples and geochemical data during the cruise. This work was supported by the National Natural Science Foundation of China (Nos. 41620104001 and 41806131) and the Scientific and Technological Innovation Project of the Qingdao National Laboratory for Marine Science and Technology (No. 2016ASKJ02).