Water and surface sediment samples from the PR to the NSCS region (18.501^° to 23.791^°N, and 112.164^° to 115.501^° E) were collected along a salinity gradient during three cruises in June 2011, January 2012, and October 2012 (Figure 1; Supplementary Table S1). Supplementary Table S1 provides additional information on eight sediment sample data from one of our previous studies (Xie et al., 2014). The uncontaminated seawater was taken from different depths and was sequentially filtered by glass-fiber filters (0.7 μm) to obtain microorganisms, stored on board in liquid nitrogen, and transferred to a laboratory freezer at −80°C for further analysis. Water samples for the determination of inorganic nutrients were collected, filtered through a 0.45 μm cellulose acetate membrane and stored at −20°C until analysis. Immediately after waterbed sampling, a centrifuge tube of 50 mL was used to seal the sediments. Then these tubes were stored on board in liquid nitrogen and transferred to a laboratory freezer at −80°C for further analysis. Before sediment collection, in situ measurements of salinity, temperature, and water depth were conducted with a YSI instrument (model YSI 650 MDS, YSI Corporation, United States). Before pore water extraction, samples were stored under −20°C until centrifugating at 4,000 rpm for further analysis. The pore water treatment method was as described in our previous reports (Xie et al., 2014).

DNA was extracted using the FastDNA SPIN Kit for Soil (MP Biomedicals, OH, USA) according to the manufacturer’s instructions and preserved at −80°C until further processing. Twenty-eight water samples and 17 sediment samples were selected for pyrosequencing targeting archaeal 16S rRNA (8 sediment data have been published in our previous article (Xie et al., 2014); detailed in Supplementary Table S1), and 20 water samples and 11 sediment samples were selected for pyrosequencing targeting the archaeal amoA gene. The Arch_344F/Arch_915R (Gantner et al., 2011), and Arch_amoAF/Arch_amoAR (Francis et al., 2005) primers were used for PCR of the archaeal 16S rRNA and AOA-amoA, respectively. Pyrosequencing was performed on the Roche GS FLX+ (454) system, raw data was extracted from 454 data formats by Mothur pipeline (version 1.29.2) (Schloss et al., 2009) and then analyzed using the QIIME2 standard pipeline (version 2020.6) (Bolyen et al., 2019). Specifically, the first step was to filter raw data to remove reads with low quality and ambiguity. The remaining high-quality sequences of 16S rRNA and amoA were clustered into operational taxonomic units (OTUs) with the command “qiime vsearch cluster-features-denovo” with 97 and 96% sequence identity thresholds, respectively. Taxonomic assignment for archaeal 16S rRNA representative sequences was obtained by the SILVA database (release 138; Quast et al., 2013). Phylogenetic clades of AOA were affiliated following past studies (Schleper et al., 2005; Alves et al., 2018).

Alpha and Beta diversity were calculated from the 16S rRNA and amoA datasets using uniform subsampled depths for 1732 and 1,283 sequences, respectively. Rarefaction curves and alpha diversity index for each sample were calculated using the diversity function within R package “vegan” (version 3.6.1; Oksanen et al., 2020). Sample clustering was employed at the archaeal OTU level and the unweighted pair-group method with arithmetic means (UPGMA) based on the Bray–Curtis similarity algorithm. Analysis of similarities (ANOSIM) was implemented using PAleontological Statistics (PAST, version 3.16) to test the differences among groups. Principal component analysis (PCA) was conducted to identify the dissimilar relationships among samples in this study based on the total archaeal community, and AOA subcommunity composition at the OTU level in the STAMP program (Parks et al., 2014). To determine the phylogenetic clades of amoA sequences, we obtained reference sequences and tree from a past study (Alves et al., 2018), performed the multiple sequences alignment with MAFFT v7 (Katoh and Standley, 2013), the phylogenetic placement of amoA sequences with EPA-ng (Barbera et al., 2019), and visualization with ITOL v5 (Letunic and Bork, 2021). Environmental heterogeneity was estimated by computing the average dissimilarity between sites (Huber et al., 2020) based on abiotic variables by using a Euclidean distance.

