Spatiotemporal Dynamics of Bacterial Taxonomic and Functional Profiles in Estuarine Intertidal Soils of China Coastal Zone

Bacteria play an important role in regulating carbon (C), nitrogen (N), and sulfur (S) in estuarine intertidal wetlands. To gain insights into the ecological and metabolic modes possessed by bacteria in estuarine intertidal wetlands, a total of 78 surface soil samples were collected from China’s coastal intertidal wetlands to examine the spatial and seasonal variations of bacterial taxonomic composition, assembly processes, and ecological system functions through shotgun metagenomic and 16S rRNA gene sequencing. Obvious spatiotemporal dynamic patterns in the bacterial community structure were identified, with more pronounced seasonal rather than spatial variations. Dispersion limitation was observed to act as a critical factor affecting community assembly, explaining approximately half of the total variation in the bacterial community. Functional bacterial community structure exhibited a more significant latitudinal change than seasonal variability, highlighting that functional stability of the bacterial communities differed with their taxonomic variability. Identification of biogeochemically related links between C, N, and S cycles in the soils showed the adaptive routed metabolism of the bacterial communities and the strong interactions between coupled metabolic pathways. Our study broadens the insights into the taxonomic and functional profiles of bacteria in China’s estuarine intertidal soils and helps us understand the effects exerted by environmental factors on the ecological health and microbial diversity of estuarine intertidal flats.


Introduction
Estuarine intertidal wetlands refer to the transition-related areas between the land and coastal seas, acting as hotspots of biogeochemical cycles [1]. In dynamic environments, bacteria are critical to the functioning of ecological systems. In particular, the functionally diverse microbial community contributes to biogeochemical transformations of elements such as carbon (C), nitrogen (N), and sulfur (S) [2]. The study of microbial spatiotemporal patterns and related environmental factors is vital to understanding the mechanisms that generate and maintain species diversity [3][4][5][6]. Previous reports demonstrated that, at a small spatial scale, the microbial community in estuarine intertidal soils changes tightly with environmental gradients, while at a large spatial scale (such as thousands of kilometers), the microbial community generally shows obvious geographic patterns, which are not only mediated by variations in regional environmental characteristics but also correlated with spatial factors [5,7,8]. One of the most common biogeographic patterns is the distance-decay relationship, which refers to decreasing community similarity with increasing geographic distance [9,10]. This relationship can be influenced by both deterministic and stochastic processes [11][12][13]. Their relative importance in governing bacterial community assembly varies across diverse habitats. For instance, stochasticity exerts a stronger effect than determinism on the bacterial community assembly process in subtropical marine ecosystems [14], whereas the deterministic process is more important in lake and agricultural soils [10,15]. It is thus important to identify the major assembly process governing the variation of bacterial communities.
Microorganisms in a specific ecological niche form complex interaction webs, and their interactions are important aspects in terms of maintaining the stability and diversity of microbial communities [16,17]. Correlation-based network analysis has been extensively utilized to infer microbial interactions in marine water [18], activated sludge [19], and grassland soil [20]; the results of which provide important information about microbial ecology. However, the spatiotemporal patterns of microbial co-occurrence relationships in the transitional zones between the land and coastal seas remain unclear, and the relationships between dynamics of microbial interaction and variations of community structure also need further exploration.
Over the past several decades, most of China's estuarine intertidal wetlands have been affected by human activity, and there has also been an increase in the frequency and severity of extreme runoff events caused by climate change [21,22]. As notable contributors to biogeochemical cycles in estuarine intertidal wetlands, bacterial communities and their metabolic pathways may have important ecological implications [23]. Compared to other terrestrial and marine ecosystems [9,10,15], we know less about the spatiotemporal dynamics of the bacterial community, species interaction, assembly processes, and functional potential in soils of estuarine intertidal wetlands. Direct sequencing of environmental DNA, avoiding the cultivation bottleneck, can potentially characterize the genetic capabilities in novel and uncultivated organisms [24,25]. Building on these metagenomic approaches, we hoped to reveal the genetic basis of microbial communities involved in biogeochemical processes, thereby deciphering the ecological functions of bacterial communities in estuarine intertidal wetlands.
Here, the spatiotemporal dynamics, co-occurrence patterns, and assembly processes of bacterial community and functional structures in China's estuarine intertidal soils were investigated. In addition, this study also focused on the microbial functions involved in element (including C, N, and S) metabolisms, and the relationships of soil environmental variables with taxonomic and functional compositions were compared to elucidate potential mechanisms driving the changes in microbial communities.

