Phytoplankton and Microzooplankton Community Structure and Assembly Mechanisms in Northwestern Pacific Ocean Estuaries with Environmental Heterogeneity and Geographic Segregation

Understanding the underlying mechanisms shaping a metacommunity is useful to management for improving the ecosystem function. The research presented in the manuscript mainly tried to address the effects of habitat geography and environmental conditions on the phytoplankton and microzooplankton communities, and the underlying mechanisms of community assembly in temperate estuaries. ABSTRACT Phytoplankton and microzooplankton are crucial players in marine ecosystems and first responders to environmental changes, but their community structures and how they are shaped by environmental conditions have rarely been studied simultaneously. In this study, we conducted an eDNA metabarcoding sequencing combined with multiple statistical methods to simultaneously analyze the phytoplankton and microzooplankton in Liaohe (LH) and Yalujiang (YLJ) estuaries. The major objective was to examine how plankton community structure and assembly mechanism may differ between two estuaries with similar latitudinal position and climate but geographical segregation and differential level of urbanization (more in LH). Clear differences in diversity and composition of phytoplankton and microzooplankton communities between LH and YLJ estuaries were observed. Richness of phytoplankton was significantly higher in LH than YLJ, while richness of microzooplankton was higher in YLJ. The magnitude of intrahabitat variations in phytoplankton communities was significantly stronger than that of microzooplankton. Some phytoplankton and microzooplankton taxa also showed interhabitat differences in their relative abundances. Phytoplankton showed a stronger geographic distance-decay of similarity than microzooplankton, while significant environmental distance-decay of similarity in microzooplankton was found in the less urbanized YLJ estuary. Community assembly of phytoplankton was, based on the neutral community models, driven primarily by stochastic processes, while deterministic processes contributed more for microzooplankton. Furthermore, we detected wider habitat niche breadths and stronger dispersal abilities in phytoplankton than in microzooplankton. These results suggest that passive dispersal shapes the phytoplankton community whereas environmental selection shapes the microzooplankton community. IMPORTANCE Understanding the underlying mechanisms shaping a metacommunity is useful to management for improving the ecosystem function. The research presented in the manuscript mainly tried to address the effects of habitat geography and environmental conditions on the phytoplankton and microzooplankton communities, and the underlying mechanisms of community assembly in temperate estuaries. In order to achieve this purpose, we developed a metabarcoding sequencing method based on 18S rRNA gene. The phytoplankton and microzooplankton communities from two estuaries with similar latitude and climatic conditions but obvious geographical segregation and significant environmental heterogeneity were investigated. The results of our study could lay a solid foundation for ascertaining phytoplankton and microzooplankton communities in estuaries with obvious environmental heterogeneity and geographic segregation and mechanisms underlying community assembly.

phytoplankton and microzooplankton communities simultaneously by metabarcoding sequencing of 18S rRNA gene in the LH and YLJ estuaries areas. Variations in assembly mechanisms of phytoplankton and microzooplankton communities between different studied areas were also evaluated. We hypothesized that the habitat heterogeneity and geographic segregation of LH and YLJ estuaries could induce distinct community assembly mechanisms and result in different plankton communities.