The neutral community model was used to evaluate the effects of stochastic processes on archaeal community assembly (Sloan et al., 2006), and nonlinear least-squares methods were used to determine the best fit between OTU occurrence frequencies and their relative abundance (Elzhov et al., 2016). Östman’s method was used to calculate the R² value, which indicates the goodness of fit of the model to datasets (Östman et al., 2010). As R² approaches 1, the community composition follows a completely stochastic process. R² can be ≤0 when the model does not match the community composition.

To weigh the relative importance between deterministic verse stochastic processes to the archaeal communities, we evaluated the deviation of each observed metric from the average of the null model (checkerboard score (C-score)) to test the overdispersion or underdispersion of archaeal communities (Stone and Roberts, 1990). The results were normalized so that they could be compared among groups using the standardized effect size (SES). SES was computed under the null model after converting the sequence table to a binary matrix with presence (1) and absence (0) (Gotelli and McCabe, 2002). The SES for the C-score was calculated by dividing the difference between the observed and the mean of the stimulated index by the standard deviation of the stimulated index (Crump et al., 2009). Overdispersion or underdispersion is defined as an SES value that is greater or lower than the expected null value. The strength of the influence of deterministic processes on communities is understood as the magnitude of SES (Swenson, 2014). The dominance of stochastic processes is indicated by SES values between −2 and 2, whereas SES values less than −2 or more than 2 show that deterministic processes (i.e., aggregation and segregation) are more important in community assembly than stochastic processes (Stone and Roberts, 1990; Mo et al., 2021). We calculated the C-score based on the sequential swap randomization algorithm and 30,000 simulations with a R-package named “EcoSimR” (R Core Team, 2013). In addition, to assess how stochastic and deterministic processes affect archaeal community assembly, build-in function “niche.width” of R package “spaa” was used to calculate Levins’ niche breadth (B) index as previous studies (Oksanen et al., 2013). OTUs with higher B values indicated a wider range of habitat niche breadth. Community-level B value (Bcom) was obtained by averaging the B values of all taxa occurring in this community (Wu et al., 2018; Jiao et al., 2020). At the community level, archaeal communities with wider niche breadths should exhibit greater metabolic flexibility than those with narrower niche breadths (Pandit et al., 2009; Wu et al., 2018; Jiao et al., 2020).

The co-occurrence networks were constructed using the Molecular Ecological Network Analysis Pipeline (MENA). Only OTUs that were found in more than two samples and with a proportion above 0.2% were retained to reduce complexity. Spearman’s correlation coefficient between OTUs was calculated, and relationships were classified as valid when the correlation coefficient |r| was greater than 0.6 with a p value less than 0.05. The network was visualized in Gephi version 0.9.2 (Bastian et al., 2009). Moreover, 1,000 Erdös-Réyni random networks, with the same number of nodes and edges as the real network, were calculated using the R package “igraph,” where each edge is assigned the same probability to each node (Erdos and Rényi, 1960). To further describe the topological parameters, the clustering coefficient, average path length, and modularity of both real and random networks were calculated and compared. Finally, we removed nodes in the static network with a random repetitive process to observe the changes in the natural connectivity of nodes to assess network robustness (Peng and Wu, 2016).

As indicated by the rarefaction curves for total archaeal communities and AOA subcommunities, all samples reached a plateau at the sequencing depth applied (Supplementary Figure S1). The number of reads were sufficient for a robust statistical analysis. According to the clustering results, all water samples were divided into two categories: the low-salinity-water group (salinity <14.58 ‰, nitrate ranged from 67.47 to 494.4 μM, silicate ranged from 1.13 to 157.64 μM, and phosphate ranged from 4.155 to 26.885 μM) and the high-salinity-water group (salinity >21.08 ‰, nitrate ranged from 0.139 to 149.223 μM, silicate ranged from 0.22 to 65.031 μM, and phosphate ranged from 0.016 to 6.481 μM; Supplementary Figure S2A, Supplementary Table S1). Similarly, the sediment samples were categorized into the low-salinity-sediment group and the high-salinity-sediment group (Supplementary Figure S2B). The dramatic environmental heterogeneity may account for the differences in archaea community structure between groups.