Field Sampling and Environmental Characteristics
The coastal wetlands of China cover an area of more than three million square kilometers across tropical, subtropical, and temperate climatic zones. In this study, 10 estuarine intertidal wetlands were totally selected along the coastal zone of China, which are located in the estuaries of Liaohe River (LH), Haihe River (HH), Yellow River (YR), Sheyanghe River (SYH), Changjiang River (CJ), Oujiang River (OJ), Minjiang River (MJ), Jiulongjiang River (JLJ), Zhujiang River (ZJ), and Beibuwan Gulf (BBW), respectively ( Figure S1), ranging from high-latitude (CJ to LH; 30°N to 41°N) to low-latitude areas (BBW to OJ; 21°N to 28°N). From March to April and August to September 2019, 3-5 sampling sites were selected in intertidal flats of each estuary, respectively. At the respective sites, soil samples (0-5-cm deep) were taken during ebb tides from fifteen to twenty sampling plots within an area of approximately 100 m 2 , homogenized to generate a composite sample, and then immediately transported to the lab at a 5 °C cooler. In total, 78 soil samples were obtained (Supplementary Table S1). One aliquot of each composite sample was stored in liquid nitrogen for microbial analyses, and a second aliquot was stored at -20 °C for analyses of sediment physicochemical parameters involving salinity, pH, exchangeable nitrite (NO 2 − -N), nitrate (NO 3 − -N), ammonium (NH 4 + -N), microbially oxidizable ferrous iron (Fe 2+ ), reducible ferric iron (Fe 3+ ), sulfide concentrations, total organic carbon (TOC), and moisture content (MC). Measurements of the respective physiochemical parameters of the soil samples were performed in triplicate, and specific analysis methods and related data are shown in Supplementary Tables S1.

DNA Extraction and 16S rRNA Gene Amplicon Sequencing
The FastDNA®SPIN Kit (MP Biochemicals, Solon, OH, USA) was used to extract total DNA from soil samples (0.5 g) following the manufacturer's protocol. Based on the 338F/806R primer pairs, bacterial 16S rRNA gene polymerase chain reaction (PCR) assays were conducted. Electrophoresis was used to detect PCR products in the bright band range of 400-450 bp in a 2% (w/v) agarose gel, and the mentioned PCR products were mixed in the identical density ratio and purified based on the Gene-JET Gel Extraction Kit (Thermo Scientific) [7]. Next, we mixed three parallel PCR products in the respective subsample for sequencing on an Illumina platform. by using the UCHIME algorithm [27]. The mentioned sequences contained the unique barcodes in the assignment to the respective sample. The UPARSE software package was adopted to analyze the sequences based on the UPARSE-OTU algorithm and UPARSE-OTUref algorithm. We endowed sequences with a similarity of 97% as identical operational taxonomic units (OTUs). The RDP classifier (version 2.2, http:// sourc eforge. net/ proje cts/ rdp-class ifier/) helped assign the sequences representing the respective OTU to the taxonomic groups against the SILVA database (release 138, http:// www. arb-silva. de) [26]. To equalize sequencing depth, the respective sample was rarefied to 10,151 reads (the minimum sequence number across overall samples). Using Mothur software (version v.1.30.1), alpha diversity indices consisting of community coverage (coverage), species diversity indices (Shannon-Weiner diversity and phylogenetic diversity), and richness (Ace and Chao I) were determined [28]. For beta diversity, principal coordinate analysis (PCoA), nonmetric multidimensional scaling analysis (NMDS), and permutational multivariate analysis of variance (PERMANOVA) were conducted by exploiting Bray-Curtis dissimilarities. ANOSIM analysis (analysis of similarities) and the partial least squares discriminant analysis (PLS-DA, a supervised analysis applied to highdimensional data) were performed to examine the differences between communities.

