Vertical diversity and association pattern of total, abundant and rare microbial communities in deep-sea sediments.

Microbial abundance and community composition in marine sediments have been widely explored. However, high-resolution vertical changes of benthic microbial diversity and co-occurrence patterns are poorly described. The ecological contributions of abundant and rare species in sediments also remain largely unknown. Here, by analysing microbial populations at 14 depth layers of 10 subseafloor sediment cores (water depth 1,250-3,530 m) obtained in the South China Sea, we provided the vertical profiles of microbial β-diversity and co-occurrence influenced by subcommunities of different abundance. These 134 sediment samples were clustered into four groups according to sediment depth (1-2, 6-10, 30-90 and 190-790 cm) with obvious shifts in microbial community compositions. The vertical succession of microorganisms was consistent with redox zonation and influenced by terrestrial inputs. Partitioning of vertical β-diversity showed extremely high species replacement between deep layers and the surface layer, indicating selection-induced loss of rare species and dispersal of dormant cells and spores. By contrast, for horizontal β-diversity, richness of rare species became increasingly significant in deep sediments. Accompanying this β-diversity profile were clear changes in the association pattern, with microorganisms being less connected in deeper sediment layers, probably reflecting reduced syntrophic interactions. Rare species accounted for an indispensable proportion in the co-occurrence network, and tended to form complex "small worlds." The rare subcommunity also responded differently to various environmental factors compared with the abundant subcommunity. Our findings expand current knowledge on vertical changes of marine benthic microbial diversity and their association patterns, emphasizing the potential roles of rare species.

INTRODUCTION

Marine sediments are the largest reservoir of organic carbon (OC) (Hedges & Keil, 1995) and are populated by a diverse microbial community (Parkes et al., 2000; Zinger et al., 2011). Bacteria and archaea in sediments play key roles in organic matter (OM) remineralization, nutrient cycling and energy transfer (Mason et al., 2009). Their metabolic activities profoundly impact the biogeochemical cycles of carbon, nitrogen and sulphur on a global scale (Berner, 1982; Nealson, 1997; Orcutt et al., 2011). As sediments accumulate over geological time, buried microbial life persists at decreasing growth rates due to the depletion of respiratory electron acceptors and the lesser reactivity of substrates with depth (D’Hondt et al., 2004, 2019; Jorgensen & Marshall, 2016), resulting in the assembly of the deep biosphere. It has been proposed that microbial communities in this deep biosphere represent descendants of surface communities that were buried in the past (Lever et al., 2015; Starnawski et al., 2017; Walsh et al., 2016). Selection is thus considered the predominant force driving vertical community succession in the sediment column, while other processes, including diversification, dispersal and drift, might impart minor changes mainly within the surface environment (Nemergut et al., 2013; Petro et al., 2017). In‐depth knowledge of microbial biodiversity and interaction patterns is important for understanding these biogeographical processes, which is particularly true for the large and complex microbial populations in marine sediments (Hoshino et al., 2020). Previous studies have discussed vertical benthic microbial diversity in varied depth ranges (Inagaki et al., 2003, 2015; Nunoura et al., 2016, 2018; Schauer et al., 2010; Starnawski et al., 2017; Vuillemin, Vargas, et al., 2020; Vuillemin et al., 2019; Walsh et al., 2016), and highly resolved vertical profiles of microbial communities in sediments will help to better depict the stratified community structure and continuous variations in microbial diversity. However, so far, the majority of these studies have often focused on specific types of sedimentary environment (Nunoura et al., 2018; Yang et al., 2020) or comparison of microbial community over large spatial or temporal scales (Böer et al., 2009; Hoshino et al., 2020; Hoshino & Inagaki, 2019). Parallel comparison of site‐to‐site microbial diversity at different depths in a less varied sediment region is scarce.

β‐Diversity of microbial communities, a long‐standing interest in the field of marine microbial ecology, indicates how microbial assemblages vary in space and along environmental gradients. It also provides insights into the mechanism driving biodiversity changes and their consequences for multiple ecosystem functions (Mori et al., 2018). β‐Diversity reflects two different phenomena: species replacement (also called turnover) and richness difference or nestedness (Baselga, 2010). By partitioning β‐diversity into replacement and richness components, we are able to interpret the underlying ecosystem processes. In marine sediments, species replacement implies the simultaneous gain and loss of species due to processes such as selection during burial or along environmental gradients, genetic diversification and historical events (Macalady et al., 2013; Martiny et al., 2006; Nemergut et al., 2013), while richness difference may reflect shifts in the relative abundance of taxa that are always present, typified by the persistence of certain surface communities in deep sediments (Legendre, 2014). Many studies have shown that although both components contribute to total β‐diversity, their relative contributions may vary across different biological groups (Baselga, 2010; Hu et al., 2019; Si et al., 2015), emphasizing the need for β‐diversity partitioning in order to discern the antithetic biological implications underlying overall community changes. Partitioning of β‐diversity gives more insights into the factors driving spatial variability in biotic communities compared with a consideration of total β‐diversity alone (Soininen et al., 2018). For example, β‐partitioning patterns of microbes in upper lake sediment clearly distinguished taxonomic changes caused by either microbial activity or those related to sediment burial (Wurzbacher et al., 2017). A nested community selected from the surface has also been proposed for the deep biosphere in lacustrine sediment based on partitioning of β‐diversity (Vuillemin et al., 2018). However, no studies have examined how replacement and richness components of microbial β‐diversity change with depth in marine sediments and their biogeographical implications.

Patterns of microbial dynamics are associated with interspecies interactions (Gorter et al., 2020). Thus, knowledge of alterations in microbial interactions underlying community changes will help to create a better understanding of microbial ecology. Co‐occurrence networks provide system‐level insights into the complex microbial interactions (Lv et al., 2019), although there is still debate over whether this method provides true and valid interactions (Liu et al., 2019). Reconstruction of microbial networks can advance our knowledge of the complex behaviours in microbial communities and predict the effects of perturbations on community dynamics (Röttjers & Faust, 2018). Previous studies have used network analysis to explore interactions between marine benthic microbes and microbe–environment relationships (Buongiorno et al., 2019; Wang et al., 2020; Zhang, Li, et al., 2020), but most of these studies focused only on shallow sediments and none of them revealed the vertical shift of microbial networks in the subseafloor. In addition, network analysis can be a powerful tool for inferring the keystone taxa in a microbial community which can drive community composition and function irrespective of their abundance (Banerjee et al., 2018). While keystone taxa underline the importance of numerically inconspicuous taxa in microbiome functioning (Banerjee et al., 2018), they are seldom discussed for marine sediments.