The composition of the archaeal communities at the phylum level is shown in Figure 2A. Euryarchaeota was the major dominant phylum in almost all low-salinity-water and low-salinity-sediment groups samples. The following were Bathyarchaeota and Woesearchaeota. Thaumarchaeota were also significant components in these two groups of samples. In the high-salinity-water group samples, Euryarchaeota and Thaumarchaeota were the predominant archaeal phyla. Bathyarchaeota and Woesearchaeota also account for a certain fraction. Compared with the other three groups, the proportion of Lokiarchaeota and Hydrothermarchaeota increased in the high-salinity-sediment group samples. Euryarchaeota, Thaumarchaeota, Bathyarchaeota, and Woesearchaeota were also dominant components in this group. The community composition of the low-salinity-water and low-salinity-sediment group samples was similar at the genus level (Figure 2B). The dominant Archaea genera were Bathyarchaeia, Woesearchaeales, Thermoprofundales, Methanosaeta, Methanoregula, Methanobacterium, Candidatus Methanoperedens, Methanolinea, and Methanosarcina. Marine Group II was dominant in most high-salinity-water group samples, except for samples PR120103_2_14 m, PR110615A_17 m, PR110617B_2 m, PR110617B_9 m, SCS_0_20 m, and SCS_1_3 m. Bathyarchaeia was dominant in sample PR120103_2_14 m, and Candidatus Nitrosopumilus was dominant in the other five samples. Marine Group III had a high proportion in some high-salinity-water group samples. Bathyarchaeia, Candidatus Nitrosopumilus, Woesearchaeales, Thermoprofundales, Nitrosopumilaceae, Lokiarchaeia, and Hydrothermarchaeales were dominant genera in the high-salinity-sediment group samples.

After filtering out the AOA OTUs of less than 1% in all samples, the others were mapped to a specific evolutionary tree (Alves et al., 2018) to determine their phylogenetic clades (Supplementary Figure S3, Supplementary Table S2). A total of 121 AOA OTUs were assigned to 20 subgroups, and sediments harbored greater clade diversity than water (Supplementary Figure S4, Supplementary Table S2). NP-α-1 was the dominant amoA type in almost all low-salinity-water group samples, except for sample PR110619C_1.5 m, whose community composition was more similar to that of the low-salinity-sediment group samples. In the low-salinity-sediment group samples, NS − δ, NP-η, and NP-γ-2.2 were dominant subgroups, and NT-α had a high abundance in sample PR120112_21S. Sample PR110617B_9 m had a different community composition from the other samples in the high-salinity-water group. In PR110617B_9 m, NP-γ-2.1 had the highest abundance, followed by NS − δ. NP-α-1 and NP-ε were dominant subgroups in the other samples. In the high-salinity-sediment group samples, NP-α-1, NP-α-2.2.4, NP-δ, NP-α-2.1, NP-α-2.2.3, and NP-θ were the predominant subgroups.

Principal component analysis showed that the community exhibited obvious group segregation characteristics. This is in line with the changes of the relative abundance of dominant archaeal clades in each group (Supplementary Figure S5A). This pattern is supported by ANOSIM analysis, where large differences between clusters were observed (Supplementary Table S3). It is worth noting that there was no significant separation between the low-salinity-water and low-salinity-sediment group samples according to the ANOSIM analysis, which was consistent with the PCA. The AOA subcommunities differed significantly among the four groups (Supplementary Table S4). The AOA composition differed significantly between the low-salinity-water and low-salinity-sediment groups (Supplementary Table S4), indicating the specificity of AOA between water and sediments. The differences in AOA composition between the high-salinity-water and low-salinity-water groups were small (Supplementary Table S4), indicating a high tolerance of AOA to salinity.