RESULTS
Annotation of phytoplankton and microzooplankton. According to the sequencing and data processing pipeline (Fig. 1b), a total of 3,881,365 sequencing reads were obtained for 33 water samples in this study (Table S1), which clustered into 635 high-quality ASVs. Among them, 76 and 149 ASVs were annotated as phytoplankton and microzooplankton, respectively, and these accounted for 5,897 to 33,809 reads for phytoplankton and 1,908 to 44,845 reads for microzooplankton (Table S1). All phytoplankton and microzooplankton ASVs were successfully annotated at the phylum level (Fig. S1). A total of 6 and 20 phytoplankton phyla and microzooplankton lineages were identified, ranging from 2 to 5 and 1 to 6 for each sample, respectively (Table S2 and S3). Pyrrophyta, Bacillariophyta, and Charophyta accounted for most of the phytoplankton ASVs (Fig. S2a), whereas Ciliophora was the dominant microzooplankton lineage (Fig. S2b). About 75% and 50% ASVs were assigned phytoplankton and microzooplankton genera, respectively, but only 17.1% phytoplankton and 25.5% microzooplankton ASVs can be assigned to a species (Fig. S1). In addition, rarefaction curves based on the detected species numbers were all close to the plateau (Fig. S3), demonstrating that the amounts of effective sequencing reads were sufficient to represent the whole phytoplankton and microzooplankton communities. Moreover, the species accumulation curves were close to a straight line (Fig. S4), indicating that the number of samples collected in this study can fully represent the phytoplankton and microzooplankton communities in LH and YLJ estuaries.
Differences in diversity and compositions of phytoplankton and microzooplankton communities between LH and YLJ estuaries. Phytoplankton communities in the LH estuary had significantly higher Chao1, Pd_faith, and Shannon indices than those in the YLJ estuary (t test, P , 0.05, Fig. 2a). In contrast, the Chao1 and Shannon indices of microzooplankton communities in YLJ were significantly higher than those in LH (t test, P , 0.05, Fig. 2a). No significant difference was found in Pielou_J indices between LH and YLJ for either phytoplankton or microzooplankton communities (t test, P . 0.05, Fig. 2a). These results indicated that there was no difference in community evenness between LH and YLJ, but the richness of phytoplankton and microzooplankton was higher in LH and YLJ, respectively.
Bray-Curtis distance was used to evaluate the variations of phytoplankton and microzooplankton communities between the LH and YLJ estuaries. Communities from different habitats were separately clustered and distributed on different sides of PC1, which explained 75% and 33% of the total variations in phytoplankton and microzooplankton, respectively (Fig. 2b). Significant between-region differences in composition of both phytoplankton and microzooplankton communities were also confirmed (Adonis test, P , 0.05). The Bray-Curtis distances of phytoplankton and microzooplankton communities in LH were significantly larger than those in YLJ (t test, P . 0.05, Fig. 2c). Moreover, intrahabitat distances of phytoplankton and microzooplankton communities in LH were similar, while in YLJ, the distances of microzooplankton communities were more pronounced relative to distances of phytoplankton communities (Fig. 2b and c). These findings indicated the differences in phytoplankton and microzooplankton compositions between the LH and YLJ estuaries and stronger variations in phytoplankton.
Geographic differences in phytoplankton and microzooplankton. For phytoplankton, Pyrrophyta and Chlorophyta were dominant phyla in the LH and YLJ estuaries, respectively, followed by Bacillariophyta (Fig. 3a). The relative abundances of Pyrrophyta, Bacillariophyta, and Charophyta were significantly higher in LH than YLJ, whereas Chlorophyta was significantly more abundant in YLJ (Wilcox Rank-Sum test, P , 0.05, Fig. 3b). At the genus level, Paragymnodinium and Ostreococcus were the most dominant genera in the LH and YLJ estuaries, respectively (Fig. S5a). A random forest model was further constructed to recognize key taxa for predicting the origins of phytoplankton communities based on the genus-level Microbiology Spectrum data sets. Five phytoplankton genera, including Ostreococcus, Bathycoccus, Micromonas, Thalassiosira, and Syndiniales Group II, were identified to distinguish the phytoplankton communities between the LH and YLJ estuaries (Fig. 3c). The random forest model based on phytoplankton genera showed 93.94% accuracy for discrimination LH and YLJ samples (Fig. S6a). For microzooplankton, Ciliophora was the dominant lineage in both the LH and YLJ estuaries, followed by Copepoda, Tetrapoda, and Hexacorallia in LH, and Trachylinae and Pteriomorphia in YLJ (Fig. 3d). Trachylinae and Sagittoidea were significantly more abundant in LH estuaries (Wilcox Rank-Sum test, P , 0.05, Fig. 3e). In contrast, significantly higher abundances of Ciliophora, Copepoda, and Descapodiformes were found in YLJ (Wilcox Rank-Sum test, P , 0.05, Fig. 3e). At the genus level, Mesodinium and Trachymedusae were the most dominant microzooplankton genera in the LH and YLJ estuaries, respectively (Fig.  S5b). Based on the Random Forest model, six microzooplankton genera were chosen, including Trachymedusae, Spirotontonia, Myrionecta, Strombidinopsis, Variatrombidium, and Tintinnidium to distinguish the microzooplankton communities between the LH and YLJ estuaries (Fig. 3f). The model recorded an accuracy of 100% for the microzooplankton community differentiation (Fig. S6b).
Distance-decay of similarity for phytoplankton and microzooplankton communities. Distance-decay of similarity of phytoplankton and microzooplankton communities based on environmental variables and geographic distances were evaluated in LH and YLJ, respectively (Fig. 4a). Significant distance-decay of similarities of geographic distances were found in phytoplankton communities from both LH and YLJ estuaries (Table S4, linear regression, P , 0.05). Among them, the slopes of the LH group were higher than that of the YLJ group (Table S4). In contrast, the distancedecay of similarities between environmental variables and phytoplankton communities were not significant (Table S4, linear regression, P . 0.05). For microzooplankton communities, both significant environmental and geographic distance-decay similarities were found in YLJ estuaries but not in LH samples (Table S4). Moreover, the power of distance-decay similarity for geographic distance was stronger than that for environmental variables in microzooplankton communities from YLJ (Table S4). VPA was further applied to quantify the relative contributions of environmental selection and passive dispersal on the variations of phytoplankton and microzooplankton communities (Fig. 4b). Regarding variations in phytoplankton communities, the pure effects of spatial variables were stronger in both LH and YLJ estuaries than those of environmental variables (16.9% versus 1.7% in LH and 22.0% versus 0.8% in YLJ). Regarding variations in microzooplankton communities, the pure effect of environmental variables increased; however, it was also lower than that of spatial variables (7.9% versus 11.3% in LH and 7.7% versus 14.5% in YLJ). These results indicated that the phytoplankton and microzooplankton communities in the LH and YLJ estuaries were strongly more governed by passive dispersal than environmental selection.
Community assembly of phytoplankton and microzooplankton. To explore the mechanisms responsible for the observed community patterns, neutral community model was performed (Fig. 4c). In the LH and YLJ estuaries, 65.7% and 77.1% phytoplankton ASVs, respectively, can be explained by the neutral community models, indicating that stochastic process dominated the assembly of phytoplankton communities. The explained ratios of neutral community models were much lower for microzooplankton communities (43.2% and 53.8% in LH and YLJ, respectively). These results indicated that deterministic process, such as environmental selection, contributed more to microzooplankton communities. Moreover, the habitat niche breadth and dispersal ability for phytoplankton communities were significantly higher than those for microzooplankton communities in both LH and YLJ estuaries (t test, P , 0.05, Fig. 4d). Phytoplankton and microzooplankton in the LH and YLJ estuaries can be separated into three fractions comprising those ASVs found above, below, and neutral partitions (Fig. 5). The proportions of phytoplankton and microzooplankton ASVs showed that neutral fraction accounted for much higher richness than the above and below fractions in both LH and YLJ estuaries. However, the abundance proportions of phytoplankton ASVs were dominant by the above fraction (67.72% in LH and 90.64% in YLJ). In contrast, although the abundance proportions of microzooplankton ASVs were dominant by the neutral fraction, the abundance proportions for above and below fractions also increased compared to those for richness (Fig. 5).
Taxa above the prediction are supposed to have higher migration ability and can Plankton Community Assembly in Estuaries Microbiology Spectrum disperse to more locations. Below the prediction are taxa with lower dispersal ability in the estuaries. For phytoplankton in LH, Charophyta, Chlorophyta, Pyrrophyta occupied the richness of the above fraction (Fig. 6a). In comparison, the abundance proportions of Chlorophyta and Pyrrophyta for the above fraction were notably higher (Fig. 6b), indicating that these two phyla contributed most to the migration of phytoplankton in LH. The same phenomenon was also observed for Chlorophyta in YLJ (Fig. 6). In addition, the abundance proportion of Pyrrophyta in YLJ estuary was also higher than richness proportion of it (Fig. 6). According to the results of microzooplankton, the migration of microzooplankton in the LH and YLJ estuaries were contributed by Ciliophora and Trachylinae, respectively (Fig. 6). Moreover, Copepoda and Hexacorallia in LH as well as Pteriomorphia in YLJ were found to be less dispersible (Fig. 6).