Microorganisms often present a skewed abundance distribution in a local community, with a few relatively abundant species coexisting alongside a large number of rare species referred to as the “rare biosphere” (Logares et al., 2014; Lynch & Neufeld, 2015; Pedrós‐Alió, 2012; Sogin et al., 2006). Abundant and rare taxa play distinct roles in an ecosystem. The abundant taxa are important for the fluxes of OM and biomass production due to their huge quantities (Pedrós‐Alió, 2012), whereas recent studies increasingly emphasize that the rare biosphere not only serves as reservoirs of genetic and functional diversity, but also has substantial functional roles (Lynch & Neufeld, 2015; Shade & Gilbert, 2015; Shade et al., 2014). Moreover, the rare subcommunities can serve as components of microbial seed banks, which may maintain a dormant state and recruit to higher density upon changes in environmental conditions, and the dispersal of microbial seed banks appears to primarily shape subsurface community compositions across oceans (Lennon & Jones, 2011). Rare microbes have also been reported to achieve essential functions in nutrient cycling, and may enhance the functionality of the abundant microbes (Bodelier et al., 2013; Pedrós‐Alió, 2007; Philippot et al., 2013), suggesting that there could be comprehensive interactions between abundant and rare species (Jousset et al., 2017). Therefore, rare microbes could be overlooked keystone species regulating the functions of diverse ecosystems. The distribution patterns of rare taxa have been studied in different environments (Jiao & Lu, 2020; Mo et al., 2018; Wu et al., 2017), with the results indicating that microbial community assembly in the rare biosphere can be driven by distinct processes and dependent on distinct environmental variables as compared with the abundant taxa. Nevertheless, the distribution of rare microbial taxa in deep marine sediments remains largely unexplored. Additionally, the moderate subcommunity has been increasingly discussed in recent studies (Bay et al., 2020; Sinkko et al., 2019; Xue et al., 2018; Zhang, Hou, et al., 2020), which reflects an intermediate state between the abundant and rare counterparts in several important aspects, including diversity and ecological processes that structure microbial community compositions. Whether the abundant, moderate and rare subcommunities in sediments show similar vertical changes or whether they respond differentially to environmental heterogeneity is yet to be answered.

The South China Sea (SCS) is one of the largest marginal seas in the world. The Shenhu Area, located on the northern slope between the Xisha Trough and Dongsha Islands of the SCS, is an important region for methane hydrate exploration (Wu et al., 2011). Organic‐rich sediments have accumulated since the late Mesozoic. The thickness of sediment is 1,000–7,000 m with 0.46%–1.9% OM (Wang et al., 2000). A previous study revealed that hydrate‐free and hydrate‐containing surface sediments in the Shenhu Area possess distinct microbial community compositions (Jiao et al., 2015). However, a comprehensive study of vertical microbial ecology in sediments near the Shenhu Area is lacking. Here, we collected 10 sediment cores ranging from 6 to 8 m and provided a highly resolved vertical profile of microbial abundance and distribution in this area, mainly aiming to evaluate (i) how microbial β‐diversity and its two distinct components change along sediment depth, (ii) what are the vertical patterns of microbial interaction networks, and (iii) the importance of the abundant and rare subcommunities regarding microbial diversity and interactions in marine sediments.

MATERIAL AND METHODS

Sample collection and environmental characterization

Ten sediment gravity cores were collected on the northern slope of the SCS (18.44–19.26°N, 114.25–115.30°E), ranging from 6.64 to 8.36 m. Water depths at these sites were 1,250–3,530 m. The geographical information of the 10 sampling sites is summarized in Figure 1a and Table S1. The sediment samples were collected by a gravity corer during cruise HYSH201901 organized by Guangzhou Marine Geological Survey from February 26 to March 6, 2019. The cores were sliced with a stainless‐steel cutter for subsampling. To reduce contamination, the outer 1 cm of sediment was sliced off with a red‐hot, sterile spatula. An unused, sterile spatula was used to carefully sample the uncontaminated centre of the remaining core sample at each depth. The sediments were transferred to sterilized plastic tubes and stored at −80°C before OM measurement and DNA extraction. We obtained 134 samples in total at 1, 2, 6, 10, 30, 50, 90, 190, 290, 390, 490, 590, 690 and 790 cm of these cores. Sediment porewater was extracted using a Rhizon sampler attached to a vacuum tube under N2 atmosphere (Seeberg‐Elverfeldt et al., 2005; Yao et al., 2014). After collection, porewater samples were further split into acidified or unacidified subsamples, and stored at 4°C for subsequent analysis of dissolved inorganic carbon (DIC), nutrients and other solutes (Zhao, Yao, et al., 2017). DIC was measured using the nondispersive infrared method with a precision of 0.1% and a detection limit of 0.2 µM (AS‐C3, Apollo SciTech) (Cai & Wang, 1998). Dissolved inorganic nutrients (,,, and ) in porewaters were determined with the standard colorimetric methods using a continuous flow analyser (Auto‐Analyzer 3; Seal Analytical) (Catalano, 1987; Worsfold et al., 2016). Major cations (Ca2+ and Mg2+) and dissolved Mn (assumed to be present as Mn2+) were determined using inductively coupled plasma optical emission spectrometry (ICP‐OES; iCAP6300, Thermo Fisher Scientific) (Zhao, Yao, et al., 2017). Major anions ( and Cl−) in porewater were analysed using ion chromatography (ICS‐3000, Dionex) on 1:800 diluted aliquots in Milli‐Q water (Yao et al., 2014). Total OC (TOC), total nitrogen (TN) and stable carbon isotope composition (δ13CTOC) in sediments were analysed after removing carbonates with HCl (6 m) on an elemental analyser (vario MICRO cube, Elementar) interfaced to a continuous flow isotope ratio mass spectrometer (Isoprime100, Isoprime Co.) as described by Yao et al. (2014).