For each sample, the diversity indices including the Shannon index and the Observed-species number, were shown in Supplementary Table S1. The Shannon and observed-species diversity indices of the high-salinity-water group samples were significantly lower than those of the other three groups, while these indices were not significantly different among the three groups (Supplementary Figure S6A). For AOA subcommunities, the Shannon and observed-species indices of the high-salinity-water and low-salinity-water groups were lower than those of the high-salinity-sediment and low-salinity-sediment groups (Supplementary Figure S6B).

The archaeal community immigration rate was estimated to be 0.1796, 0.4475, 0.3886, and 0.5878 in the high-salinity-water, low-salinity-water, high-salinity-sediment, and low-salinity-sediment groups, respectively (Figure 3). In addition, the archaeal community immigration rates across all water and all sediment were smaller than the corresponding high-salinity and low-salinity groups (Figure 3), indicating that the species dispersal within groups was stronger than that among groups.

The high-salinity-water, low-salinity-water, high-salinity-sediment, and low-salinity-sediment groups contained 11, 75, 60, and 73 OTUs, respectively, which occurred in over 70% of their samples (Figure 3). However, only 9 and 9 OTUs occurred in more than 70% of all-water and all-sediment samples (Figure 3), indicating that most archaeal species can only be found in specific groups and few generalists. To judge whether the high occurrence frequency of OTUs was due to differences in the number of samples between groups, we plotted bubble charts for OTUs that occurred in more than 70% of all-water and all-sediment samples. The results showed that only Ca. Nitrosopumilus (OTU7 and OTU37) were able to occur in both the high-salinity-water and low-salinity-water group samples (Supplementary Figure S7A), indicating the high tolerance of Ca. Nitrosopumilus to salinity. However, MGII and MGIII only appeared in the high-salinity-water group, indicating that they were only adapted to the high-salinity-water environment, and their high occurrence frequency was due to the larger number of samples in the high-salinity-water group. In all sediment samples, MBGD (OTU101) and Bathyarchaeia (OTU448 and OTU72) were present in both the high-salinity-sediment and low-salinity-sediment samples (Supplementary Figure S7B). The high frequency of other OTUs was due to the larger number of samples in the high-salinity-sediment group. We randomly selected the same number of high-salinity and low-salinity samples to draw the neutral community model, and the results were consistent with the bubble charts (Supplementary Figure S8), indicating that Ca. Nitrosopumilus in water, and MBGD and Bathyarchaeia in sediment were archaeal taxa that occurred with high frequency across drastic salinity gradients from the PR to NSCS.

We constructed archaeal co-occurrence networks in water and sediments (Supplementary Figure S9), and analyzed the high-salinity-water, low-salinity-water, high-salinity-sediment, and low-salinity-sediment subnetworks (Figure 4, Supplementary Table S5). The topological properties of these subnetworks varied significantly with salinity. The proportion of negative links in the low-salinity environments was higher than that in the high-salinity environments (Figure 4, Supplementary Table S5), which may indicate that more archaeal taxa in the low-salinity environments tend to adopt a competitive survival style. The eight largest modules accounted for 64.91, 83.09, 55.44, and 91.64% of the high-salinity-water, low-salinity-water, high-salinity-sediment, and low-salinity-sediment subnetworks, respectively (Figure 4). And the network diameter and average path length of the low-salinity groups were lower than those of the high-salinity groups, while the density and average clustering coefficient were higher than those of the high-salinity groups (Supplementary Table S5), indicating that the archaeal community in the low-salinity groups had closer interactions than that in the high-salinity groups. Furthermore, power-law distributions rather than Poisson distributions were observed for all network degrees (Supplementary Table S5). Accordingly, the network presented a structure of scale-free and non-random distribution. The observed network parameters (i.e., average path length, average clustering coefficient, and modularity index) were greater than the corresponding Erdös-Réyni random networks (Supplementary Table S5), indicating a “small world” nature and modular structure.

