In June 2022, 22 stations were selected from inshore to ocean near the northern Hangzhou Bay ( Figure 1 ). Surface water samples at 0.5 m depth were collected at each site. Approximately 500 ml of water sample from each site was filtered onto a 0.2 μm polycarbonate membrane (Millipore, USA), and immediately stored in a tank filled with liquid nitrogen before DNA extraction. We imported the global positioning system (GPS) coordinates that were recorded at every sampling point into the NOAA website (http://www.nhc.noaa.gov/gccalc.shtml) for calculating the pairwise geographic distance between samples.

Table 1 lists the water chemical parameters, including salinity, pH, exchangeable nitrite (NO₂ ^–N), nitrate (NO₃ ^–N), ammonium (NH₄ ⁺-N), Fe²⁺ and Fe³⁺ concentration, total organic carbon (TOC), total phosphorus (TP), total nitrogen (TN) and Sulfide concentration. A portable water quality analyzer (HQ 40d, HACH, USA) assisted in measuring the pH and salinity. The continuous-flow injection analysis was used to measure the NH₄ ⁺, NO₃ ^-, and NO₂ ^- concentration (Skalar Analytical SAN+ +, Netherlands) (Hu et al., 2017; Lin and Lin, 2022). Measurement of Fe²⁺ and Fe³⁺ concentration followed the previous literature (Zheng et al., 2019) together with their references (Roden and Lovley, 1993). The methylene blue spectrophotometric method assisted in analyzing the sample sulfide (Cline, 1969). Vario EL CN Elemental Analyzer (Elementar, Germany) assisted in measuring the TOC and TN (Zhang et al., 2015). And the ascorbic acid-molybdate blue method for assaying TP (Hou et al., 2013). All of these physiochemical parameters were in triplicate.

Plankton-archaea microorganisms were collected by filtering about 500 mL of water samples with a 0.2µm pore size polycarbonate membrane (Millipore, USA). Power Soil DNA Isolation Kit (MO BIO, USA) assisted in extracting community DNA from polycarbonate membranes in line with the protocol of the manufacturer. The primer set 524F/958R was used for archaeal 16S rRNA gene amplification (Shen et al., 2017), modifying the forward primer to contain a unique 8 nt barcode at the 5’ end. Gene amplification was conducted in a 20-μL reaction system containing 4 μL of FastPfu Buffer (5×), 2 μL of dNTP mix (2.5 mM), 0.8 μL of each primer (5 μM), 0.4 μL of Fastpfu polymerase, 10 ng of template DNA, and 0.2 μL of BSA. The PCR parameters were 96°C for 5 min, followed by 35 cycles of 96°C for 35 s, 56°C for 30 s, and 72°C for 40 s, with a final extension at 72°C for 5 min. We employed electrophoresis to identify the PCR products, which were located in the distinct band range of 400-450 base pairs (bp), within a 2% (w/v) agarose gel. These PCR products were then combined in equal proportions and subjected to purification using the GeneJET Gel Extraction Kit from Thermo Scientific (Jiao et al., 2016). Subsequently, four replicate PCR products derived from a single DNA sample were consolidated and purified using a DNA gel extraction kit from Axygen, USA, following the manufacturer’s guidelines. The concentration and purity of the extracted DNA were ascertained utilizing a NanoDrop spectrophotometer. Subsequently, we mixed three parallel PCR products in each subsample to sequence on an Illumina platform.

FLASH (V1.2.7, http://ccb.jhu.edu/software/FLASH/) assisted in merging the paired-end reads, and the quality of reads was filtered based on the literature (Caporaso et al., 2011). USEARCH software was applied to the detection and removal of acquired sequences following the UCHIME algorithm (Edgar et al., 2011). These sequences had their barcodes during the assignation to each sample. UPARSE software package was applied to the sequence analysis by the UPARSE-OTU algorithm. We assigned sequences with a similarity of 97% to the same OTU. The RDP classifier assisted in assigning the sequences representing each OTU to the taxonomic groups (Caporaso et al., 2011).

Null model analysis was carried out using the framework described by Stegen et al. (2013), to classify community pairs into underlying drivers of deterministic process and stochastic process based on the β-nearest taxon index (βNTI). Deterministic processes including heterogeneous selection (βNTI<−2) and homogeneous selection (βNTI> 2), and the stochastic processes including homogenizing dispersal (mass effect), dispersal limitation, and drift were distinguished by the Raup-Crick metric (RCbray) as described by Liu et al. (2020). The null model expectation was generated using 999 randomizations. Network assisted in exploring the co-occurrence pattern exhibited by microbial taxa, selecting OTUs with top 500 relative abundances. A Spearman’s correlation between 2 OTUs presented statistical robustness with a coefficient over 0.6 and a P-value less than 0.01 (Barberán et al., 2012; Jiao et al., 2016). The interactive platform Gephi assisted in visualizing the network (Bastian et al., 2009).

QIIME (http://qiime.org/index.html) was employed for calculating the Alpha and Beta diversity, and the minimum sequence number used for normalizing the sequencing depth difference. The potential for the cycling of C, N, and P elements within microbial communities was forecasted based on the composition of the 16S rRNA gene library, using the FAPROTAX software version 1.2.1. Differences in environmental variations, community diversity, abundances, and functional profiles among microbial groups were compared via the Kruskal-Wallis H test and One-way analysis of variance (ANOVA), respectively. We confirmed the relation of geographical coordinates or the environmental factor (the Euclid distance based on environmental factors) with the sample ordination on the taxonomic (Bray-Curtis). The Mantel test assisted in checking the effect of environmental and spatial dissimilarity on taxonomy. Nonmetric multidimensional scaling (NMDS) was used to evaluate the structure and function differences of microbial communities (Caporaso et al., 2010; Legendre and Legendre, 1998).

Two-way analysis assisted in comparing the spatial variations of environmental variables. Based on the values of lengths of gradient, we applied to the redundancy analysis (RDA) examining the relationships between microbial taxonomic and functional structure and environmental indices in Canoco (version 4.5) software, respectively (Danovaro and Gambi, 2002). The linear trend together with the principle coordinates of neighbor matrices (PCNM) procedure assisted in deriving the spatial variable from the geographic coordinates (Griffith and Peresneto, 2006), thereby minimizing the association of spatial distance with environmental factors, and capturing all spatial scales detected in the dataset. Then adjusted R² was used in the analysis of variation partitioning, to determine the different effects of spatial factor and environmental factors on the microbial community variation, as well as the combined effect of the two factors.

Except as described above, the rest of the data analysis processing is done in R environment (http://www.r-project.org). The raw reads of high-throughput sequences have been deposited in the NCBI SRA database.

There was an obvious salinity gradient along the transect of samples ( Table 1 ). The salinity distribution was considered to divide these study sites into high-salinity regions (22.60–35.50 ppt, S16 to S22), mid-salinity regions (11.30–19.70 ppt, S10 to S15), and low-salinity regions (1.60-8.90 ppt, S1 to S9). All samples were alkalescent and pH ranged from 8.00 to 8.73. Water nutrients NH₄ ⁺ showed significant variations among three salinity gradients (two-way ANOVA, P = 0.02). The values of Fe²⁺ and Fe³⁺ were highest at low-salinity regions (two-way ANOVA, P < 0.05). Concentrations of water sulfide also showed substantial variability across the salinity zones (two-way ANOVA, P < 0.001), with the highest levels observed in the high-salinity areas. Nevertheless, the variations in other water environmental parameters did not reach statistical significance, implying that the observed changes in water geochemical characteristics are likely due to specific biogeochemical processes rather than water salinity shifts.

The Miseq sequencing platform assisted in identifying the structure of the archaeal community regarding all water samples. We obtained a data set containing 1551070 quality sequences from the 22 samples. And the total OTU number was 7456 defined by 97% sequence similarity. Table 2 lists the archaeal community abundance and diversity in sediments along the salinity gradient from inshore to ocean near northern Hangzhou Bay. All of the coverage of the archaeal sequencing library reached 100%, indicating that the primers used in the study effectively reflect the diversity of archaeal populations in the samples. The Shannon indices were in the range of 0.130 to 3.714, and the highest values were seen at mid-salinity sites S14. Chao 1 and ACE indices ranged from 253.240 to 1089.500 and 262.860 to 1683.239, respectively. Moreover, the highest values of these two richness indices appear at mid-salinity sites S10. According to the results of the alpha diversities index, although the difference was not significant, archaeal communities at high-salinity regions presented relatively higher diversity, with the average values of the Shannon index in high-, middle- and low-salinity regions were 1.740, 1.040, and 1.470, respectively. However, the richness was highest in high-salinity regions. The average values of Chao 1 and ACE indexes are the highest in the mid-salinity regions, which are 665.909 and 875.850 respectively (one-way ANOVA, P < 0.05).

Most sequences (99.84%) were included in the phyla of Woesearchaeota, Thaumarchaeota, Euryarchaeota, and Crenarchaeota (relative abundance > 2%), occupying 99.28% of the total sequences ( Figure 2A ). As revealed by NMDS analysis, the planktonic archaeal community distribution exhibited a distinct geographical specificity along the salinity gradient. According to the first two axes of NMDS, there were 3 groups (low-, middle- and high-salinity regions) of microbial assemblages along the salinity gradient ( Figure 2B ). An analysis of similarity made a further conformation on these patterns, ANOSIM result indicating the vital significance of sampling sites for determining the community composition (the mantel test, R = 0.792, P < 0.001). Based on all these results, archaeal composition presented a dynamic spatial pattern along the salinity gradient.

Figure 3 displays the impact of environmental factors on the structure of archaeal communities along various salinity gradients. The clusters based on the water environment factor in the right diagram showed a similar cluster to the microbial community structure in the left diagram, confirming the obvious association between them (the mantel test, R = 0.792, P < 0.001) ( Figure 3A ). RDA analysis assisted in exploring the effect of environmental factors on archaeal communities ( Figure 3B ). The first 2 RDA axes accounted for 43.96% of the cumulative variation regarding the relation between them. In which, microbial community distribution and structure were affected by salinity (P = 0.001), NH₄ ⁺ (P = 0.003), TP concentration (P = 0.008), pH (P = 0.009), TN (P = 0.011) as well as spatial factors (PCNM 1 and PCNM 2, P < 0.05). The variance partitioning analysis (VPA) was further used to quantify the effect of environmental and spatial factors on the variation of the archaeal community. The effect of the two factors accounted for 38.66% of the community variation. The pure effect of environmental factors explained 12.11% of the variation. While the PCNM purely explained 6.62% of the variation. Also, the combined effect of PCNM and environmental variables accounted for 19.93% of the variation, demonstrating the mutual dependence of the two factors.

To explore community assembly mechanisms under the observed geographic pattern, the relative roles of niche and neutral processes in community assembly were analyzed. Across all samples, stochastic processes (-2< βNTI <2) explained a slightly higher proportion than deterministic processes (2< |βNTI|) of the archaeal community variation, and the balance of these ecological processes varied greatly across the salinity gradient ( Figure 4A ). Our result shows that the deterministic processes explained a lower proportion (16.67%) of the archaeal community variation in middle-salinity regions than in the low- and high-salinity regions (34.57% and 36.74%, respectively). In addition, the effect of homogenous selection (βNTI> 2) on community assembly in the middle-salinity areas was 0%. For the stochastic processes, the result showed that drift was the most important process, accounting for 33.33%, 50.00%, and 34.69% of the community variation in high-, middle- and low-salinity regions, respectively. Then was followed by dispersal limitation, which respectively explained 9.88%, 27.78%, and 24.49% of the community variation in high-, middle- and low-salinity regions ( Figure 4B ).

The co-occurrence network results indicate that compared to areas with low and medium salinity, the archaeal co-occurrence networks in high salinity areas have significantly higher avgK and avgCC ( Figure 5A ). The results of the node-level topological properties of the co-occurrence network show that archaeal communities have a significantly higher node degree, node CC, and degree centrality indices, as well as a notably lower betweenness centrality index in high salinity areas ( Figure 5B ), suggesting that the archaeal community in high salinity areas has more complex interaction relationships. Furthermore, our results show that from areas with low salinity to those with high salinity, there is a clear gradient change in the network-level properties (such as avgK and avgCC), node-level properties (such as node degree, node CC, and betweenness centrality), and modularity of the archaeal symbiotic network, implying a key regulatory role of water salinity on the co-occur relationships of archaeal communities.

The prediction results from FAPROTAX indicate that the main metabolic mode of the archaeal community in the study area is aerobic ammonia oxidation, followed by nitrification, methanogenesis, and hydrogenotrophic methanogenesis. Both aerobic ammonia oxidation and nitrification processes are predicted to be present in all samples, suggesting a key impact of archaeal communities on the biogeochemical nitrogen cycle ( Figure 6A ). Additionally, microbial groups involved in methanol oxidation, anammox, nitrogen fixation, aromatic compound degradation, and iron respiration have been detected ( Figure 6A ). Beyond these, enzymatic processes such as chitinolysis, cellulolysis, xylanolysis, and ureolysis are also present ( Figure 6A ; Supplementary Table S1 ). This indicates that the archaeal communities in the study area are involved in a variety of biogeochemical cycles or metabolic processes of organic and inorganic substances. Among the nutritional types, chemoheterotrophy is the dominant type of microorganism, with high abundance in all three salinity zones and peaking in the medium salinity area; photoautotrophs and photoheterotrophs are also distributed, mainly in the medium salinity area ( Figure 6A ; Supplementary Table S1 ). Furthermore, from low to high salinity areas, functions such as aerobic ammonia oxidation and nitrification significantly increase, while functions like Methylotrophy, Methanogenesis by reduction of methyl compounds with H₂, and Aerobic chemoheterotrophy show a significant decrease ( Figure 6B ).

The NMDS results indicate that the functional characteristics of the archaeal community do not show significant variation across the salinity gradient. (ANOSIM, R = 0.012, P = 0.241) ( Figure 6C ), CCA was used to explore environmental and biological effects on the observed functional pattern ( Figure 6D ). The CCA result showed that the environmental factors including salinity and pH were the most powerful predictors separating functional patterns along the salinity gradients (P < 0.05) ( Figure 6E ). In addition, according to fitting analyses, significant correlations were confirmed between salinity and functional β diversity (P < 0.001, R² = 0.34), as well as microbial richness (Chao 1 index) and functional β diversity (P = 0.029, R² = 0.28).

The archaeal communities were under the dominance of the phyla of Woesearchaeota, Thaumarchaeota, Euryarchaeota, and Crenarchaeota. This pattern presents a similarity to the patterns of other marine waters in China. It has been reported the temporal distribution of archaeal communities in the Changjiang Estuary and the adjacent East China Sea was examined, revealing that Euryarchaeota and Crenarchaeota were the predominant groups, with salinity identified as a key factor influencing their distribution (Liu et al., 2011). Our findings align with their observation of salinity’s impact; However, we also detected a more diverse range of archaeal phyla, suggesting a broader response of archaeal communities to the complex interplay of environmental factors in our study area.

It is noteworthy that the high abundance of Woesearchaeota in our study is a striking finding that contrasts with some previous reports. However, recent research has begun to uncover the ecological preferences and adaptive strategies of this archaeal group, which may explain their success in specific environments (Liu et al., 2021). In the context of our study, the high salinity and nutrient-rich conditions of the coastal waters, influenced by the estuarine environment, may favor the proliferation of Woesearchaeota. Similar findings have been reported, such as in the Qiantang River estuaries where Woesearchaeota were found to be dominant, suggesting that the dynamic physicochemical characteristics of the estuarine mixing zone may provide a suitable habitat for this group (Liu et al., 2020). Another possible explanation is that Woesearchaeota widely distributed in terrestrial ecosystems (Wei et al., 2016). The runoff from rivers may carry a large number of Woesearchaeota cells/species into the estuarine water bodies, promoting the high enrichment of this species. Thaumarchaeota and Euryarchaeota dominate in marine water bodies far from the estuary, where the sources of chemical energy differ from those in estuarine water bodies, such as the lack of available organic matter, which may lead to the adaptation and dominance of specific microbial groups. Thaumarchaeota and Euryarchaeota can utilize dissolved organic carbon through chemosynthesis or obtain energy through ammonia oxidation, ensuring their survival and proliferation in the marine environment, thereby occupying a dominant position in marine water bodies far from the estuary. Furthermore, Thaumarchaeota plays an important role in the marine nitrogen cycle, capable of performing ammonia oxidation in the deep waters of the ocean. The distribution pattern of Thaumarchaeota in the deep sea in this study is consistent with previous studies (Baker et al., 2020), indicating that they play a significant role in the biogeochemical cycles of the deep sea. Additionally, we found that Thaumarchaeota is also abundant in some water bodies from the estuarine area, which may be related to the high nitrogen load in the estuarine region.

From the Northern Hangzhou Bay to the East China Sea, the spatial differentiation exhibited by samples presented an obvious correlation with the salinity, and it seems that these factors pose a dramatic impact on archaeal community composition. As shown in several studies of coastal zone, salinity mainly affects the microbial diversity variation (Fortunato et al., 2012; Wang et al., 2015; Liu et al., 2018). Salinity, as an environmental factor, mainly determines the community composition in soil, sediment, freshwater, and marine, and exerts a more important effect compared with pH, temperature, etc (Lozupone and Knight, 2007; Liu et al., 2015). Besides, regarding the coastal environment, salinity is likely to cause a density gradient that separates the residential microbial community in water masses (Fortunato et al., 2012). Here, we also found the importance of nutrient-related variables (primarily TN, TP, and NH₄ ⁺) and pH also affected the archaeal community. Nutrients remarkably affect the composition and growth of bacteria in the environment, including nitrogen, carbon, and iron in the bioavailable form, which may be fully utilized by microorganisms and closely related to the composition of the archaeal community (He et al., 2016). In conclusion, a deterministic process (environment filtering) assisted in shaping the archaeal community from Hangzhou bay to the adjacent sea.

In river estuary and its adjacent areas, water flows together with tidal activity mainly dominate the microbiota spread as well as the microbial transportation to the benthic zone (Augspurger et al., 2010). Therefore, because the dispersal rate (mass effect) is high, the archaeal composition may be randomly distributed within this highly connective ecological system. On that account, geographic distance may slightly affect the taxon diversification in the zone. Nonetheless, based on the null model, a stronger role of stochastic processes rather than deterministic processes in governing the archaeal community assembly was detected. In the middle-salinity samples, the null model revealed that stochastic processes were overwhelming. In contrast in low- and high-salinity sites, the influence of stochastic processes was weak. This was consistent with the result of CCA, where the spacial factor (PCNM) plays a more important role in structuring the microbial community in middle-salinity sites than in other regions. A potential explanation is that middle-salinity regions have lower regional connectivity, Due to the interactive effects between rivers and oceans (Zinger et al., 2014), in the present study, the water bodies in the medium-salinity area exhibit less restricted diffusion and mixed archaeal communities from both rivers and oceans, allowing the stochastic process to dominate the community assembly.

Here we found that the stochasticity of the archaeal community was mainly controlled by dispersal limitation and drift. The result is consistent with Wu and Huang (2019) reported that dispersal limitation and drift accounted for ~50% of the protist community turnover in the South China Sea and homogeneous dispersal explained 0%. In addition, these indicators have obvious dynamic changes along the salinity gradient. The influences of dispersal limitation had the lowest proportion in low salinity sites, this might be due to the strong runoff effect of the Qiantang River, allowing archaeal communities to disperse randomly. Therefore, environmental gradients (such as salinity gradient) and environmental conditions (such as hydrological conditions) might determine the balance of determinism and stochasticity, and thus construct the geographical pattern of archaeal community structure from Northern Hangzhou Bay to the East China Sea.

A microbial community network sees the occurrence of complicated interactions, which could be reflected in the co-occurrence networks to some extent (Faust and Raes, 2012; Liu et al., 2019). In this study, the archaeal communities exhibited a connectedness as well as a non-random co-occurrence. In addition, the network structure was significantly subjected to salinity, which further demonstrated the impact of the salinity (or other environmental factors of water) on the community structure, possibly because of the co-occurrence of dominant taxa with other taxa which had the same requirement for habitat (Delgado-Baquerizo et al., 2018). Functions possessed by nodes in different modules were different (Newman, 2006), and taxas strongly linked in ecology could cluster as one module. The network characteristics of archaeal communities within the same module exhibit noticeable variations across different salinity zones, thus, our results suggest that community compositions and ecological functions might be interconnected along the salinity gradient of the water. Also, it is necessary to pay attention to the case that the estimated co-occurrence pattern variation obtained by using a system method based on topology was incapable of reflecting the actual correlation among taxa (Ma et al., 2016).

We further reconstructed microbial biogeochemical processes in this salinity gradient region based on the FAPROTAX predictions. The ubiquity of aerobic ammonia oxidation and nitrification highlights the central role of archaea in the nitrogen cycle. This finding is consistent with previous studies on the importance of archaea in nitrogen transformation processes, further emphasizing their role in the global nitrogen cycle (Baker et al., 2020). Particularly in the medium salinity zone, the abundance of chemoheterotrophic microorganisms reaches the highest, which may be related to the specific environmental conditions of this area (such as salinity, temperature, and the availability of organic matter) as well as the ecological characteristics of the archaeal community (such as a lower deterministic process). This diversity and abundance distribution of trophic types may have profound effects on the function and stability of the ecosystem, as they affect energy flow and material cycling (Escalas et al., 2019; Lu et al., 2023). Moreover, from low to high salinity regions, the enhancement of aerobic ammonia oxidation and nitrification functions may be related to the increase in environmental stress, such as the rise in salinity that may promote the microbial adaptability of these processes. In contrast, the reduction of methylotrophy and methanogenesis may reflect the decreased energy efficiency of these metabolic pathways under high salinity conditions. This trend is significant for understanding the adaptive strategies of microbial communities in extreme environments. It is also worth noting that the detection of various enzymatic processes, such as chitinolysis, cellulolysis, and xylanolysis, indicates that archaeal communities play a key role in the decomposition of complex organic matter. These processes are crucial for maintaining the ecosystem’s carbon cycle and may respond to global climate change.

In summary, this study provides an in-depth understanding of the metabolic patterns of archaeal communities in the study area and emphasizes their key role in biogeochemical cycles. These findings lay the foundation for further exploration of the functions of archaea in ecosystem services and their responses to environmental changes. Additionally, in this study, spatial gradients were not found in archaeal community functional characteristics. This suggested that the geochemical functions of the archaeal communities are not influenced by the gradients of environmental factors in water. Furthermore, considering that nitrogen cycling pathways in sediments are mediated by a series of microorganisms occupying different taxonomic levels (Marcel et al., 2018), thus the taxonomic characteristics change might affect the potential of nitrogen cycling functions. Here, our results showed that the archaeal diversity characteristics (e.g. Chao 1 index) affect the spatial pattern of functional characteristics, except with the environmental factors. Suggested the unneglectable influences of microorganisms in the key element transformation in marine ecosystems. Furthermore, although further research is needed, we suggest that the balance of niche and neutral processes within microbial communities might influence their functional characteristics.

To sum up, based on the study results, the spatial distribution exhibited by archaeal community and function is subject to the direct control of environmental factors and the indirect control of the stochastic genetic drift and dispersal limitation. Future study is suggested to conduct sampling in time series and at least in the seasonal scale, for predicting dynamic archaeal assembly patterns on the time scale. It is necessary to integrate more detailed physiochemical information as well as hydrologic conditions with the functional information and taxonomy of microbe, for better characterizing the biogeographic community distribution as well as the functional structure. Although the predicted microbial function based on 16S rRNA do not reflect true abundance and are not always expressed in proportion, the annotated gene abundance seems to demonstrate the potential biogeochemical process instead of the microbiome-mediated actual biogeochemical process. Multi-omic profiling is needed to gain holistic insights into the genetic potential and main elements (eg: nitrogen, sulfur, and carbon) metabolic coupling of marine water microbial communities in future research. Anyway, our study provides an overall insight into archaeal diversity and functionality in this coastal ecosystem. This might provide essential basic data for the protection and restoration of estuarine and seawater ecosystems.