FIGURE 1

Site locations, environmental parameters, microbial abundance and communities in SCS sediments. (a) Sampling map of the northern slope of the SCS. (b) Depth profiles of key chemical characteristics. (c) Vertical changes of bacterial and archaeal 16S rRNA gene abundance. (d) Microbial community compositions at the phylum level and at class level of Proteobacteria. (e) db‐RDA illustrating the relationship between the microbial community at the OTU level and significantly influential environmental variables

DNA extraction, sequencing and reads processing

Genomic DNA was extracted from 0.25 g of sediment (wet weight) using the Power Soil DNA Isolation Kit (MoBio Laboratories) and a FastPrep‐24 cell disrupter (MP Biomedicals) according to the manufacturers’ instructions. The quality and quantity of the extracted DNA were measured by a Nanodrop ND‐2000 spectrophotometer (Thermo Fisher Scientific). The primer set of 515FmodF (5′‐GTGYCAGCMGCCGCGGTAA‐3′) and 806RmodR (5′‐GGACTACNVGGGTWTCTAAT‐3′) was used to amplify the V4 hypervariable region of the 16S rRNA gene from bacteria and archaea (Walters et al., 2016). These primers were modified according to Parada et al. (2016) and Apprill et al. (2015), which improves detection of Thaumarchaeota and clade SAR11 without degrading performance on taxa already amplified effectively by the original primer set. Polymerase chain reactions (PCRs) were conducted as follows: 3 min for denaturation at 95°C; 29 cycles of 30 s at 95°C, 30 s for annealing at 55°C, 45 s for elongation at 72°C, and 10 min for a final extension at 72°C. PCRs were conducted in triplicate 20 µl mixture comprising 0.4 µl of FastPfu Polymerase, 0.8 µl of each primer (5 μm), 2 µl of 2.5 mm deoxyribonucleoside triphosphates (dNTPs), 4 µl of 5× FastPfu Buffer, and 10 ng of template DNA. The PCR products were purified from 2% agarose gels using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences) and further quantified using QuantiFluor‐ST (Promega) following the manufacturers’ protocol. Purified amplicons were pooled in equimolar ratio and paired‐end sequenced (2 × 300) on an Illumina Miseq PE300 platform (Illumina) at Majorbio Bio‐Pharm Technology.