In the water network, archaeal communities were dominated by taxa that preferred high salinity within modules III and VII, and by taxa that preferred low salinity within modules I, II, IV, V, VI, and VII (Supplementary Figures 9B,C). Furthermore, MGII, MGIII, and Candidatus Nitrosopumilus exhibited the highest degrees of centrality in the modules for high salinity levels. In the modules for low salinity levels, however, Woesearchaeales, Bathyarchaeia, Methanosaeta, and Methanobacterium exhibited the highest degrees of centrality (Supplementary Figure S9A, Supplementary Table S6). Similarly, in the sediment network, archaeal communities were dominated by taxa preferring high salinity in modules I, IV, V, and VII (Supplementary Figures S9E,F), with Bathyarchaeia, Woesearchaeales, and Nitrosopumilaceae exhibiting the highest degrees of centrality (Supplementary Figure S9D, Supplementary Table S7); while archaeal communities were dominated by taxa preferring low salinity in modules II, III, VI, and VII. Bathyarchaeia, Woesearchaeales, and Methanosaeta exhibited the highest degrees of centrality. Finally, we compared the network stability of different archaeal subnetworks. The low-salinity subnetworks had higher natural connectivity than the high-salinity subnetworks (Supplementary Figure S10), indicating that the low-salinity subnetworks had stronger robustness.

The standardized effect size (SES) was calculated based on OTU abundance to investigate the relative contributions of deterministic and stochastic processes to archaeal community assembly (Figure 5). The SES values in all four groups were greater than 2, the threshold for stochastic and deterministic processes, indicating that deterministic processes contributed more to archaeal community assembly than stochastic processes from the PR to the NSCS (Figure 5). Consistent with this, the relationship between OTU occurrence frequency and relative abundance can only be characterized by the neutral community model to a small extent (Figure 3), indicating that stochastic processes accounted for a smaller proportion of archaeal community assembly. For the high-salinity-sediment, low-salinity-sediment, high-salinity-water, and low-salinity-water groups, stochastic processes explained 39.8, 23.8, 42.5, and 14.3% of the community variation, respectively, showing a trend of decreasing gradually with the decrease in salinity. This indicated that stochastic processes contributed more in the high-salinity than low-salinity environments from the PR to the NSCS. Furthermore, stochastic processes showed lower contributions to archaeal community assembly in all-water and all-sediment samples (Figure 3). Additionally, in the neutral community model with the same number of high-salinity and low-salinity samples, the results suggested that stochastic processes contributed less to the assembly of archaeal community. In planktonic archaeal communities, the neutral community model (NCM) cannot even describe the archaeal community composition (Supplementary Figure S8). All archaeal communities had wider niche breadths at the high salinity than at the low salinity levels (Supplementary Table S8).

This study revealed detailed archaeal community composition across four environments from the PR to the NSCS by high-throughput sequencing analysis. Although the archaeal communities were sampled across time and seasons, ANOSIM analysis revealed that archaeal communities differed significantly among habitats in the current study (Supplementary Table S3). Because the archaeal community samples in the four groups were segregated in PCA (Supplementary Figure S5A), we surmised that seasonal differences within groups were unlikely to be more influential than differences among groups in this study. This result confirmed previous findings that benthic communities of bacterial (Böer et al., 2009) and microbial eukaryotes (Gong et al., 2015) varied with sediment depth or area rather than the season.

Archaeal community structure showed a distinct pattern of heterogeneity (Figure 2). Previous reports have indicated that the structure of microbial communities varies spatially from estuaries to the open ocean (Herlemann et al., 2011; Liu et al., 2014). It can also be seen here that the composition of planktonic and benthic archaeal communities differed considerably from the PR to the NSCS. We observed a high abundance of methanogenic archaea in the low-salinity environments (Figure 2B). Salinity is reportedly the most important abiotic factor governing the distribution of methanogenic communities worldwide and determines which lineages may survive to a large extent (Wen et al., 2017). Increasing salinity is reported to restrain hydrogenotrophic methanogens while promoting acetoclastic methanogenesis (Liu et al., 2016). Additionally, the activity of methanogens that utilize acetate or H₂ was highly impacted by pH. Low pH, for example, makes acetoclastic methanogenesis difficult because it lowers acetate dissociation (Megonigal et al., 2004; Kotsyurbenko et al., 2007). MGII and MGIII were typically found in seawater and were almost absent in freshwater environments, suggesting that suitable salinity is necessary for their survival. Bathyarchaeota were widespread in the low-salinity environments and the high-salinity-sediment samples, and were the dominant taxa in these environments (Figure 2). As one of the most important and active microbial groups in freshwater and marine sediments, the relative abundance of Bathyarchaeota accounted for 36 ± 22% of the total archaea abundance (Fillol et al., 2016). Previous studies have shown that different Bathyarchaeota subgroups are indicators for different environments: subgroups 1 and 8 in marine environments, and subgroups 5 and 11 in freshwater environments (Fillol et al., 2016; Wang et al., 2020).