Metagenome Sequencing and Statistical Analysis
The extracted DNA of the samples was pooled together (as one sample for metagenome sequencing) from one estuarine intertidal wetland in the respective season. Overall, 20 DNA samples were collected from 10 estuarine intertidal wetlands during the field surveys of two seasons. For the respective sample, we placed 120 ng of genomic DNA in 60 μl of buffer (10 mM Tris HCl, pH 8.0) into 1.5-ml RNase/DNasefree microcentrifuge tubes with low binding. To analyze quality, a NanoPhotometer spectrophotometer (IMPLEN, CA, USA) and a Qubit dsDNA Assay Kit in a Qubit 2.0 fluorometer (Life Technologies, CA, USA) were adopted. The NEBNext UltraTM DNA Library Prep Kit for Illumina (NEB, USA) assisted in preparing the DNA library following the recommendation of the manufacturer. For the respective sample, the attribute sequence was added under an index code. The AMPure XP system helped purify the samples, and an Agilent 2100 Bioanalyzer (Agilent Technologies, CA) helped determine the libraries. A cBot Cluster Generation System was adopted to generate clusters, and then, paired-end reads (PE150) were generated with an Illumina HiSeq2500 platform to generate a 10-Gb metagenomic dataset.
For quality filtering, the resulting paired-end sequences were analyzed on the free online platform of the Majorbio Cloud Platform (www. major bio. com). The paired-end reads with low quality (quality score ≤ 38; base N > 10 bp, and overlap length between adapter and reads > 15 bp) were filtered. Megahit (https:// github. com/ voutcn/ megah it) was adopted for assembling the sequences based on the methods of succinct de Bruijn graph [29]. Contigs with lengths over 300 bp were reserved and produced based on scaffold fragmentation for subsequent study. For assembled metagenomes, MetaGeneMark v.2.10 helped predict the open reading frames (ORFs) [30], CD-HIT v.4.5.8 was used to build a nonredundant gene catalog (90% nucleotide identity), and SOAPaligner v.2.21 was adopted for comparison with the nonredundant gene set [31,32]. Based on the reads per kilobase per million mapped reads (RPKM) values, gene abundance was calculated [33]. Next, BLASTP (BLAST Version 2.2.28 +) with the help of the KEGG database (Kyoto Encyclopedia of Genes and Genomes, http:// www. genome. jp/ kegg/) and NR database (e-value 1e-5) was used to functionally and taxonomically annotate C, N, S, and methane (CH 4 ) metabolisms [34]. PCoA and NMDS were used to visualize the pairwise Bray-Curtis dissimilarities between functional and taxonomic communities. To explore the diversity of genes in C metabolism, high-quality reads were searched against the Carbohydrate Active Enzyme Database (CAZy, http:// www. cazy. org). All the raw reads of 16S genes and shotgun metagenomic sequences were deposited in the NCBI SRA under accession identification number PRJNA755846.

Data Analyses
Principal component analysis (PCA) was performed to detect spatial and temporal distributions of physicochemical parameters. Procrustes analysis was conducted to assess congruence between ordinations of microbial community structure and soil environmental parameters by testing the significance of the Procrustes statistics with 999 permutations [35]. Spearman rank correlations were employed to test the effect exerted by environmental properties on the top 50 relative abundances of OTUs. In addition, redundancy analysis (RDA) was used to examine the relationships between microbial community structure and environmental indices. The Mantel test was adopted to test the correlation between two Bray-Curtis distance matrices. Multivariate association with linear models (MaAsLin) was adopted to identify the effects exerted by specific environments or climatic factors on bacterial dynamics [36]. To compare the degree of functional variability versus taxonomic variability within CH 4 , S, and N metabolisms, we examined the coefficients of variation (CV, the standard deviation divided by the mean) of relative abundances (percentage of community abundance based on RPKM result) in each functional and taxonomic level.
The distance-decay rate exhibited by the bacterial communities referred to the slope obtained from a linear least-square regression regarding the association of geographic distance with 1 3 bacterial pairwise Bray-Curtis dissimilarities. The principal coordinates of neighbor matrices (PCNM) procedure assisted in deriving the variable of spatial information from the geographic coordinates [37]. Environmental factors and PCNMs that were significant in accounting for bacterial changes in the RDA were used for further variation partitioning analysis (VPA) to address the relative roles of these two factors, as well as their combined effect. Welch's t-test was used to identify the differences in relative abundance, bacterial taxonomy, and functions between sampling seasons and sites (confidence interval method) [38].
The network assisted in exploring the co-occurrence modes exhibited by bacteria, with the selection of OTUs of the top 100 relative abundances. Spearman's correlation between the OTUs presented statistical robustness with a coefficient greater than 0.6 and a P-value less than 0.01 [7,39]. The interactive platform Gephi assisted in visualizing the network [40]. We also compared the topology of a real network of functional structures with a random network to test the interaction of functional genes. We built 10,000 Erdös-Rényi random networks, in which each edge exhibits the identical probability of being assigned to any node. Null model analysis was carried out using the framework described by Stegen et al. [41,42] based on the β-nearest taxon index (βNTI), classifying the community into underlying drivers of deterministic and stochastic processes. Stochastic processes consisting of homogenizing dispersal (mass effect), dispersion limit, and drift were distinguished by the Raup-Crick metric (RCbray) as described by Liu et al. [10]. The null model expectation was generated using 999 randomizations.

Bacterial Community and Composition Based on 16S rRNA Genes
For all soil samples, a total of 4,183,137 high-quality bacterial 16S rRNA gene sequences and 18,152 OTUs were obtained by high-throughput sequencing. Seven phyla, Proteobacteria, Bacteroidetes, Chloroflexi, Actinobacteria, Acidobacteria, Firmicutes, and Cyanobacteria (relative abundance > 3.50%), accounted for 88.34% of the total classified ( Figure S2). The Good's coverage values varied between 85.21 and 94.81% across samples, demonstrating that the sequencing effort recovered most of the local species diversity. Among the sampling areas, the average Chao I index ranged from 2585.84 to 4452.51. The average Shannon and PD indices varied from 4.60 to 6.63 and from 130.96 to 199.33, respectively (Table S2). The bacterial α-diversity indices (except PD) tended to be higher for soils from low-latitude sites than for those from high-latitude sites and higher in August to September than in March to April (Table 1). However, for most of the sampling sites in this study, the seasonal changes in bacterial α-diversity from the same estuary were not significant (Fig. 1).
The PCoA and NMDS analysis of the bacterial communities at the OTUs level showed more significantly seasonal separations than spatial separations (ANOSIM, P = 0.001) (Fig. 2). However, although not obviously significant, a structural difference was found between high-and low-latitude regions for soil bacteria (PERMANOVA, P = 0.001) ( Figure S3). Furthermore, significant corresponding relationships were identified between the bacterial composition and soil physicochemical parameters (PCA plots) ( Figure S4) (Mantel test, r = 0.3108, P = 0.001). The correlations of environmental parameters and dominant bacterial OTUs in the heatmap indicated that the dominant bacterial taxa were significantly associated with soil physicochemical properties ( Figure S5). Similarly, RDA results demonstrated that all environmental factors measured except sulfide impacted the soil bacterial community structure. In addition, bacterial community structural changes were mainly dependent on the relative abundances of Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria (Fig. 3).

Bacterial Co-occurrence Patterns
A bacterial co-occurrence network was inferred, consisting of the top 100 OTUs (nodes) and 487 correlations (edges) in August to September and 943 edges from March to April. In comparison, the structural properties, including the modularity index (MD), graph density (GD), average degree (AD), and clustering coefficient (CC), were greater in March to April and exhibited significantly closer interconnections than those in August to September ( Table 2, Fig. 4a, b, and Figure S6ab). It was also observed that there were no significant spatial differences in bacterial co-occurrence patterns between the high-and low-latitude sites (Fig. 4b, c and Figure S6cd).

Bacterial Community Assembly Processes
To quantify the contributing effect exerted by distance (as a proxy for dispersal barriers) and soil characteristics on bacterial community changes, VPA was performed by using the RDA model (Fig. 5c). The VPA diagram appeared to show a higher effect of spatial distance, which accounted for 14.14% of the variation, than that of environment factors (9.16%) (Fig. 5a). However, 62.91% of the structural changes could not be explained in this analysis. We further explored the correlation between the bacterial similarity matrix and the mean distance of the sampling sites. One significant negative distance decay of similarity correlation was identified in the linear regression between sampling sites (Fig. 5b). This result indicates the importance of spatial factors in bacterial community variation. The null model analysis suggests that spatial factors play a relatively dominant role in both March to April and August to September (Fig. 5d).
A significantly higher effect of stochastic processes was found in the bacterial community assembly, and dispersal limitation (DL) accounted for 74.15% and 46.75% of the total variation in August to September and March to April, respectively. This was followed by variable selection (VS), which explained 21.70% in August to September and 33.03% in March to April, and homogeneous selection (HS), drift (DF), and homogenizing dispersal (HD), which explained 0.14%, 3.45%, and 0.56% of the total variation in August to September and 15.92%, 3.69%, and 0.61% of the total variation in March to April, respectively.    carbohydrate metabolism (12,607,312 reads), and energy metabolism (10,020,358 reads) ( Figure S7). In this study, the functional community structure based on metagenome annotation was different from the taxonomic community structure based on the 16S rRNA sequence. The functional gene (KEGG level 3) compositions were gathered by complying with geographical distribution (PERMANOVA, P = 0.014) instead of sampling seasons (PERMANOVA, P = 0.339) (Fig. 6). Similar to the 16S rRNA sequence results, the bacterial community structure also had obvious spatiotemporal variation based on the metagenome sequences (PER-MANOVA, P < 0.05) ( Figure S8). This result implied that the taxonomic changes in bacterial communities in estuarine intertidal soils along China's coastal zone could not affect the functional potential. At the same time, we noted obvious latitudinal variations in the functional similarity of bacteria ( Figure S9), suggesting the influence of latitudinal gradients of environmental factors (or obvious climatic gradients) on  the similarity of functional genes. Pearson correlations indicated that mean annual temperature (MAT) as a climate factor was significantly correlated with the functional potential of bacterial communities in estuarine intertidal soils along China's coastal zone (MaAslin, P = 0.0008) ( Figure S10).

Functional Profiles of Main Element Metabolisms
Genes involved in C, N, and S metabolism-related processes were identified. In the study area, the relative abundance of functional genes encoding carbohydrate-active enzymes (CAZymes) was relatively more conserved across seasons (PERMANOVA, P = 0.355) than spatial areas (PER-MANOVA, P = 0.049) (Fig. 7a). In comparison, although their difference was not significant, the functional genes involved in glycosyl transferases (GT) were more abundant in low-latitude areas and in August to September, and the glycoside hydrolase (GH) genes were more abundant in high-latitude areas and in March to April (Fig. 7b). KOs annotated to CH 4 , S, and N metabolisms showed a relatively conservative distribution ( Figure S11). Their taxonomic composition (family levels) in individual functional groups was highly variable compared with their functional structure (KO level) ( Figure S12). The CVs within KO abundances (0.07-0.53, 0.10-1.20, and 0.06-0.79 for CH4, S, and N, respectively) were typically much lower than taxonomic proportions at the family level (0.08-3.05, 0.14-3.13, and 0.10-1.84 for CH 4 , S, and N, respectively) within each functional group. The network based on CH 4 , N, and S metabolism genes was more significantly clustered than random graphs. This result suggests a close collaboration between functional genes (Table S3). The potential metabolic profile maps based on the relative abundance of key genes showed that CH 4 metabolism was likely a vital connection of the C cycle with the cycles of S and N; in the putative coupling mechanisms among CH 4 , N, and S cycles (Fig. 7c), formate and l-serine linked the CH 4 cycle to the N and S cycles. The genes in the CH 4 cycle were largely involved in acetyl-CoA synthetase (K01895, acs), as well as the later conversion of formate into carbon dioxide (K00123, fdoG) and pyruvate metabolism (K01007, pps) (Fig. 7c). However, the relative abundance of genes involved in methanogenesis (mcr genes) and methanotrophy (pmo and mmo genes) (Fig. 7c) was very low. Moreover, in the top 15 pathway modules (module level) of CH 4 , N, and S metabolisms, there were 7 modules with significant spatial differences between high-and low-latitude sites, while only 2 modules had significant seasonal differences (Fig. 7d). This comparison demonstrates that the key element cycle of bacteria had obvious seasonal conservatism in function across the study area. For bacteria participating in vital element cycles, 415 families were annotated in CH 4 metabolisms; 381 and 364 families were annotated in N and S metabolisms, respectively. The composition similarity of taxonomic structure was congruent with their functional structure (Mantel test, r = 0.5643, P = 0.001) (Fig. 8). However, there was a decoupled relationship between the taxonomic and functional α-diversity (based on Shannon index) (r = 0.182, P = 0.06).

Discussion
Determining the biogeographic distribution, assembly mechanisms and ecological functions of microorganisms and their responses to environmental variations is an important issue in microbial ecology. In this study, the bacterial communities in soils of China's estuarine intertidal wetlands covering a temperate to tropical climate gradient were extensively studied. The datasets here help explore the temporal and spatial patterns of bacterial systems in China's estuarine intertidal soils. Among the major species in the bacterial community, Proteobacteria showed absolute superiority among all phyla. This group accounted for approximately half of this bacterial community, and its wide occurrence has also been suggested in various ecological systems (e.g., rivers and marine habitats) and is believed to be the most dominant population in soils globally [43]. Another dominant phylum, Bacteroidetes, showed a more obviously spatial variation, which might be due to their physiological characteristics.
Bacteroidetes are probably one of the vital groups degrading high-molecular-weight organic substances [44]. The occurrence of Bacteroidetes at sites SYH, YR and HH might be attributed to high hydrocarbon-contaminated environments in North China's coast zone during the establishment of industrial cities in North China [45,46]. For most estuary sites, soil bacteria exhibited little difference in α-diversities between sampling seasons (Fig. 1), although the bacterial species composition had obvious separation (Fig. 2). This revealed that the variation of soil erosion and deposition during two sampling seasons along China's coastal zone generally could not significantly influence bacterial α-diversity, but this different mode caused obvious variation of the bacterial community structure (known as β-diversity). In fact, at most sampling sites, strong hydrodynamic forces caused the mixing of surface soils in both vertical and horizontal directions and thus led to the significant variation of species compositions, especially the major species (e.g., phyla of Proteobacteria and Bacteroidetes) which significantly impacted bacterial β-diversity (Fig. 3b). However, the seasonal conservation of α-diversity was probably the result of bacterial long-term environmental adaptation that finally exhibits this equilibrium [9,47].
An obvious community separation according to sampling seasons rather than sampling sites (Fig. 2) was inconsistent with previous studies in that bacterial and archaeal groups varied more significantly with depth or area of soils compared with seasonal variation [10,48]. The present study found that the composition of bacteria was associated with relevant changes in soil physicochemical characteristics ( Figure S4). Results of RDA and Procrustes analysis showed that the significant changes in soil environmental indices in different seasons might significantly impact the composition of the indigenous bacterial communities (Fig. 3a and Figure S4), demonstrating that deterministic processes (filtering by soil environments) acted to shape the bacterial communities in estuarine intertidal soils [49,50]. However, a large fraction of the variations remained unexplained, implying that other reasons (e.g., stochasticity processes in bacterial systems or unmeasured environmental parameters) could also affect the changes in their compositions.
A noticeable negative correlation of bacterial similarity and spatial distance was found, demonstrating the spatial variation in the bacterial community configuration of China's estuarine intertidal soils (Fig. 5b). To understand the mechanisms of microbial community diversity and biogeographic patterns, Stegen et al. [41,42] formulated a framework integrating a null model and phylogenetic data for discerning ecological processes from different aspects. At this large spatial scale along China's estuarine intertidal wetlands, stochastic processes (spatial impact) were found to have a more significant effect on soil bacteria, compared with deterministic processes (environmental impact) (Fig. 5d). Similarly, the VPA results also showed that spatial factors had a slightly greater impact. However, it should be emphasized that previous studies have found that VPA does not correctly predict the environmental and spatial effects on microbial community variations [10,51]. In fact, it has also been reported that microbial assemblages in the river and intertidal soils are more strongly governed by stochastic processes [52]. In the targeted estuarine intertidal soils, physical disturbances, such as groundwater discharge and tidal flow, frequently occur. For this reason, indigenous microorganisms are likely to have evolved to endure the environmental dynamics, which have a weak response to deterministic factors. Indeed, this could also explain the lower effect of environmental changes in the bacterial community composition (Fig. 5a, c and Fig. 3a).
For the co-occurrence network analysis, seasonal dynamics of bacterial community compositions were investigated ( Fig. 4 and Figure S6). The outcomes indicated distinct season-related patterns of their correlations, and the bacterial taxa were more linked to each other in March to April than in August to September. This result was consistent with the observation of archaeal community composition in the marginal sea soils [10]. An explanatory remark is likely that the low spatial connectivity of estuarine intertidal areas in August to September is primarily promoted through high river discharge (generating various solitary patches), making sampling zones more geographically divided and inducing a more dispersive co-occurring mode. Indeed, we found that geographical factors (PCNMs) more obviously impacted the bacterial structures in August to September than in March to April (Fig. 5c). As revealed by the null model, compared with March to April, the geographical isolation (weight of DL) in stochastic processes was enhanced but homogeneous processes (weight of HS) was weakened from August to September, suggesting that the bacterial habitat was more fragmented in August to September than in March to April (Fig. 5d). Therefore, we concluded that although seasonal environmental variations affected the structural differences of bacterial communities, their communities in estuarine intertidal wetlands of China seemed more regulated by seasonal changes in river runoff discharge. Compared with environmental and spatial factors, the dynamics of runoff exerted a more significant impact on the distribution patterns of microbial communities in estuarine intertidal soils along the coastal zone of China.
The bacterial communities exhibited spatiotemporal changes in taxonomic composition, network topological structure, and assembly processes in China's estuarine intertidal soils, raising the question of whether the change in the microbial communities affected their ecological function. Here, we found that the functional composition of communities (functional genes annotated in KEGG) was highly conserved across the sampling seasons (Fig. 6), which was decoupled from their taxonomic composition. This decoupling is common in natural environments and in microcosmic systems, and a consistent decoupling process of functional and taxonomic compositions was found in bioreactor experiments, such as the methane-producing bioreactors in which a great change in taxonomic composition gradually coincided with stable bioreaction performance [53,54]. One study on marine microbes showed a strong link between taxonomic community composition and function across different temporal and spatial scales [55], which was based on metagenomes (all reads) instead of exploiting merely known genes (orthologous genes annotated). Furthermore, robust correlation processes between compositions of functional and taxonomic groups across a robust redox gradient or a single biotope environment have been identified [56]. In fact, in this study, whether based on nonredundant gene sets (all reads) or KEGG databases (annotated orthologous genes), a similar distribution with significant seasonal conservatism was observed (Mantel test, r = 0.2423, P = 0.013, Figure S13). This result means that there was high functional stability for bacteria in estuarine intertidal soils, regardless of which functional database was used. The data here supported the point of redundancy in microbial systems, and the coexistence of multiple distinct taxa performing identical biochemical activities or various microbial species could perform the identical set of enzymatic reactions [57].
The strong seasonal variation in environmental factors contrasted with the composition of bacterial functions, suggesting that the measured physicochemical indices may not be the main factors affecting this functional distribution, which was also supported by statistical analyses ( Figure S10). In addition, we suggest that the ecological processes in the bacterial system cannot explain the functional variation with latitude because of the strong stability of bacterial functional potentials. Metabolic genes largely appeared early in Earth's history, and as presented, they have propagated to multiple clades during geological periods [58]. Although the environment might be selected for specific functions [59], we assumed that the composition of the function bacteria was more likely to be regulated in long term by soil types or climate. Here, we found that the mean annual temperature (MAT) strongly affected bacterial function levels ( Figure S10), suggesting that the geographical pattern of bacterial function in estuarine intertidal flats of China was the result of long-term bacterial adaptation. Estuarine intertidal ecosystems in this study integrated the changes in the surrounding landscape and atmosphere, and the results highlighted that bacterial functional profiles from high to low latitudes were vulnerable to climate change.
Furthermore, this study focused on the functional characteristics of the key element (C, N, and S) cycles, as well as their metabolic coupling. Here, we found stable proportions based on the abundance of metagenome sequences in CAZymes. Organic matter contained in the surface soils primarily originated from plant polymers (e.g., hemicellulose and cellulose). The degradation of these polysaccharides to oligomerization and monomer sugars was critical to the microbial decomposition of soil organic substances [60]. The distributions of organic matter in estuarine intertidal wetland soils of China displayed a declining trend from high-to low-latitude sites, and the advantageous environmental factors for microbial growth in low-latitude areas were likely to facilitate the comparatively abundant property exhibited by genes that decompose organic matter [61][62][63]. However, there was no significant spatiotemporal variation in CAZyme-affiliated gene abundances except for the presence of spatial clustering (Fig. 7a, b). This result might be due to the diverse and rich C pools in China's estuarine intertidal areas [61].
With the degradation of organic C, together with the consumption of oxygen, various advantageous anaerobic respiration procedures were performed by soil microorganisms, which contributed to CH 4 -, N-, and S-transformation 1 3 processes. For the CH 4 cycle, in this study, a high abundance of the acs and fdoG genes (K01895 and K00123) showed the underlying superiority of the conversion processes of formate to acetyl-CoA and further to acetate or pyruvate metabolic processes. In addition, a high occurrence of the glyA and pps genes (K00600 and K01007) involved in l-serine generation was found ( Figure S11). It is noteworthy that these top genes in CH 4 metabolic processes were distributed primarily in the channels under the coupling of other metabolic processes (e.g., S metabolic process, N metabolic process and C-fixing process), revealing that methane metabolism may be a possible bridge connecting the C cycle to N and S cycles in estuarine intertidal soils. Previous studies indicate that in coastal ecosystems with the riverine input of considerable organic detritus and significant primary production, microbial sulfate reduction is considered a dominant pathway mineralizing organic C and is likely to contribute to over half of the total mineralization [64,65]. As exhibited in the metagenome data, sulfate reduction in China's estuarine intertidal soils was focused mainly from sulfate to adenylyl sulfate (APS) and finally to sulfide. Next, it flowed into sulfate assimilation as an amino acid or to polysulfide existing in soils. Notably, polysulfides refer to the only bioavailable sulfate pool in organisms [66]. For this reason, cysteine and polysulfide biosynthesizing potential boosted estuarine intertidal ecological systems to reduce the adverse effects exerted by sulfide on other organisms and offer them bioavailable and flexible sulfate substances [67,68]. For the N-cycle, the large genomic potential for nitrite/nitrate transformation in the N-cycle pathways may be explained by reliably high loads of N in such severe eutrophic zones. Notably, a higher abundance of the gene in the process of dissimilatory nitrate reduction to ammonium (DNRA) may act as a possible factor resulting in frequent ecological risks in estuarine intertidal ecosystems (e.g., serious eutrophication and harmful algal blooms) [69,70]. Furthermore, network analysis showed that the tight coupling between functional genes of the key element cycles could ensure the highly efficient functioning of CH 4 , S, and N metabolic pathways. In brief, the coupling of biogeochemically related processes displayed the flexible and active properties of metabolisms exhibited by the bacteria in the S, N, and C cycles, promoting the high efficient transformation processes of elements in soils capable of ensuring high productivity and ecological health in estuarine intertidal systems.
Nonsignificant seasonal changes in taxonomic compositions of identified CH 4 -, N-, and S-cycling bacteria were found in the soils, which were consistent with the composition of each predicted protein (functional gene) (Fig. 8). Hence, this study emphasized that the community of CH 4 -, N-, and S-cycling bacteria in estuarine intertidal soils depended on the stability of the functional community composition instead of dynamic environmental factors. However, this result did not mean that there was no functional redundancy in a specific functional group. In this research, a more variable taxonomic composition across sampling sites suggested functional redundancy for functional communities or at least exhibited partial redundancy, in contrast to the comparatively conserved functional composition ( Figure S9 and Figure S10) [55]. Numerous microorganisms able to oxidize hydrogen could jointly exist in groundwater, and considerable oxygenic photoautotrophs could jointly exist on the ocean surface [71,72]. Taxonomic variability should not imply differences in community function, and it is generally determined by various approaches to foraging, attachment to particles through biofilm formation, employment of various nitrogen pools, and resistance to a range of antibiotics [55]. Thus, the high diversity of functional bacteria in estuarine intertidal soils could significantly boost element cycles in terms of ecological function.

Conclusions
This study demonstrated obvious spatiotemporal patterns of the taxonomic composition of bacteria in China's estuarine intertidal soils, and the spatial distance exerted a more significant role in structuring the bacterial spatial distributions than soil environmental factors. In addition, the quantitative assessment of ecological processes revealed that dispersal limitation, variable selection, and homogeneous selection exerted important roles in bacterial community assembly processes. Moreover, such scenarios experienced obvious seasonal variations, due to the seasonal differences in river discharges. The links between taxonomic and functional diversities of soil bacteria demonstrated that taxonomic variability should not imply differences in bacterial community function, and there was high functional stability in the estuarine intertidal wetland microbial system. Our results also suggested that the possible transformation processes of C, N, and S were provided by highly active and adaptive bacterial metabolism channels, exhibiting notable significance to maintain the functional stability of bacteria. Our study provided important ecological insights about bacteria biogeographic patterns in China's estuarine intertidal soils and their driving factors. Moreover, the integration of bacterial taxonomic and functional data provided an important implication for the ecological functions of bacteria in estuarine intertidal wetlands.