All samples were collected from the Qi’ao Island (113°36′∼113°39′E, 22°23′∼22°27′N) located on the western side of the mouth of the Pearl River Estuary in Guangdong province, Southern China. The annual average temperature and precipitation are 22.4°C and 1,964 mm, respectively. The area of this island is 24 km² and the sediment is mainly composed of clayey silt (Wu et al., 2017). We collected five field samples from inshore sites (randomly named CM, SAP, and KO sites) to offshore sites (randomly named MF, SA sites) in July 2018 (Supplementary Figure S1). In each sites, we established six about 2 m × 2 m mudflats and a sediment core was sampled in each plot using a custom-made sampler (8 cm in diameter, 10 cm in depth). In total, 30 samples (5 sites × 6 replicates) were collected, then all samples were sealed in sterile bags and transported to the laboratory in ice boxes within 12 h. Each sample was divided into two sub-samples: one was kept at 4°C for sediment properties analysis, and the other was stored at –80°C for DNA extraction within 24 h.

The sediment properties were measured as previously described (Yu et al., 2020). Briefly, the temperature in situ was determined by a temperature meter (PHBJ-260, INESA instrument, China). The fresh sediment (2.0 g) was centrifuged at 3500 rpm for 10 min. The supernatant was collected to measure the contents of ammonium nitrogen (NH₄⁺-N), nitrate nitrogen (NO₃^–-N), and nitrite nitrogen (NO₂^–-N) by a continuous flow analyzer (SKALAR San++, Skalar Co., Netherlands) (Hong et al., 2019). The water content was measured after drying 5.0 g fresh sediments at 105°C to a constant weight. The sediments before measuring pH, salinity, total carbon (TC) and total nitrogen (TN) were dried at 65°C to a constant weight. TC and TN were detected with 40 mg dry sediment using an elemental analyzer (Vario TOC, Elemental, Germany). The pH was determined by a pH meter (SevenCompact210, Mettler-Toledo, United States) with water and sediment ratio of 2.5:1. Salinity was measured by a salinity meter (EUTECH SALT6+, Thermo Scientific, United States) with water and sediment ratio of 5:1.

Microbial community DNA was extracted from the sediment using a modified sodium dodecyl sulfate extraction method (Zhou et al., 1996) combined with a Power Soil DNA kit (Mo Bio Laboratories, Carlsbad, CA, United States) following the manufacturer’s instructions. Briefly, sediment samples had been fully ground with liquid nitrogen before DNA extraction. Then, DNA was extracted from 0.5 g sediment by using the same method described above. The quality of DNA was determined by a NanoDrop spectrophotometer (NanoDrop One, Thermo Fisher Scientific, United States) with ratios of 260/280 and 260/230, and then diluted into 2 ng μL^–1 for subsequent PCR amplification.

The primers 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) were used for amplifying of V3-V4 hypervariable region of bacterial 16S rRNA gene with an expected fragment size of 468 bp (Wang et al., 2019), and the primers Arch349F (5′-GYGCASCAGKCGMGAAW-3′) and Arch806R (5′-GGACTACVSGGGTATCTAAT-3′) were used to amplify V3-V4 hypervariable region of archaeal 16S rRNA gene with a correct size of 458 bp (Frank et al., 2013). Amplicon sequencing was performed with an Illumina HiSeq platform 2500 (Illumina, Inc., CA, United States) (2 × 250 paired ends) in Biomarker Technologies Corporation, Beijing, China. Raw 16S rRNA gene sequences data were deposited in the NCBI Sequence Read Archives¹ with accession numbers PRJNA666027 and PRJNA666234, respectively.

The trimming process of rRNA gene amplicon sequences were conducted by a publicly available pipeline² (Yu et al., 2021b). In brief, each sample was assigned based on unique barcode and the primers were removed using Cutadapt (Martin, 2011). Low-quality sequences (length < 350 bp, with the ambiguous base “N,” and average base quality score < 30) and poorly overlapped were filtered by Trimmomatic (Bolger et al., 2014). Potential chimeric reads were identified and removed from the dataset using UCHIME in reference database mode (Edgar, 2013). The high-quality sequences were clustered into operational taxonomic units (OTUs) at a 97% identity threshold by UPARSE (Edgar, 2013) and singleton OTUs were not involved in subsequent analysis. Resampling was performed for bacterial sequences (53,166) and archaeal sequences (76,548) to make sure that the sequencing depth is same for all the samples (Supplementary Table S1). Taxonomy was identified by the Ribosomal Database Project classifier at a confidence level of 80% (Wang et al., 2007; Zhong et al., 2017).