DISCUSSION
Phytoplankton and microzooplankton are foundational components in aquatic Plankton Community Assembly in Estuaries Microbiology Spectrum ecosystems where they play a key role in energy flows and biogeochemical cycles. Ecological processes that shape the phytoplankton and microzooplankton communities are among most longstanding issues in microbial ecology. However, the driving factors and assembly mechanisms of sympatric phytoplankton and microzooplankton communities remain poorly understood, especially in vulnerable aquatic ecosystems affected by strong natural and anthropogenic disturbances. Recent development of high-throughput sequencing based on eDNA metabarcoding provides a feasible option for monitoring the microbial communities with fast, standardized, and comprehensive features (29). The sequencing data were used to investigate the diversity, composition, and assembly mechanisms of both phytoplankton and microzooplankton communities in estuarine areas of main rivers in northeast China with remarkable environmental heterogeneity and geographic segregation. Despite the successful application of eDNA metabarcoding technology in the present study, it still has some shortfalls. The lack of a thorough barcode data set to base the eDNA analysis is the biggest obstacle for application the eDNA metabarcoding technology (30), which could lead to false negative and a few wrong taxonomic results (31). In addition, further comparison between the results from classical analysis with the use of microscopes could benefit to evaluate the accuracy of eDNA metabarcoding for microplankton (28).
Effects of environmental heterogeneity on the phytoplankton and microzooplankton communities in estuarine ecosystems. According to the results of alpha diversity indices, the richness of phytoplankton was significantly higher in the LH estuary, while that of microzooplankton was higher in YLJ (Fig. 2a). Phytoplankton generally thrives in nutrient-rich environments (29), to which the LH estuary belongs (Fig. S7). Positive