Raw reads were quality‐filtered by trimmomatic (Bolger et al., 2014) and merged by Fast Length Adjustment of Short reads (flash 1.2.11, https://ccb.jhu.edu/software/FLASH/) with the following criteria: (i) the reads were truncated at any site receiving an average quality score <20 over a 50‐bp sliding window; (ii) sequences whose overlap was greater than 10 bp were merged if there were no more than two nucleotide mismatches; (iii) sequences of each sample were separated by barcoding (exact matches) and primers (allowing two nucleotide mismatches); and (iv) reads containing ambiguous bases were deleted. Operational taxonomic units (OTUs) were clustered with a 97% similarity cutoff using uparse (Edgar, 2013) (version 7.0.1090, http://drive5.com/uparse/). Singleton OTUs that may represent sequencing errors were removed before downstream analyses. The taxonomy of each 16S rRNA gene sequence was analysed via the Ribosomal Database Project (RDP) classifier algorithm (version 2.11, http://rdp.cme.msu.edu/) (Wang et al., 2007) against the silva database (Release132, http://www.arb‐silva.de). To equalize sequencing depth, each sample was rarefied to 33,686 reads (the lowest sequence number across all samples) for further analysis. The raw data were submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database (accession no.: PRJNA606622).

Quantitative PCR

Quantitative PCR was performed to quantify the abundance of total bacteria and archaea using primer sets specific for the bacterial 16S rRNA gene (Eub338F/Eub518R) (Yin et al., 2013) and archaeal 16S rRNA gene (967F/1060R), respectively (Cadillo‐Quiroz et al., 2006). A 20‐µl mixture contained 10 µl of SYBR Premix ExTaq II (2×), 0.4 µl of ROX Reference Dye II (50×) (TaKaRa), 0.2 µl of primers for each gene (10 µm), and 2 µl of template. The thermal cycling steps for the bacterial 16S rRNA gene consisted of an initial denaturation at 95°C for 5 min, 40 cycles of 95°C for 30 s, 53°C for 1 min and 72°C for 15 s, and a final extension at 72°C for 10 min. The thermal cycling steps for the archaeal 16S rRNA gene consisted of an initial denaturation at 95°C for 30 s, 40 cycles of 95°C for 5 s, 50°C for 30 s and 72°C for 30 s, and a final extension at 72°C for 5 min. All assays were conducted in triplicate, with negative controls, using an ABI 7500 Real‐Time PCR System (Applied Biosystems). Standard curves were constructed by PCR amplifying a 10‐fold serial dilution of plasmids containing target gene fragments. The amplification curves showed good linear relationships (R 2 > .999) and the amplification efficiencies were 0.90 and 0.95 for the bacterial and archaeal 16S rRNA gene, respectively.

Data analysis and statistics

The alpha diversity indices, including Good's coverage, Chao 1, Shannon, Shannon‐even, and phylogenetic diversity (PD), were calculated using mothur (version 1.30.1) (Schloss et al., 2009). Analyses of similarities (ANOSIM) was performed with the “vegan” package in R. A linear discriminant analysis effect size (LEfSe) was applied at the genus level (all‐against‐all, LDA > 3.5; http://huttenhower.sph.harvard.edu/galaxy/) to identify the discriminant bacterial and archaeal clades at different depths.

We defined “abundant” OTUs as those having relative abundances above 0.1% of the total sequences, “moderate” OTUs as those having relative abundances between 0.01% and 0.1%, and “rare” OTUs as those having relative abundances below 0.01% (singletons, which may be caused by sequencing errors, were discarded) (Jiao & Lu, 2020). Nonmetric multidimensional scaling (NMDS) was performed based on Bray–Curtis distances with the “vegan” package in r. Tests of homogeneity of dispersions (“betadisper” in the “vegan” package) followed by ANOVAs were performed to compare the multivariate dispersions. Differences of within‐group dissimilarity were determined using the Kruskal–Wallis test followed by a post‐hoc analysis in originpro 2020 (OriginLab). Total β‐diversity and its partitioning into replacement and richness components (β repl and β rich) were calculated based on the Jaccard dissimilarity via the BAT package in r (Cardoso et al., 2015). We defined vertical β‐diversity as β‐diversity between the sum of OTUs in the surface layer (1 cm) and each deeper layer, and horizontal β‐diversity as β‐diversity calculated among samples from the same depth layer.

Co‐occurrence networks were constructed for all samples and four depth groups, respectively, using the “igraph” and “Hmisc” libraries in r. To reduce complexity, only OTUs occurring in more than 25% of samples were retained. The pairwise Spearman's correlations between OTUs were calculated, with a correlation coefficient > |.7| and a p value <.01 (Benjamini and Hochberg‐adjusted) being considered as a valid relationship. The network‐level (mean node degree, clustering coefficient, average path length, modularity, density, diameter, betweenness centralization and degree centralization) and node‐level (degree, transitivity, betweenness centrality and closeness centrality) topological features of a network were calculated. gephi (version 0.8.2 beta, http://gephi.org) was used for network visualization. Keystone species in each co‐occurrence network were defined as those with a high degree of ranking in the top 20% among all nodes, and low betweenness centrality, ranking in the bottom 20% among nodes with high degree (Berry & Widder, 2014).

We implemented the Mantel test using the “vegan” package in R to identify significant environmental factors correlated with the variations of microbial communities. The relative importance of factors in explaining community variation was investigated by variation partitioning based on distance‐based redundancy analysis (db‐RDA) of Hellinger‐transformed OTU abundance data using the “vegan” package in R.

RESULTS

General environmental characterization

Detailed sediment and porewater environmental parameters are shown in Figure 1b. In general, DIC, and showed a similar tendency to increase with depth, especially below 90 cm. and Mn2+ were also higher in deeper sediments than in surface layers (p < .001 in Mann–Whitney U tests). Conversely, higher NO3 − and NO2 − concentrations were detected in the top 10 cm (p < .001 in Mann–Whitney U tests). Ca2+ and exhibited a decreasing trend with depth. TOC and TN contents were constant with depth, but abrupt increases of TOC and C/N were observed at 50, 190 and 490 cm in several sites. δ13C was lower at 50 cm in most sites, in concert with the increases of TOC and C/N at this depth, indicating that there might be distinct sources of OC input compared with other depths.

General microbial abundance and α‐diversity in SCS sediment

Quantitative PCR showed that the 16S rRNA gene abundances of both bacteria and archaea were highest in surface sediments (1–2 cm; 108–109 copies g–1 for bacteria, 107–108 copies g–1 for archaea), and then declined gradually with depth by 1.5–3 orders of magnitude (Figure 1c). Differences were also observed among the studied sites but were less distinguishable than those among depths (ANOSIM global R: .0831, p = .001). Across all samples, the abundance of the bacterial 16S rRNA gene was higher (0.2–2.0 orders of magnitude) than that of archaea.

A total of 4,513,924 reads were obtained after quality control, discard of singletons and rarefaction, and were clustered into 29,520 OTUs at a 97% similarity level (25,896 bacterial OTUs and 3,539 archaeal OTUs). We further divided the OTUs into abundant, moderate and rare OTUs depending on their relative abundance. In the whole data set, 147 OTUs and 878 OTUs, representing 57.0% and 24.1% of all sequences, belonged to abundant and moderate OTUs, respectively; 28,495 OTUs contributing 18.9% of all sequences were classified as rare OTUs.

Good's coverage values ranged between 95.8% and 99.9% across all samples, indicating that sequences generated from these samples could represent most of the microbial community in the studied sites. All α‐diversity indices, including Chao 1, Shannon, Shannon evenness and PD, fluctuated with a similar trend with depth (Figure S1). Chao 1 and PD indices exhibited an evident median increase in the top 10 cm and a subsequent decrease to the deeper layers (290–490 cm). Shannon and Shannon evennness indices were slightly lower at 50–490 cm compared with those in the surface and deepest sediments.

Vertical and horizontal β‐diversity patterns

NMDS analysis based on the Bray–Curtis distance revealed that the sediment samples clustered more strongly by depth than by site. In general, samples were divided into four groups (Figure 2a) at the OTU level: 1–2, 6–10, 30–90 and 190–790 cm, for the whole data set (p = .001, R 2 = .36 in PERMANOVA) as well as the abundant (p = .001, R 2 = .38 in PERMANOVA), moderate (p = .001, R 2 = .40 in PERMANOVA), and rare subcommunities (p =.001, R 2 = .22 in PERMANOVA). However, the dispersion of these groups also differed significantly (“betadisper,” ANOVA, F = 16.0, p < .001). Rare OTUs best resembled the segregation of samples of the entire community. In general, the dissimilarity between samples shown by the rare subcommunity was highest compared with the abundant and moderate subcommunities at all depths, with least variation within each group (Figure 2b). Comparison of community variance at each depth group showed depth‐related horizontal heterogeneity (Figure 2b). For the entire community, within‐depth distances increased from 1–2, 6–10 to 30–90 cm, and then decreased at 190–790 cm (p < .001, χ2 = 444.1 in Kruskal–Wallis test). Similar patterns were observed for abundant (p < .001, χ2 = 274.4 in Kruskal–Wallis test) and moderate subcommunities (p < .001, χ2 = 571.9 in Kruskal–Wallis test), whereas the within‐depth distances increased with depth and were highest at 190–790 cm for rare OTUs (p < .001, χ2 = 761.5 in Kruskal–Wallis test).

FIGURE 2

Comparison of microbial communities at the OTU level in SCS sediments. (a) Nonmetric multidimensional scaling ordination of total, abundant, moderate and rare subcommunities based on the Bray–Curtis dissimilarity. (b) Comparison of Bray–Curtis dissimilarity of microbial OTUs within different depth groups for each subcommunity

To explore the extent of community turnover with depth, we compared the deeper layers with the surface layer (1 cm). The β‐diversity based on the entire community increased with depth and approached 1 (i.e., no shared taxa) in the deepest layer (790 cm) (Figure 3a). Upon partitioning the overall β‐diversity into taxonomic richness (β rich) and replacement effects (β repl), taxonomic replacement increased in the top 90 cm and was consistently high at 190–790 cm (Figure 3a). In contrast, the effect of richness was high at 6–30 cm and then decreased with depth, accounting for only minor proportions at 190–790 cm. For abundant and moderate subcommunities, replacement between surface and deeper layers significantly decreased compared with that of the entire community, while the increased and decreased trends for replacement and richness effects, respectively, were similar to the entire community. In addition, vertical changes in replacement and richness effects for the rare subcommunity exactly reflected that of the entire community.

FIGURE 3

β‐Diversity and β‐diversity partitioning patterns for the entire community and subcommunities. (a) Vertical β‐diversity between the surface (1 cm) and each deeper layer. (b) Horizontal β‐diversity among samples from the same depth layer

We further analysed horizontal β‐diversity within samples from the same sediment depth, and the results showed that total horizontal β‐diversity increased in the top 30 cm and remained at a similar level in deeper layers (Figure 3b). Replacement effects were higher than richness effects in each layer above 190 cm, and showed a decreasing trend with sediment depth. Richness effects gradually increased from the surface to 190 cm (Figure 3b). Although there was no significant change of total horizontal β‐diversity at 290–790 cm, large variations of replacement and richness effects were observed between different sites. The patterns exhibited by the entire community exactly matched those of the rare subcommunity, with slightly higher horizontal distance among rare OTUs. The horizontal distances in abundant and moderate subcommunities were lower compared with that of the rare subcommunity, and fluctuated in different layers.

Microbial community compositions at different depths

In total, 80 phyla were observed in the SCS sediment samples. The most dominant phylum was Proteobacteria accounting for 47.25% of all sequences, followed by Chloroflexi (25.12%), Planctomycetes (8.58%) and Actinobacteria (6.62%) (Figure 1d). Dominant archaeal phyla included Thaumarchaeota (6.17%), Crenarchaeota (1.99%), Nanoarchaeaeota (1.62%), Hadesarchaeaeota (0.95%) and Euryarchaeota (0.83%). The relative abundance of different microbial groups changed significantly with depth (Figure 1d). Using LEfSe analysis, the surface sediments (1–2 cm) were found to be characterized by a high abundance of Nitrososphaeria, Dadabacteriia, Nitrospinia (within the phylum Nitrospirae), and several other groups (Figure S2). Lokiarchaeia, within the phylum Asgardaeota, was a distinct archaeal group enriched at 30–90 cm, while Hadesarchaeaeota was enriched in deeper layers (190–790 cm). Bacterial groups which were significantly more abundant at 190–790 cm contained several orders from Dehalococcoidia and Aminicenantia, as well as other groups within the Chloroflexi and Deltaproteobacteria.

The abundant OTUs were primarily affiliated with Proteobacteria (38.10%), Chloroflexi (22.45%) and Thaumarchaeota (8.84%), whereas the proportions of Proteobacteria was 23.01% and 17.09% in moderate and rare OTUs, respectively. Chloroflexi (24.06%) and Planctomycetes (10.93%) were two other primary taxa of moderate OTUs. Many more taxonomic groups were identified in rare OTUs, with Proteobacteria (17.09%), Planctomycetes (11.61%) and Chloroflexi (9.77%) being the major groups, as well as a high proportion of Nanoarchaeaeota (9.08%) and Patescibacteria (7.60%). The rare OTUs also contained considerably more unclassified taxa (7.45%) compared with those in abundant and moderate OTUs (0.68% and 2.05%, respectively). The proportions of subcommunities were similar in different depth groups: abundant OTUs, 0.46%–0.59%; moderate OTUs, 2.44%–3.27%; rare OTUs, 96.14%–96.98%.

Microbial co‐occurrence networks in SCS sediment

A co‐occurrence network was built based on correlation relationships across all samples from different depths (Figure 4a). The resulting network consisted of 1,468 nodes linked by 158,298 edges (Table S2). The entire network was clearly parsed into five major modules, of which modules 1, 2 and 3, respectively, accounted for 42.9%, 38.8% and 13.2% of the whole network. Most of the OTUs from module 1, the largest module, and module 3 had higher relative abundances in surface and shallow sediments, whereas the majority of the OTUs from module 2 had higher relative abundances in deep sediments (Figure 4b). The OTUs from Proteobacteria, Chloroflexi and Planctomycetes displayed wide correlations with others (Figure 4a). A much higher number of strong positive correlations (92.4%) were observed than negative ones (7.6%), with negative correlations mainly distributed between OTUs in module 2.

FIGURE 4

Co‐occurrence networks of the microbial community in all sediment samples and its properties. (a) Co‐occurrence networks with OTUs coloured by modularity, taxonomy and subcommunity, successively. Each connection shown has a correlation coefficient >|.7| and a p value <.01. The size of each node is proportional to the number of connections. (b) Numbers of sequences at a certain depth for each OTU in different modules. (c) Comparison of node‐level topological features among subcommunities. A, abundant; M, moderate; R, rare

By partitioning nodes into subcommunities, rare OTUs accounted for nearly the same proportions as moderate OTUs in the network (Figure 4a). We further compared node‐level topological features of abundant, moderate and rare subcommunities (Figure 4c). Abundant OTUs showed the highest degree level among the three subcommunities (p < .05, χ2 = 10.6 in Kruskal–Wallis test). Clustering coefficients increased significantly from abundant and moderate OTUs to rare OTUs (p < .001, χ2 = 186.5 in Kruskal–Wallis test), while the opposite pattern was observed for closeness (p < .001, χ2 = 57.4 in Kruskal–Wallis test) and betweenness centrality (p < .001, χ2 = 205.2 in Kruskal–Wallis test). No significant differences in eigenvector centrality values among the three subcommunities were found (data not shown).

We identified 59 keystone OTUs in total with high degree (>477) and low betweenness centrality (<311), including 41 rare OTUs, 16 moderate OTUs, and only two abundant OTUs, which may play key roles in maintaining the structure and function of microbial communities. The majority of keystone OTUs belong to Proteobacteria (25), Planctomycetes (13) and Acidobacteria (9). Other keystone OTUs belong to Bacteroidetes, Chloroflexi, Nitrospirae, Acidobacteria, Patescibacteria and Dadabacteria (Table S3).

Vertical changes of co‐occurrence networks

To study the shift of microbial co‐occurrence patterns along depth, we built separate co‐occurrence networks for samples from 1–2, 6–10, 30–90 and 190–790 cm (Figure 5). The complexity of the microbial network at 190–790 cm decreased significantly (Figure 5; Table S2). Both the number of OTUs involved in the network and the number of correlations decreased, as did the diameter, average degree and average path length. The power law distributions of networks in the shallow sediments (1–2 and 6–10 cm) conformed to scale‐free network topology, with a power‐law exponent γ approaching 2 (1.840 and 1.955, respectively), while networks of deep sediments (30–90 and 190–790 cm) less resembled the scale‐free networks (γ = 1.381 and 1.316, respectively) (Figure S3). However, although there were fewer nodes in the network at 30–90 cm compared with 6–10 cm, the edge number increased drastically at 30–90 cm compared with other depths. Nodes with a high degree in the 30–90 cm network were mainly Chloroflexi, with 4,474 edges in total connected to this phylum. Rare OTUs accounted for the majority of nodes (>60%) in the networks of the 1–2, 6–10 and 30–90 cm networks, and decreased to 28% in the 190–790 cm network (Figure S4), whereas the proportions of rare OTUs as keystone nodes were higher in layers other than 1–2 cm. Keystone OTUs identified at 1–2, 6–10 and 30–90 cm mainly belonged to Proteobacteria, Chloroflexi, and Planctomycetes, while keystone OTUs in the 190–790 cm layer were primarily Bacteroidetes and Firmicutes.

FIGURE 5

Co‐occurrence networks of the microbial community in each depth group. Each connection shown has a correlation coefficient > |.7| and a p value <.01. The size of each node is proportional to the number of connections. The OTUs were coloured by subcommunity and taxonomy, respectively

The influence of environmental factor on community variation

Environmental factors with significant influence on sediment communities as revealed by db‐RDA (p < .05) are shown in Figure 1e. DIC, NH4 +, PO4 3− and SiO4 2− significantly correlated with microbial community structure at 190–790 cm, while and were correlated with samples in the top 10 cm. Samples from 30–90 cm appeared to be separated from other samples by TOC. The Mantel test results indicated that changes of total, abundant, moderate and rare communities were correlated with DIC, , , , , , , Cl−, TN and C/N (Table S4). In addition, DIC, , and TN showed stronger correlations with rare OTUs than abundant OTUs, while and conversely showed stronger correlations with abundant OTUs than rare OTUs (Table S4). The α‐diversity of the entire community was significantly influenced by DIC, , , , , Ca2+, TOC, TN, C/N, and δ13C.

DISCUSSION

Knowledge of microbial community composition, diversity and association patterns is important for a better understanding of how environmental variables and biogeographical processes influence subsurface biospheres in marine sediments. We found that microbes in deeper sediment layers exhibited a distinct β‐diversity partitioning and interaction pattern compared with those in shallower sediments. Selection was the major process determining the community assembly in sediments, while other processes, such as diversification and dispersal from terrestrial inputs, may exert impacts on the rare subcommunity. In addition, we provided evidence that the rare subcommunity played a non‐negligible role in shaping microbial diversity and interactions in marine sediments.

Microbial communities shaped by geochemical zonation and palaeoenvironmental factors in the SCS sediments

SCS is one of the largest marginal seas in the world and is abundant in biotic resources. The microbial cell abundance in surface sediment of the northern slope of the SCS has been reported to vary between 106 and 109 cells g–1 sediment (Cui et al., 2019; Jiang et al., 2007; Jiao et al., 2015). The bacterial and archaeal 16S rRNA gene abundance in surface sediments of this study was 107–109 and 106–108 copies g–1, respectively, which was one to two orders of magnitude lower compared with those in sites located closer to the Pearl River Estuary (Dang et al., 2013) with high OM from terrestrial inputs. Due to fierce selection processes in the subseafloor sediment, bacterial and archaeal abundance decreased logarithmically with depth. The α‐diversity also showed a decreasing trend from 10 to 690 cm, in accordance with previously reported distribution patterns of microbial abundance and diversity (D’Hondt et al., 2019; Kallmeyer et al., 2012).

Successive community shifts were observed with increasing sediment depth, consistent with redox zones characterized by aerobic respiration, denitrification, reduction of oxidized manganese and sulphate (Figure 1b). was rapidly depleted by abundant ammonia‐oxidizing Thaumarchaeota in the top 10 cm (Vuillemin et al., 2019), and re‐accumulated in deeper layers. The ammonia‐oxidizing process led to high concentrations of and in shallow sediments, which might support denitrification in a wide diversity of bacteria (especially Proteobacteria) (Devol, 2015). The increased Mn2+ at 30 cm indicated the activity of Mn‐oxidizing communities, possibly accompanied by Fe oxidizers. Concentrations of started to decrease below 190 cm, due to increased sulphate‐reducing bacteria including Desulfuromonadales, Desulfovibrionales, Desulfobacterales and Desulfarculales in Deltaproteobacteria. However, these sulphate‐reducing bacteria only accounted for a small proportion in Deltaproteobacteria, outcompeted by fermentative members of Chloroflexi (Dehalococcoides), Atribacteria and Aerophobetes (Kerrigan et al., 2020; Vuillemin, Vargas, et al., 2020; Wang & Li, 2016). Fermentation of organic acid was also achieved by several archaeal groups enriched in deep sediments, including Lokiarchaeota (Orsi et al., 2020) and Bathyarchaeota (He et al., 2016). The oxidation of OM through sulphate reduction and fermentation was considered the main mechanism contributing to DIC production in deep sediments (Valdemarsen & Kristensen, 2010). Methanotrophic archaea (ANME) were not found in these samples, indicating that the sulphate–methane transition zone (SMTZ) and methanogenic zone in SCS sediments were much deeper (Graw et al., 2018).

Although the vertical changes of microbial communities corresponded with the succession of predominant terminal electron‐accepting pathways with sediment depth, terrestrial inputs also play a role in shaping microbial communities, as evidenced by TOC, C/N and δ13C. Surrounded by the Asian continent and Taiwan to the north (Figure 1a), 80% of the total SCS surface sediments were composed of river‐borne terrigenous sediments, especially in the northern slope (Liu et al., 2008). The sedimentation rate in the study area reached 0.77–8.12 cm ka–1 as reported by Schönfeld and Kudrass (1993) and 2–12 cm ka–1 by Shao et al. (2001), forming organic‐rich sediment with TOC remaining at a relatively high level (around 1% except for 50 cm) with depth. There was a sudden increase in TOC and change in δ13C percentage at 50 cm in most sites, reflecting terrestrial inputs affected by historical climate oscillations such as the glacial and interglacial cycles (Lai et al., 2020; Zhao, Liu, et al., 2017). Distinct TOC, C/N and δ13C profiles of SCS005 were observed at 50–490 cm and this may explain the peculiar microbial community compositions in deep sediments of SCS005 featured by large proportions of Actinobacteria, which were considered to be dispersed from different sources of terrestrial input compared with other sites.

Therefore, the microbial community in SCS sediment has resulted from a combination of palaeoenvironmental factors and substrate depletion with depth (Figure 1e). However, few studies have differentiated abundant and rare taxa from the total microbial community when assessing its response to environmental factors. We found that although abundant, moderate and rare subcommunities were all significantly correlated with the measured factors (Table S4), the correlation coefficients varied between certain factors and the subcommunities. These results indicated that the abundant, moderate and rare subcommunities may be differentially influenced by environmental factors.

β‐Diversity patterns indicate distinct processes shaping abundant and rare subcommunities in sediments

Microbial communities in sediments are considered to derive from selection from within the surface communities, as well as dispersal of microorganisms in seawater and the terrestrial environment. To explore the detailed succession of abundant and rare microorganisms with sediment depth and the processes behind these changes, we partitioned the β‐diversity into richness and replacement components. This can help to distinguish taxonomic changes that result from loss of surface species during burial processes and microbial diversification, and that from selection of surface populations. In a simplified, competing structural model proposed by Wurzbacher et al. (2017), the richness effect represents a community exclusively consisting of species filtered from the surface taxa, resulting in the transition from a complex and rich community at the surface to an increasingly depauperate community in deeper sediments. The replacement effect instead resembles a structured community of niche specialists at various layers that are well adapted to the specific environmental conditions. Our results showed clear changes of these distinct kinds of effects in marine sediments. For the abundant subcommunity, the richness effect consistently accounted for a high proportion of β‐diversity, indicating that the abundant species in deep sediment are a subset of surface populations primarily resulted from the process of selection. This was consistent with previous reports that selection is the major process determining the assembly of abundant communities (Starnawski et al., 2017). However, the observation of a decreasing trend in the richness effect in deepest layers (690 and 790 cm) for both abundant and moderate subcommunities indicated that the filtering effect might eventually be weakened as sediment depth increased. For the rare subcommunity, by contrast, the replacement effect showed overwhelming dominance below 190 cm, indicating loss of surface taxa and existence of novel species compared with those in the surface. Rare microorganisms are susceptible to local elimination through environmental or biological processes that select for the survival of abundant organisms (Pedrós‐Alió, 2012). During burial, many rare species tended to be filtered out due to their low abundance, resulting in the high replacement effect when compared with the surface. Meanwhile, distinct rare communities derived from terrestrial input and persisting in deep sediments for a long time as inactive seed banks might also be a possible explanation for the elevated replacement effect. Other processes, such as diversification, could occasionally occur in deep sediments, but may not be the dominant determinant for species replacement.

The β‐partitioning patterns above 30 cm in SCS sediments were highly similar to those reported by Wurzbacher et al. (2017) in lake sediment, featuring a consistently high replacement effect and increased richness effect. However, in contrast to their assumption, our results from deeper sediments revealed that the richness component did not show a further increase in a linear manner. In fact, the richness effect was abated below 30 cm and replacement took over in deeper layers, and such changes were predominantly determined by the rare subcommunity. Given the different environmental conditions such as salinity and sulphate concentrations, as well as sedimentation process in marine and freshwater sediments, whether this changing pattern of β‐diversity components can be applied to distinct sediment types is yet to be proven.

Although these SCS samples were collected on a small regional scale, site‐to‐site community variations were observed. As indicated by the increased horizontal β‐diversity, microbial communities at the same depth exhibited greater differences in deep sediments than in shallow sediments, in line with the results shown between two coastal habitats of the Mediterranean Sea (Luna et al., 2013). Despite this moderate increase in total horizontal β‐diversity with depth, the underlying changes in the two β‐diversity components were more pronounced and showed opposing trends, especially for the rare subcommunity. The replacement effect was stronger than the richness effect especially in the shallow sediments, reflecting species loss and gain by stochastic processes such as dispersal from seawater columns and distinct terrestrial inputs. However, during the burial processes over time, these communities have been confined to specific populations that can survive in the deep seafloor, resulting in an increasing contribution of the richness effect to horizontal β‐diversity. These results are the first to reveal the detailed contributions of spatial turnover (replacement) and nestedness (richness) underlying the total amount of β‐diversity at certain depths of marine sediments.

Co‐occurrence patterns indicated reduced microbial interactions in the deep seafloor and the importance of the rare subcommunity

Co‐occurrence networks can potentially provide novel and unique perspectives on microbial interactions and ecological assembly rules beyond those of diversity and composition. The correlation‐based network was divided into three major modules with species occupying either shallower (oxic) or deeper (anoxic) sediments, which clearly reflected niche differentiation related to sedimentary redox environment. Each module comprised a group of diverse microbial taxa, with large numbers of Proteobacteria, Chloroflexi and Planctomycetes in Modules 1, 2 and 3, respectively. A similar network pattern has been reported for the South Yellow Sea (SYS) and East China Sea (ECS) sediment samples from the top 30 cm (Qiao et al., 2018), in which a separation of microbes from shallower and deeper sediments into two modules was observed. Here, we created a co‐occurrence network involving deeper subseafloor samples and more OTUs (especially rare OTUs), which may represent a general model for microbial co‐occurrence patterns in marine sediments.

The topological features of the network revealed that the abundant OTUs had more connections with other OTUs and were more frequently located in central positions of the network compared with rare OTUs located relatively distantly (Ma et al., 2016), emphasizing the vital role of abundant OTUs in microbial networks (Zhang, Hou, et al., 2020). However, rare OTUs accounted for a considerable proportion of nodes (45.4%) and keystone OTUs (69.5%), indicating that rare OTUs constituted an indispensable part in the microbial network. Previous network analysis using selected OTUs based on relative abundance or read number to reduce complexity might have neglected the impacts of these rare species (Buongiorno et al., 2019; Liu et al., 2014, 2020). A higher clustering coefficient of rare OTUs indicated that they were more likely to present in an interconnected “small world,” emphasizing that rare OTUs contributed to the complexity of networks (Masuda et al., 2018). Rare OTUs with high degree were primarily found in Module 1, suggesting that rare OTUs might be more important in supporting community function and stability in shallow sediments than in deep sediments. Among them, OTUs belonging to Acidobacteria and Phycisphaerae in Planctomycetes are involved in the decomposition of various biopolymers (Hoshino et al., 2020; Quaiser et al., 2008), and some of these OTUs were identified as keystone OTUs (Table S3). These rare keystone taxa were considered to exert a pronounced influence by participating in specific single steps of broader processes such as OM decomposition (Banerjee et al., 2018). However, the function and metabolic features of most of these keystone OTUs remain largely unknown, such as the uncultured NB1‐j in Deltaproteobacteria. Their potential roles in sediment microbial networks have yet to be revealed. It was noticeable that the roles of rare OTUs should be interpreted with caution, as they might represent dead or dormant cells and extracellular DNA (cysts, spores and inactive seed banks) (Vuillemin et al., 2017), especially those belonging to Bacteroidetes, Firmicutes, Actinobacteria and Planctomycetes. They might have no functional roles in the network but only coexist with other species (Lennon & Jones, 2011), and their co‐occurrence reflected distal inflows from the Pearl River down to the continental slope. Despite this, the importance of rare taxa has been emphasized in the networks of various ecosystems including soil, rhizosphere and reservoir systems (Lupatini et al., 2014; Shi et al., 2016; Xue et al., 2018), and here we provide evidence suggesting a high proportion of rare species in sediment microbial networks and some of them may be involved in OM decomposition.

While the network of samples from all depths largely reflects shared niche preferences due to geochemical zonation in sediment, networks of separate depth groups may better represent microbial interactions at specific depths. Networks with reduced complexity in deep sediment implied that microbes in the deep subseafloor may interact less with each other, and this might be due to (i) the reduction of microbial abundance and the loss of many species compared with the surface environment; and (ii) the slowed growth rate and metabolism of these microorganisms with some of them entering the dormancy state to maintain survival (Jorgensen & Marshall, 2016). For instance, although Dehalococcoidia were highly enriched in deep sediment, the correlations involving Dehalococcoidia in deep sediments were less than those in the surface (189 and 238 edges, respectively). In the surface, Dehalococcoidia had more correlations with other microbial clades, such as members of the Gammaproteobacteria, Gemmatimonadetes and Phycisphaerae, and such corrections across different clades decreased in deep sediments. However, associations between Dehalococcoidia and sulfate‐reducing bacteria in Deltaproteobacteria were identified in deep sediments, suggesting that sulphate‐reducing activity may enhance dechlorination by Dehalococcoidia (May et al., 2008). Other studies have suggested that carbon inputs and shifts in carbon resource increased network complexity (Farrer et al., 2019; Shi et al., 2016), which could possibly explain the drastic increase of interactions in the network of 30–90 cm, as increased TOC and distinct δ13C was detected at ~50 cm. These results indicated that terrestrial inputs may also influence the structure of microbial co‐occurrence networks.

CONCLUSIONS

Through a comprehensive investigation of the microbial communities in surface and subseafloor SCS sediments, we showed distinct patterns of β‐diversity partitioning and co‐occurrence along sediment depth down to ~8 m. Generally, selection was the dominant processes determining the β‐diversity partitioning patterns. The high replacement effect between deep sediment and the surface layers may suggest loss of rare species during burial as well as dispersal of species with terrestrial inputs. Co‐occurrence network analysis further showed that microbes in the deep subseafloor tended to survive independently with less syntrophic relationships, due to the reduced abundance and metabolic activity. However, rare species accounted for a large proportion in the network, with several keystone taxa that may involve in OM decomposition. Our results expand current knowledge on the ecological significance of rare species related to microbial diversity and interactions in marine sediments.

AUTHOR CONTRIBUTIONS

X.‐H.Z., J.L. and Y.Z. conceived the study. Y.Z. and J.L. analysed the data and wrote the manuscript. P.Y. analysed the environmental factors. C.S. performed DNA extraction and qPCR experiments. X.S. performed sampling and S.L. provided the opportunity for sampling. All authors edited and approved the final manuscript.

ETHICAL APPROVAL

Not applicable.

CONSENT FOR PUBLICATION

Not applicable.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.