To understand possible interactions among different taxa of microbial community and identify potential keystone taxa, the networks of bacterial and archaeal communities were constructed by the Molecular Ecological Network Analysis (MENA) pipeline³ (Deng et al., 2012). In each network, the size of each node is proportional to the number of connections (i.e., degree). Nodes represented individual microbial taxa and edges indicate the pairwise correlations between nodes. The connectivity of each node was estimated by its within-module connectivity (Zi) and among-module connectivity (Pi). The nodes with either a high value of Zi or Pi were identified as potential keystone taxa, including module hubs (Zi ≥ 2.5, Pi < 0.62), connectors (Zi < 2.5, Pi ≥ 0.62), and network hubs (Zi ≥ 2.5, Pi ≥ 0.62) (Guimera and Amaral, 2005). The final visualization of networks was generated by Gephi software (Bastian et al., 2009).

Structural equation modeling (SEM) is theory-oriented and capable of evaluating network of causal hypotheses via testing simultaneous influences rather than individual (bivariate) causes (Grace et al., 2012; Eisenhauer et al., 2015). A large amount of studies confirm that it is a powerful tool to reveal the synergistic relationships between multiple environmental factors and the diversity or keystone taxa of microbial communities (Chen et al., 2019; Mamet et al., 2019; Li et al., 2020). To untangle the relationship between the two most important sediment parameters (pH and salinity), keystone taxa and microbial diversity, we conducted SEM analysis in AMOS 22.0 (AMOS IBM, United States), and the fitness was evaluated by a non-significant chi-square test (χ², P > 0.05), the goodness-of-fit index (GFI) and the root mean square error of approximation (RMSEA) (Byrne, 2016).

To quantitatively evaluate the contribution of major ecological processes to bacterial and archaeal community assembly, we used a null-model-based framework (Stegen et al., 2013; Dini-Andreote et al., 2015). The dynamics of phylogenetic and taxonomic diversity were measured with the null model-based beta nearest taxon indices (βNTI) and Bray–Curtis-based Raup–Crick (RC_Bray) metrics, respectively (Stegen et al., 2013). The βNTI > 2 was interpreted as heterogeneous selection, and βNTI < –2 was denoted as homogeneous selection. When the | βNTI| < 2 and RC_Bray < –0.95, the community assembly is assembled primarily by the homogenizing dispersal. Whereas | βNTI| < 2 and RC_Bray > 0.95, the community is mainly driven by the dispersal limitation. When the | βNTI| < 2 and | RC_Bray| < 0.95, that means no dominant assembly process (Stegen et al., 2013). All above ecological processes analyses were performed in a publicly available pipeline (see footnote 2).

Before statistical analysis, normal distribution was examined and non-normal data were transformed to normal ones using Blom’s formula under equivalent command Rank Cases (Li et al., 2016). The α-diversity was evaluated by Chao1 and Shannon indices, and the β-diversity was estimated by Bray–Curtis distances. The community dissimilarity was visualized using principal co-ordinates analysis (PCoA) based on the Bray–Curtis distances, and the statistical significance was tested by permutational multivariate analysis (PERMANOVA) and analysis of similarities (ANOSIM) using adonis and anosim functions in the R package vegan (Yu et al., 2021a). The relationships between sediment parameters and microbial communities were analyzed by Mantel tests and redundancy analysis (RDA). Variation partition analysis (VPA) was used to quantify the relative contributions of the physicochemical factors to bacterial and archaeal communities with the vegan package (Yu et al., 2021b). Differences among sites were evaluated by one-way analysis of variance (ANOVA) in SPSS 18.0 software.