Thaumarchaeota was suggested to be a salinity-tolerant species (Xie et al., 2014). Based on amoA genotypes suggested by Alves et al. (2018), we determined the AOA subcommunity composition from the PR to the NSCS in the current study (Supplementary Figure S4). NS-δ is abundant in the low-salinity sediment samples. Previous studies have indicated that clade NS-δ accounts for a significant proportion of AOA communities in salt lakes and freshwater habitats worldwide (Alves et al., 2018). Considering its global environmental distribution, NS-δ represents either versatile microorganisms or ubiquitous soil microorganisms that are frequently transferred to adjacent rivers. A systematic analysis of AOA on a global scale found that clades NP-α, NP-ε, NP-θ, and NP-ζ occur primarily in marine environments, NP-γ occurs primarily in freshwater and estuarine-coastal environments, and NP-η occurs more frequently in freshwater and soils-sediments (Alves et al., 2018), consistent with the results in the low-salinity environments in the current study (Supplementary Figure S4). The composition of AOA clades in seawater and sediment is almost different: most AOA in seawater belong to clades NP-ε-2 and NP-α-2.2.2.1, while sediment harbor greater clade diversity, including NP-δ, NP-γ, NP-θ, and NP-α (Durbin and Teske, 2010; Alves et al., 2018; Zou et al., 2020a), which is consistent with the results in the high-salinity environments presented here (Supplementary Figure S4). Overall, our results showed that the change in environmental conditions from the PR to the NSCS controlled the distribution of archaeal communities.

The archaeal community immigration rates in the low-salinity environments were significantly higher than those in the high-salinity environments (Figure 3), suggesting the higher dispersal ability of most archaeal taxa in the upper PR than in the NSCS. These results may be attributed to the faster flow in the upper PRE, which facilitated the dispersal of archaeal communities in the water and surface sediment. Additionally, the species dispersal of archaea among groups was smaller than that within groups, indicating that the majority of archaea had difficulty surviving across drastic salinity gradients. This is consistent with the results of the neutral community model, in which the majority of archaea only occurred in specific groups. Only Ca. Nitrosopumilus (OTU7 and OTU37) in water, and Bathyarchaeia (OTU448 and OTU72) in sediments were able to occur in both low-salinity and high-salinity environments. Previous studies have shown that salinity is an important environmental factor limiting the distribution of archaea (Webster et al., 2015), Thaumarchaeota was suggested to be a salinity-tolerant species (Xie et al., 2014), and Bathyarchaeia is widely present in various environments with high phylogenetic diversity and abundance (Fillol et al., 2016; Wang et al., 2020).

Rather than growing in isolation, microbes in natural ecosystems prioritize forming complex interaction networks (Banerjee et al., 2019). Both the composition and dynamics of microbial communities are strongly influenced by microbial interactions, and co-occurrence networks can reveal interactions between species, such as competition and cooperation (Berry and Widder, 2014; Lima-Mendez et al., 2015; Wei et al., 2016). The mostly positive correlations in the high-salinity subnetworks identified in this study (Figure 4, Supplementary Table S5) suggest that cooperation is more prevalent than competition in archaeal communities (Berry and Widder, 2014). While this phenomenon is common in natural ecosystems, it is perhaps not surprising since many microbes rely heavily on cross-feeding, co-aggregation, co-colonization, or niche overlap and construction (Faust and Raes, 2012; Faust et al., 2015; Shi et al., 2019). A previous study that focused on the co-occurrence networks under the influence of graphene oxide and different temperatures found that when facing the enhancement of heterogeneity selection, the strategy of the bacterial community was to enhance the positive correlation and shared niche, while the strategy of the archaeal community was to enhance the negative correlation and competition (Liao et al., 2022), which is consistent with more negative correlations in low-salinity environments with higher environmental heterogeneity (Figure 4).