Plankton Community Assembly in Estuaries
Microbiology Spectrum correlations between phytoplankton richness and some environmental variables, such as COD and phosphate, observed in this study further confirmed the effects of nutrients on phytoplankton communities (Fig. S8a). Some phytoplankton species can be used as a food source for microzooplankton (32), thus phytoplankton and microzooplankton abundances in the same ecosystem usually show the opposite trends (33). Promotion of microzooplankton can lead to a sustained decrease in phytoplankton, despite high nutrient levels (34)(35)(36). Besides, environments with high microzooplankton abundances always had a higher level of ammonia nitrogen due to the limited phytoplankton (37,38), which can consume almost all the ammonia nitrogen when phytoplankton form blooms (39). This phenomenon was consistent with the higher ammonia concentration and microzooplankton richness in the YLJ estuary compared to the LH samples (Fig. S8b). Moreover, multiple heavy metals were found to be positively correlated with the phytoplankton richness in the LH and YLJ estuaries (Fig. S8). In a previous study at the Bohai Sea, effects of heavy metal pollution on phytoplankton communities and hence resulting in variations of phytoplankton communities had also been proposed (17). Furthermore, the PCoA illustrated that phytoplankton and microzooplankton communities from the same habitat clustered together (Fig. 2b), indicating distinct community composition from different habitats. Environmental heterogeneity and geographic segregation may explain this phenomenon. First, this result might have resulted from the differences in environmental conditions between different habitats (40). Based on the same samples from this study, we have shown significant environmental heterogeneity between LH and YLJ estuaries (23). Significantly higher concentrations of PHs and multiple heavy metals were detected in the LH estuary compared to the YLJ (t test, P , 0.05, Fig. S7). In contrast, the YLJ samples possessed significantly higher ammonium concentrations (t test, P , 0.05, Fig. S7). In addition, water temperature and salinity were two and six degrees lower in the YLJ estuary compared to LH (t test, P , 0.05, Fig. S7). Second, geographical segregation could result in the presence of regional species pool and hence different community composition (41). Only 14.5% phytoplankton ASVs and 16.1% microzooplankton ASVs were found to be shared between the LH and YLJ estuaries (Fig. S9), revealing huge differences in regional species pools between these two habitats.
Passive dispersal shaping phytoplankton community but environmental selection shaping microzooplankton community. Environmental and spatial factors are two major factors to influence the community assembly in natural environments (42). In the present study, several environmental variables were identified as having a significant effect on the phytoplankton and microzooplankton compositions by canonical correlation analysis (Table S5). Meanwhile, significant geographic distance-decay of similarities for phytoplankton and microzooplankton communities were also observed in both LH and YLJ estuaries (Fig. 4a). Although significant effects of environmental and spatial factors on the phytoplankton and microzooplankton communities were uncovered, results of VPA in this study indicated that most of the community variation remained unexplained by them in both estuaries (Fig. 4b). VPA is widely used in ecological research to determine the contribution of different factors on community variations; however, it will underestimate the role of environmental factors and should be applied together with other approaches (43). In this study, the neutral community model was used to adjust VPA and infer the influences of ecological processes.
Our results clearly support the prominent role of stochastic processes in shaping the phytoplankton community assembly of both LH and YLJ estuaries (Fig. 4c). According to the ecological theory, stochastic processes shaping the community assembly can be divided into dispersal and ecological drift (44). Based on the neutral theory, community similarity was predicted to decrease along spatial distances due to dispersal limitation when stochastic processes govern the community assembly (45). The significant and strong distance-decay pattern of phytoplankton communities in the LH and YLJ estuaries is indicative that dispersal was important (Fig. 4a). Comparatively, the relative importance of stochastic processes to shape the microzooplankton communities was much lower than that for phytoplankton based on the neutral community models (Fig. 4c). These results suggest deterministic selection by environmental factors shaped the microzooplankton communities due to the preferences and fitness of taxa (46), which were consistent with the higher contribution of environmental variables to microzooplankton communities shown by distance-decay of similarity and VPA ( Fig. 4a and b). Regarding the community immigration rate, the m values of neutral community models for phytoplankton in the LH and YLJ estuaries were higher than the microzooplankton (0.509 versus 0.054 in LH and 0.753 versus 0.464 in YLJ, Fig. 4c). These results were consistent with the average shared proportion of taxa in the LH and YLJ estuaries (Fig. 4e) and indicated the higher dispersal ability of phytoplankton than microzooplankton. Moreover, wider niche breadth means more metabolic adaptation and less environmental selection for species in a community (47). Significantly wider niche breadth of phytoplankton than microzooplankton in LH and YLJ estuaries further confirmed the importance of environmental selection on microzooplankton community assembly (Fig. 4d). Token together, our results revealed that passive dispersal shaping phytoplankton community but environmental selection shaping microzooplankton community in estuarine ecosystems.
Difference between neutral and nonneutral partitions. The neutral community model separated taxa into neutral and nonneutral partitions, which were possibly attributable to the different dispersal rates among taxa (43). Chlorophyta had the highest dispersal ability in phytoplankton communities from both the LH and YLJ estuaries (Fig. 6), which was also more abundant in YLJ estuary (Fig. 3b). Chlorophyta has been reported to have higher dispersal ability because its spores and resistant vegetative fragments facilitate dispersion and desiccation tolerance (48). In addition, higher ammonia concentrations in YLJ could be the reason for Chlorophyta to bloom in this area (49). In addition, Pyrrophyta were more abundant in LH and played a vital role in phytoplankton dispersal (Fig. 3b and 6). Pyrrophyta typically have two flagella that allow cells to move actively (50). Relative water temperature and salinity in the LH estuary could be the reason of higher abundances of Pyrrophyta (51). Besides, petroleum pollution, one of the main contaminations in LH estuary, has been shown to cause the initiation of Pyrrophyta blooms (52). For microzooplankton, Ciliophora were found to be the main mobile microzooplankton in the LH estuary (Fig. 6), which were also more abundant in LH estuary compared to YLJ samples (Fig. 3e). In accordance with our results, previous studies have revealed that Ciliophora were positively correlated with concentrations of nutrients (53,54). These suggest that Ciliophora in estuaries could actively move to locations with high nutrients. Trachylinae was the microzooplankton lineage with the highest degree of dispersal in the YLJ estuary (Fig. 6). At the landward side of the YLJ estuary, many aquaculture ponds are in operation, and jellyfish, a subclass of Trachylinae, is one of the main cultivated organisms (55). Larvae entering the estuarine area through water exchange may be the reason for the higher abundance of Trachylinae in the YLJ estuary. At the same time, jellyfish is a highly mobile taxa and can move to areas farther offshore (56).
Conclusions. Our study demonstrated the feasibility of eDNA meabarcoding sequencing to simultaneously capture the diversity and composition of phytoplankton and microzooplankton in the aquatic environments. Based on the diversity comparisons for the LH and YLJ estuaries, richness of phytoplankton was significantly higher in LH than in YLJ, while microzooplankton richness was higher in YLJ. Remarkable differences of phytoplankton and microzooplankton compositions between the LH and YLJ estuaries were observed, which could be due to environmental heterogeneity and geographic segregation between these two regions. Neutral community models showed that stochastic processes dominated the phytoplankton community assembly in both estuaries. In contrast, the relative importance of stochastic processes in microzooplankton community assembly was relatively low. Furthermore, results of distance-decay similarity and VPA indicated that passive dispersal shaped phytoplankton community but environmental selection shaped microzooplankton community in the two estuaries. Phytoplankton and microzooplankton taxa that did not conform to the neutral distribution were found in LH and YLJ estuaries, which usually have active mobility and might be selected by local environmental conditions. Based on the investigation of model estuarine systems, our results reveal the differences in planktonic communities between estuaries at similar latitudes with similar climate conditions but with different geomorphologies and anthropogenic impacts. Furthermore, our in-depth ecological model analyses shed lights on the underlying assembly mechanisms of planktonic communities, highlighting the importance of environmental heterogeneity and geographic segregation on planktonic community assemblages.