To determined variations of the environmental factors, the sediment physicochemical parameters (pH, salinity, temperature and moisture content) and nutrient contents (TC, TN, NH₄⁺-N, NO₂^–-N, and NO₃^–-N) from offshore sites (MF, SA) and offshore sites (CM, SAP, and KO) were analyzed (Supplementary Figure S1). The results showed that pH was significantly (P < 0.05) higher in MF sites than in CM and KO sites, while the salinity was significantly (P < 0.05) higher in MF, SA sites than in CM and SAP sites (Table 1), indicating a significant increased sediment pH and salinity in offshore sites. For nutrient contents, the concentrations of TC, TN, and NH₄⁺-N were significantly (P < 0.05) higher in CM, SAP, and KO sites than in MF and SA sites, indicating a high nutrient level in inshore sites.

We obtained a total of 2,828,732 and 3,827,293 high-quality bacterial and archaeal sequences, respectively (Supplementary Table S1). The Chao1 value and Shannon indices of both bacterial and archaeal communities were highest in KO sites, followed by MF, SA, CM sites, and lowest in SAP sites, while there were no significant (P < 0.05) differences between MF, SA, and CM sites (Figures 1A–D). The β-diversity (Bray–Curtis dissimilarity) within bacterial communities was significantly (P < 0.05) lower in MF sites than other sites (Figure 1E), and the β-diversity of archaeal communities was significantly (P < 0.05) lower in MF, SA, and CM sites than in SAP and KO sites (Figure 1F). The results indicated that α- and β-diversity of prokaryotic communities generally did not change significantly (P > 0.05) between offshore sites, but changed significantly (P < 0.05) among inshore sites.

Principal coordinates analysis (PCoA) showed that bacterial and archaeal communities were well separated across five sites (Figures 2A,B). Two non-parametric tests (Adonis and Anosim) showed consistently significant (P < 0.001) variations in bacterial and archaeal communities among five sites, indicating that sediment factors shifted the structure of prokaryotic communities. The bacterial communities were dominated by phyla Proteobacteria (36.4%), Chloroflexi (21.1%) and Acidobacteria (7.8%) (Supplementary Figure S2A). The majority of archaea belonged to phyla Crenarchaeota (61.3%), Euryarchaeota (32.4%), and Parvarchaeota (6.0%) (Supplementary Figure S2B). The dominant bacterial genera included Nitrospira (1.3%), Rhodoplanes (0.8%), and Desulfococcus (0.6%) (Supplementary Figure S3A). The relative abundance of Nitrospira, Rhodoplanes, and Desulfococcus was significantly (P < 0.05) higher in MF, SAP, and KO sites, respectively, than other sites (Figure 2C). For archaeal communities, the majority of classes belonged to Bathyarchaeota (34.3%), Methanomicrobia (27.9%), Thaumarchaeota (25.7%), Parvarchaea (5.8%), and Thermoplasmata (3.4%) (Supplementary Figure S3B). Specifically, the order of Methanosarcina in SA and CM sites was significantly (P < 0.05) higher than that in MF, SAP, and KO sites (Supplementary Table S2). The top 11 genera belonged to methanogens, including Methanosaeta (19.3%), Candidatus Methanoregula (2.7%), and Methanolinea (2.7%) (Supplementary Figure S3B). The relative abundance of Methanosaeta and Candidatus Methanoregula was higher in SA and CM sites than other three sites (Figure 2D), indicating that methanogens dominated in SA and CM sites.