Our results showed that network modules corresponded well to groups (Supplementary Figure S9). This suggested that the modular structure or properties of planktonic and benthic archaeal communities were sensitive to changes in environmental conditions. Studies have shown that salinity and other environmental factors (i.e., pH, temperature, and nutrients) could change bacterial and eukaryotic plankton co-occurrence networks (Ji et al., 2019; Liu et al., 2019). Environmental heterogeneity can induce microbial modularity (Röttjers and Faust, 2018), explaining why these modules predominated across different groups. Modularity could be used to indicate both niche differentiation and competitive/cooperative relationships, leading to non-random network structures that ultimately increase the ecological network complexity (Olesen et al., 2007).

For low-salinity subnetworks, the average clustering coefficient in the network topology parameters was higher, implying the higher complexity of these subnetworks. Highly complex networks tend to be more stable due to network buffering (Landi et al., 2018), so the low-salinity subnetworks may be more stable, while the high-salinity subnetworks may be less stable, which was also supported by Supplementary Figure S10. This pattern may be closely related to archaeal community diversity and richness, as higher diversity and richness may lead to more complex network structures (Chen and Wen, 2021). Changes in environmental heterogeneity may be the most important factor responsible for different patterns of microbial interactions among groups, as the complexity and connectivity of networks often had positive correlations with environmental heterogeneity (Mougi and Kondoh, 2012). Along the upper reaches of the Pearl River to the estuary (low-salinity groups), fresh water with low salinity and high-concentration nutrients and organic matter characteristics meets seawater with high salinity and low-concentration nutrients and organic matter, to form a highly environmental heterogeneity. And this heterogeneity decreased along the estuary to the shelf (high-salinity groups) (Supplementary Figure S12). Microbial communities were more likely to aggregate and form more connected networks in more heterogeneous environments (Morriën et al., 2017; Li et al., 2019). Conversely, in more homogeneous environments, microbial communities tended to occupy similar biological niches, resulting in fewer interactions (Cherif and Loreau, 2007; Freilich et al., 2011; Lu et al., 2022).

Since nodes with a higher degree of centrality are more likely to be important to preserve taxa coexistence, MGII, MGIII, and Candidatus Nitrosopumilus may play a key role in maintaining taxa coexistence in the high-salinity-water subnetwork, while Woesearchaeales, Bathyarchaeia, Methanosaeta, and Methanobacterium are more important in maintaining coexistence in the low-salinity-water subnetwork (Supplementary Table S6). Similarly, in the high-salinity-sediment subnetwork, Bathyarchaeia, Woesearchaeales, and Nitrosopumilaceae may play a key role in maintaining taxa coexistence, while Bathyarchaeia, Woesearchaeales, and Methanosaeta are more important in maintaining coexistence in the low-salinity-sediment subnetwork (Supplementary Table S7).