MATERIALS AND METHODS
Sample collection and measurements of environmental variables. A total of 15 and 18 surface water samples were collected from LH and YLJ estuaries, respectively, in July 2019 (Fig. 1a). For each sample, approximately 2 L of surface water (;0.5 m depth under the water) was collected with a water sampler (Wuhan Shuitiandi Instruments, Wuhan, China) and temporarily held on ice for less than 6 h before filtration. For DNA samples, 1 L of each sample was filtered through 0.2-mm pore size polycarbonate membranes (142 mm diameter, Millipore, USA) with a peristaltic pump. The filters were immediately preserved in liquid nitrogen, transported to laboratory, and stored at 280°C until DNA extraction. The remaining 1 L of water was placed in a clean and sterile glass container for environmental variables analyses in the laboratory. A total of 15 environmental variables, including water temperature (Temp), salinity, dissolved oxygen (DO), chemical oxygen demand (COD), petroleum hydrocarbons (PHs), nitrate, nitrite, ammonium, phosphate, and six heavy metals (Cu, Pb, Zn, Cd, Cr, and As) were measured at each sampling station as previously described (23).
DNA extraction and sequencing. Total DNA for each water sample was extracted from the filters using a PowerWater DNA isolation kit (Qiagen, CA, USA) according to the manufacturer's instructions. Agarose gel electrophoresis (1% concentration) was used to detect the success of DNA extraction. DNA qualities were measured using a NanoDrop ND-1000 Spectrophotometer (NanoDrop, USA). The 18S rRNA V9 region from each extracted DNA were amplified using primers 1380F-1510R (1380F: TCCCTGCCHTTTG TACACAC, 1510R: CCTTCYGCAGGTTCACCTAC) (57). PCR amplification and sequencing were done as previously described (23). Raw data have been deposited to the US National Center for Biotechnology Information (NCBI) under the BioProject number PRJNA782015.
Data processing. Reads with an average Phred score (Q score) lower than 20, ambiguous bases, a homopolymer greater than 6, and mismatches in primers were deleted (58). Remaining high-quality reads were ascribed to the samples based on their unique barcodes at the end of reverse primers. These reads were then assigned to amplicon sequence variants (ASVs) with the Quantitative Insights Into Microbial Ecology 2 (QIIME2) through the Divisive Amplicon Denoising Algorithm 2 (DADA2) method (59). Singletons (number of a specific ASV = 1) were abandoned to improve the efficiency of data analysis. Taxonomic annotation of ASVs was performed using the QIIME2 pipeline based on the SILVA 138 database (60). ASVs belonging to phytoplankton and microzooplankton were binned separately using an in-house R script, and then the taxonomic names at phylum level were rematched. Finally, the phytoplankton and microzooplankton ASV abundance tables were normalized using a standard number of reads according to the sample with the least number of reads (5,897 and 3,966 tags for phytoplankton and microzooplankton, respectively).
Statistical analysis. All statistics analyses were performed using R v4.0.2 platform. Rarefaction and species cumulative curves of sequenced data sets were extrapolated by "iNEXT" package. Four alpha diversity indices of phytoplankton and microzooplankton communities, including Chao1, Shannon, Pielou_J, and Pd_faith, were calculated by "vegan" package. Variations in the composition of phytoplankton and microzooplankton were evaluated by principal coordinate analysis (PCoA) and the Adonis test based on the Bray-Curtis distance using the "vegan" package. Differences in alpha diversity indices, and community distances between samples from LH and YLJ estuaries were analyzed using a Student's t test (t.test function). Differences in relative abundances of phytoplankton and microzooplankton taxa between LH and YLJ estuaries were analyzed using the Wilcox Rank-Sum test (wilcox.test function). Results of differentially abundant taxa and relative abundances of dominant taxa of phytoplankton and microzooplankton communities were visualized by "ggplot2" package. Furthermore, Random Forest analyses were performed using the "RandomForest" package to recognize key phytoplankton and microzooplankton genera for distinguishing the origins of samples. The relative abundances of these key taxa and the accuracy of the Random Forest model were displayed using the "ggplot2" and "pheatmap" packages, respectively.
Euclidean distance of environmental variables ("vegan" package) and geographic distance ("geosphere" package) between any two sampling stations were calculated, and then linear regression (lm function) was applied to evaluate the distance-decay of similarities of phytoplankton and microzooplankton communities. In addition, variation partitioning analysis (VPA) was completed by "vegan" package to determine the relative importance of environmental and spatial factors in shaping phytoplankton and microzooplankton communities in the LH and YLJ estuaries. Moreover, a neutral community model was used to determine the potential importance of stochastic process on assembly of phytoplankton and microzooplankton communities. In this model m is an estimate of dispersal between communities and R2 represents the ratio of contribution by stochastic processes (61). The ASVs from each data set were subsequently separated into three partitions depending on whether they occurred more frequently than (above partition), less frequently than (below partition), or within (neutral partition) the 95% confidence interval of the neutral community model predictions. The richness and composition of these separated partitions were further compared, and the results were displayed by "ggplot2" package.
Furthermore, niche breadth and dispersal ability are critical traits that influence the relative importance of environmental selection and passive dispersal (47,62). Niche breadth of phytoplankton and microzooplankton communities in the LH and YLJ estuaries was measured by the "spaa" package based on Levins' niche breadth index. The average B-values from all taxa in a single community (Bcom) as an Plankton Community Assembly in Estuaries Microbiology Spectrum indicator of habitat niche breadth at the community level. A wider niche breadth of an organism group means more metabolically flexible at the community level (47). To estimate dispersal ability of each taxon (passively driven by water movements), we calculated the pairwise shared proportion of sequence numbers, and used the average shared proportion as a proxy for dispersal. A higher shared proportion of sequences indicates a more successful passive transport by currents (63). The Student's t test was performed to evaluate the differences in niche breadth and dispersal ability between the phytoplankton and microzooplankton communities in the LH and YLJ estuaries, respectively. Data availability. All raw sequences of sediment bacterial communities studied in this study have been submitted to the NCBI Sequence Read Archive (SRA) database under the BioProject number PRJNA782015.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 2 MB.

ACKNOWLEDGMENTS
This study was supported by the National Key Research and Development Program of China (2020YFA0607600).
We have no conflicts of interest to declare.