To identify potential interactions and keystone taxa of microbial communities in coastal wetland sediments, we constructed co-occurrence networks of bacterial and archaeal communities, respectively (Supplementary Figure S4 and Supplementary Table S3). The proportions of positive correlations between nodes in both bacterial and archaeal networks were 64.7 and 99.3%, respectively (Figures 3A,C), indicating that cooperation was fundamental in bacterial and archaeal communities. The microbial community co-occurrence networks showed strong positive and negative links with salinity and pH, respectively (Figures 3B,D). For bacteria, we detected 16 keystone OTUs and the most dominant taxa (relative abundance > 0.1%) were derived from Gammaproteobacteria (OTU5 and OTU124) and Betaproteobacteria (OTU12 and OTU82) (Supplementary Figure S4A). The relative abundances of the four OTUs were significantly (P < 0.05) higher in MF and SA sites than in CM, SAP and KO sites (Supplementary Figure S4C). For archaea, we identified four keystone OTUs, including two module hubs belonging to class Bathyarchaeota (OTU7884) and Heimdallarchaeota (OTU28), and two connectors derived from Bathyarchaeota (OTU37 and OTU100) (Supplementary Figure S4B). The relative abundances of these four OTUs were significantly (P < 0.05) higher in MF and SA sites than in CM and SAP sites (Supplementary Figure S4D). Especially, the top four keystone OTUs of bacteria were significantly (P < 0.05) positive correlation with pH and salinity, and the top four keystone OTUs of archaea were significantly (P < 0.05) positive correlation with salinity (Supplementary Figure S5). Collectively, the bacterial and archaeal keystone taxa were mainly distributed in offshore sites and had close links with pH and salinity.

To reveal the relationship between sediment factors and the microbial community structure, their correlations were analyzed by Mantel tests (Supplementary Table S4) and Pairwise Spearman’s correlation (Supplementary Figure S6). The factors with significant (P < 0.05) influence were further assessed by RDA. The results showed that pH, TC and NH₄⁺ were the major factors to explain bacterial community variations (Figure 4A), while pH were strongly correlated with archaeal communities (Figure 4B). In addition, the salinity showed more stronger effects on the prokaryotic community in offshore sites than that in inshore sites (Figures 4A,B). VPA further indicated that sediment pH and salinity clearly contributed to a large proportion of variation to bacterial (12.3%, 40.6%) and archaeal (22.7%, 33.8%) community, respectively (Figures 4C,D). Also, linear regression analysis showed that the α-diversity of bacterial and archaeal communities was significantly (P < 0.05) positively associated with salinity (Supplementary Figures S7A,B). The β-diversity showed a significantly (P < 0.05) positive correlation with pH and salinity (Supplementary Figure S7C). In addition, pH and salinity showed strong effects on the relative aboundance of bacterial and archaeal taxa. For example, pH was significantly (P < 0.05) associated with the top four bacterial phyla (Supplementary Figure S8A), the top three archaeal classes (Supplementary Figure S8B), and the top four archaeal genera (Supplementary Figure S9B). Collectively, these results indicated that sediment pH and salinity could influence the diversity, composition and structure of bacterial and archaeal communities.

To better understand the effects of pH and salinity on the diversity and keystone taxa of microbial communities, we further constructed SEM. For bacterial communities, the result showed that sediment pH and salinity had significant (P < 0.01) positive effects on bacterial keystone taxa (Figure 5A). Also, salinity had significant positive effects on the α- and β-diversity of bacterial community (standardized effects = 0.47, P < 0.01; 0.38, P < 0.05). However, sediment pH could affect bacterial keystone taxa to indirectly inhibit the β-diversity of bacterial community. For archaeal communities, sediment salinity showed significantly (P < 0.01) positive effects on the α-diversity and keystone taxa (Figure 5B). Moreover, salinity had indirect effects on the α-diversity of archaeal community through influencing their keystone taxa (P < 0.001). Sediment pH had significant negative effects on the α-diversity of archaeal communities (standardized effects = –0.25, P < 0.05). Overall, these results indicated that sediment pH and salinity could affect the diversity and keystone taxa of bacterial and archaeal communities.

The abundance-weighted beta nearest taxon index (βNTI) showed that selection was one of essential processes in bacterial and archaeal communities assembly (Figure 6A). For bacterial community, the modified Raup–Crick (RC_bray) played a principal role in inshore and offshore sites (Figure 6B). However, heterogeneous selection (21.6%) and homogeneous selection (6.0%) were higher in offshore sites than that in inshore sites, indicating that deterministic processes were more important in offshore sites. The relative contribution of stochastic process to the archaeal community was the highest in all sites (93.4%) (Figure 6C). The proportion of homogeneous selection to the archaeal community was the highest in inshore (95.6%) and offshore sites (81.9%), while dispersal limitation and homogenizing dispersal almost could be ignored. Therefore, stochastic process was major driving factor in the assembly of prokaryotic communities in offshore sites.