The C-score analysis supported that deterministic processes played a relatively more important role than stochastic processes, whether in the water or sediments from the PR to the NSCS (Figure 5). Previous studies have indicated that geographic scales and environmental gradients largely explain the balance between deterministic and stochastic processes (Hanson et al., 2012). In general, at relatively large spatial scales (~12,000 km), deterministic (environmental) factors have been found to exert a greater impact on benthic microbial communities than stochastic (spatial) factors (Zhang et al., 2018), while at smaller spatial scales (~20 to 1,500 km), the contribution of stochastic processes was greater (Chen et al., 2017; Liu et al., 2020; Lu et al., 2022). However, in the current study with a small spatial scale (~600 km), drastic changes in environmental factors were found in both low-salinity and high-salinity environments, resulting in a greater contribution of deterministic processes than stochastic processes at these sites, which may imply the importance of environmental factors for the construction of archaeal communities from the PR to the NSCS (Xie et al., 2014, 2018; Zou et al., 2020b). Consistent with our results, a previous study found that deterministic processes had a stronger effect on soil bacterial communities in a 1,092 km range, which was mainly driven by a large pH gradient (Shi et al., 2018). In a study about the prokaryote community assembly from three parallel ~100 km transition sections in the lower Pearl River Estuary, a greater role of stochasticity to the benthic prokaryote community assembly was observed (Lu et al., 2022). However, it should be noted that the relative contribution of stochastic processes to the microbial community assembly was only 54.02%, which may be because the objects of this study were prokaryotes, and the smaller spatial scale led to less environmental heterogeneity compared with our study. Additionally, the C-score revealed that the SES values increased with decreasing salinity, implying the enhanced importance of deterministic processes in archaeal communities, showing the relative contribution of deterministic processes to low-salinity samples was higher than that of high-salinity samples, possibly due to more heterogeneous environment attributed to the steeper environmental gradients in the low-salinity environments (Supplementary Figure S12), as well as the more spatially inconsistent and complex hydrologic movement (Stegen et al., 2015).

The high-salinity environments had wider archaeal community niche breadths, while the low-salinity environments had narrower community niche breadths (Supplementary Table S8), indicating that the high-salinity environments had more generalists that can adapt to a wide range of environmental niches, while the low-salinity environments had more specialists. This result may be attributed to the higher nutrient and organic matter content in the low-salinity environments, which promoted the specialization of the archaeal niche, resulting in narrower habitat niche breadths. Conversely, in the high-salinity environments, the limited availability of resources induced the generalization of the archaeal niche, leading to wider habitat niche breadths (Pandit et al., 2009). A balance between determinism and stochasticity is determined by the resistance of microorganisms to environmental changes, along with spatial scales and gradients in the environment (Liu et al., 2020). Microorganisms that adapt to high-salinity environments have evolved “salt in” and “salt out” strategies to regulate the cytoplasm to osmotic pressure (Oren, 2011). In high-salinity environments with less environmental heterogeneity or fewer competitive interactions among environmental generalists, the relative contribution of stochastic processes may increase (Jiao et al., 2020). Seasonal changes in pollutant concentrations brought about by terrigenous input may put archaea under great selection pressure from the upper Pearl River to the estuary (Yin, 2002; Chen et al., 2009; Lu et al., 2018). The results of the neutral community models indicated that the stochastic processes had a low explanation for the archaeal community composition from the PR to the NSCS (Figure 3) or even could not explain it (Supplementary Figure S8). The difference in the explanation of the archaeal community between the two neutral community models may be due to more samples in the high-salinity groups, which increased the proportion of stochastic processes. Additionally, the archaeal communities exhibited narrower niches in the low-salinity environments than that in the high-salinity environments (Supplementary Table S8), implying that the assembly of archaeal communities under low-salinity conditions was more affected by deterministic processes, possibly because deterministic processes generally have greater effects on habitat specialists with narrow niche breadths than on generalists with wide niche breadths (Pandit et al., 2009; Wu et al., 2018).

Although the C-score analysis and NCM well described the archaeal community assembly patterns in the current study, it is difficult to infer the specific effects of environmental and spatial variables on archaeal community structure (Chisholm and Pacala, 2010). In addition, the impact of unmeasured environmental factors and species interactions on archaeal community structure needs to be considered (Vellend et al., 2014; Zhou and Ning, 2017). Because estuarine ecosystems are overly dynamic, some key environmental factors may change in a stochastic manner, and our sampling involved different seasons, which may mask the main ecological patterns. To fully understand the archaeal community assembly mechanisms in subtropical estuaries, we need further experimental programs that take space and time into consideration.