Structure and Functional Diversity of Surface Bacterioplankton Communities in an Overwintering Habitat for Large Yellow Croaker, Pseudosciaena crocea, of the Southern East China Sea

Given the biogeochemical functions of marine microbes, the structure and functional diversity of bacterioplankton communities could be regarded as bioindicators that reflect environmental quality. The habitat of the large yellow croaker (Pseudosciaena crocea) is suffering from widespread degradation, but the study of the croaker habitat is still insufficient. Here, we used 16S rRNA gene sequencing and functional prediction to investigate the surface bacterioplankton communities in an overwintering ground of P. crocea during the overwintering period (November 2018 to April 2019). Our results showed that the taxa and functions of the environmental bacterioplankton community exhibited obvious seasonality. Specifically, the bacterioplankton communities in autumn were characterized by animal parasites or symbionts, such as Acinetobacter baumannii and Moraxella, while organic substance–degrading bacteria, such as Rhodococcus, Stenotrophomonas, and Alcanivorax, dominated in spring. Moreover, phylogenetic-based mean nearest taxon distance (MNTD) analyses indicated that dispersal limitation was the most important process governing the spatiotemporal assembly of the bacterioplankton communities. The distance decay of similarity also suggested the impact of dispersal limitation on the generation of biogeographical patterns. Furthermore, variation partitioning analysis (VPA) and a partial Mantel test revealed that environmental filters, such as nutrients, temperature, and salinity, had relatively weak deterministic effects on the bacterioplankton community assembly. Overall, these findings provide a basis for understanding the theory of surface microbial community assembly in the overwintering habitat of the large yellow croaker.


INTRODUCTION
Abundant fisheries resources are indispensable cornerstones for the sustainable development of the coastal economy (Tursi et al., 2015;Keen et al., 2018). However, as China's most productive fishing area, the East China Sea has been suffering from resource degradation in recent decades due to increasing human impacts, such as overexploitation, pollutant discharge, and climate change (Wang et al., 2012). The Chinese government has implemented a series of strict fishing moratoriums and stock management systems since 1995, and the current status of most natural resources has generally improved. Despite all these programs, the wild stocks of large yellow croaker (Pseudosciaena crocea), which is known as China's national fish, have not recovered (Liu and de Mitcheson, 2008). Evidence strongly suggests that coastal pollution, habitat degradation, and possible ecosystem changes may have compromised the potential for stock recovery (Liu and de Mitcheson, 2008). However, the habitat (including the spawning, nursery, and overwintering grounds and migration channels) of the large yellow croaker has not been sufficiently researched. According to inferred migratory routes, croaker aggregate in overwintering grounds from October to April each year (Liu and de Mitcheson, 2008;Xu and Chen, 2011;Chen and Xu, 2012). Because the low temperature may suppress croaker immune functions by influencing multiple physiological processes and thereby make the croaker susceptible to the surrounding environment (Yang et al., 2015;Liu et al., 2019), the habitat conditions of their overwintering grounds are vital for recovering their stock and supporting this marine resource.
Water microorganisms, especially planktonic bacteria, are one of the largest reserves of biodiversity and promote various ecological processes in marine ecosystems (Sunagawa et al., 2015). On the one hand, bacteria have long been used as sensitive and aggregative indicators of changes in environmental conditions due to their wide distribution, short generation time, and diverse functions (Dai et al., 2017;Yang et al., 2018). On the other hand, the bacterioplankton community includes a large number of probiotic organisms and pathogens that can directly affect fish health (Xiong et al., 2015;Hu et al., 2017). Effects of bacterioplankton on fish could also be indirect through the food chain (Kalcheva et al., 2008). For example, a microbial loop mediated by bacteria could affect the upper elements of the food chain of the croaker. Therefore, an accurate representation of the overwintering ground's current conditions can be systematically derived by observing the marine bacterioplankton community. Given the important roles of the bacterioplankton community in the assessment of the habitat environment, it is necessary to reveal the fundamental processes that underlie biogeographic and ecological patterns of bacterioplankton diversity. However, this field of study remains a major challenge in bacterial ecology (Jiao et al., 2019). Historically, bacterial community assembly has been studied from a deterministic perspective, for example, environmental filtering and biotic interactions (Dini-Andreote et al., 2015). However, an increasing number of recent studies support the idea that stochasticity plays a predominant role in some microbial systems (Zhou et al., 2013;Chen W. D. et al., 2019), especially in the biogeographical patterns of marine bacterioplankton systems (Giovannoni and Nemergut, 2014). In fact, community assembly is affected by deterministic and stochastic processes simultaneously, but the relative effects of these two processes vary among different ecosystems. Various methods have been developed to separate the relative effects of deterministic and stochastic processes. Among them, a phylogenetic framework that uses information about phylogenetic relatedness and phylogenetic dissimilarity derived from phylogenetic trees has been widely used (Stegen et al., 2013. Using this framework will help us to reveal the potential factors underlying the spatiotemporal patterns of bacterial communities in the overwintering grounds for croaker. In addition, the focus of bacterial ecology is turning from its initial focus on pattern description toward functional profiling. The Functional Annotation of Prokaryotic Taxa (FAPROTAX), which focuses on the biogeochemical cycle functions of microorganisms, is commonly used (Louca et al., 2016). This method depends on prokaryotic functional data from culturable bacteria and species classification (OTU table) to obtain community function information and provides opportunities to understand the differences in the structure and function of bacterioplankton communities .
In the present study, we investigated the surface bacterioplankton communities in the overwintering grounds of P. crocea during the overwintering period using high-throughput sequencing of 16S rRNA genes, aiming at (1) revealing the spatiotemporal patterns of taxonomic and functional profiles of bacterioplankton communities in the overwintering grounds and (2) understanding the processes that underlie the observed patterns in species abundances across space and time. These results are essential for understanding microbial community assembly theory in this marine ecosystem and adding to the data on the environmental microbiome in the surface water of croaker habitat.

Study Site, Sample Collection, and Environmental Variables Analyses
Surface water samples (∼2 m) were collected in sterile containers from the overwintering ground for P. crocea (Figure 1), which is located in the East China Sea during 17 November-5 December 2018 (November) and 7-30 April 2019 (April). In total, we obtained 50 water samples (25 sites × 2 time points). The samples were stored in an icebox (4 • C) and immediately transported to the refrigerator at −20 • C. Environmental variables of water samples, including water temperature (WT), salinity (SAL), orthophosphate (PO 4 ), nitrite (NO 2 ), nitrate (NO 3 ), ammonium (NH 4 ), silicate (SiO 3 ), total nitrogen (TN), and total phosphorus (TP), were analyzed following the standard methods (AQSIQ, 2007) (Supplementary Table S1). Approximately 500 mL of a water sample from each site was filtered-sterilized onto one 0.2-µm pore-sized polycarbonate membrane (47 mm diameter, Millipore, Boston, MA, United States) on the sampling day. The membranes were immediately frozen at −80 • C until they were needed. FIGURE 1 | Study area and location of sampling stations in the overwintering ground. To keep the data comparable with future studies, we keep the discontinuous numbers for site names according to the routine-monitoring project. The sampling sites remain the same for both seasons, except for November at B2 and F2 and April at B1 and F1.

DNA Extraction, Bacterial 16S Amplification, and Illumina HiSeq Sequencing
Bacterioplankton DNA was extracted using a DNA extraction kit (MinkaGene Water DNA Kit). A NanoDrop One spectrophotometer was used to measure the DNA concentration and purity (Thermo Fisher Scientific, MA, United States). The V4 regions of the bacterial 16S rRNA genes were amplified (50 µL reaction volume) by thermocycling: 5 min at 94 • C for initialization; 30 cycles of 30 s denaturation at 94 • C, 30 s annealing at 52 • C, and 30 s extension at 72 • C, followed by 10 min final elongation at 72 • C using the primer sets 515F (5'-GTGCCAGCMGCCGCGGTAA-3') and 806R (5'-GGACTACNNGGGTATCTAAT-3') with barcodes (Parada et al., 2015). One pooled sample was combined through an equimolar amount of each PCR product from each sample and sequenced using an Illumina HiSeq 2500 platform to generate 250 bp paired-end reads (Guangdong Magigene Biotechnology Co., Ltd., Guangzhou, China).

Sequencing Data Processing
Sequenced paired-end reads were processed using USEARCH (Version 11) (Edgar, 2010). The paired reads were merged with the fastq_mergepairs function and filtered with a "maxee" value of 1.0 (Edgar, 2010). The unique sequence reads were obtained with the fastx_uniques function and then denoised using the UNOISE3 (unoise_alpha = 2 and minsize = 4 as per default settings) algorithm to correct errors and remove chimeras (Edgar, 2016). This generated a list of zero-radius operational taxonomic units (ZOTUs). A ZOTU abundance table was produced by mapping the total reads to the list of ZOTUs with the -otutab function with an identity threshold of 0.97. The representative sequences for each ZOTU were assigned to taxonomic groups using the RDP classifier within the SILVA database (release 138) clustered at 99% similarity. Raw sequences have been deposited in the NCBI Sequence Read Archive (SRA) database under the BioProject number PRJNA604077 and the accession number SRP247264. Functional profiles of the bacterial community were annotated using FAPROTAX 1.2.2 using python collapse_table.py with the rarefied ZOTU table (Louca et al., 2016).

Bacterioplankton Diversity and Composition Analyses
All statistical analyses were performed in the R software environment (version 3.5.2) unless otherwise noted. To improve the normality and homoscedasticity of the data, all of the biotic data were transformed by the Hellinger transformation, and the abiotic data were normalized by a chord transformation using the decostand function (Legendre and Gallagher, 2001). The Bray-Curtis dissimilarity matrix for bacterioplankton communities was generated to identify the relationships among samples with the vegdist function. Principal coordinate analysis (PCoA) and cluster analysis were used to visualize the overall structure of the bacterioplankton communities with the hclust (method = "ward.D") function. Permutational multivariate analysis of variance (PERMANOVA) was performed to quantitatively evaluate the contributions of seasons, latitude, and offshore distances to the variation in community structure with the adonis function. Analysis of similarity (ANOSIM) was performed to test whether there was a statistically significant difference in the community composition between the two seasons with the anosim function from the vegan package. To acquire the biomarker functions for each season, we used the randomForest (importance = TRUE, proximity = TRUE) function to generate a classification model for November and April. The biomarker functions were identified by the most important season-discriminating features using tenfold cross-validation implemented with the rfcv function .

Ecological Process Estimation and Environmental Driver Analyses
The quantification of the effects of ecological processes on bacterioplankton community assembly was performed according to Stegen et al. (2013) based on the turnover in the phylogenetic community composition and species composition. First, the observed nearest-phylogenetic-neighbor distance (βMNTD obs ) as well as the null expectation of the mean nearest-phylogeneticneighbor distance (βMNTD exp ) were used to calculate β-nearest taxon index (βNTI) metrics, which quantify the phylogenetic turnover between communities . The relative contributions of heterogeneous and homogeneous selection were estimated as the fractions of pairwise βNTI values that were >+2 and <-2, respectively . Subsequently, the outcome of the βNTI analyses in combination with a second null model referred to as a Bray-Curtis-based Raup-Crick (RC bray ) (Chase et al., 2011) model was used to estimate the relative contributions of stochastic processes (Stegen et al., 2013). Specifically, a fraction of pairwise |βNTI| <2 but RC bray <-0.95 or >+0.95 indicated the relative contributions of homogenizing dispersal or dispersal limitation, respectively; a fraction of pairwise |βNTI| <2 and |RC bray | <0.95 indicated that the community assembly was not dominated by any single process .
Linear regressions were performed on Bray-Curtis dissimilarity (dependent variable) versus geographical distance (independent variable) through distance-decay analysis to quantify the distance decay in the bacterioplankton community . The relative contributions of environmental and spatial variables were evaluated by variation partitioning analysis (VPA). First, the spatial variables were derived from geographic coordinates using both linear trend and principal coordinates of neighbor matrices (PCNM) procedures (Griffith and Peres-Neto, 2006) to capture all spatial scales with the PCNM function. Second, the environmental, trend, and PCNM variables that explained a significant (P < 0.05) amount of variation in bacterioplankton community composition were selected by a forward selection model with the forward.sel function (Blanchet et al., 2008). The significance of the selected variables was assessed with 1000 permutations by using the rda and anova.cca functions from the vegan package. Finally, variance partitioning was performed using the varpart function from the vegan package on the Hellinger-transformed bacterial data, the forward-selected environmental variables, and the forward-selected spatial variables, and the data were tested for significance with 999 permutations (Meot et al., 1998). Pairwise comparisons of environmental variables were performed using the corrplot function from the corrplot package. Partial Mantel correlations between environmental variables and bacterioplankton community composition controlled for geographic distance were computed using the mantel.partial function from the corrplot package.

Spatiotemporal Variation and Patterns of Bacterioplankton Taxa and Function
To visualize the overall taxonomic and functional structure of bacterioplankton communities and their underlying driving forces, we used PCoA, Cluster in combination with PERMANOVA on Bray-Curtis dissimilarities among the samples. We found that samples in the overwintering ground display a temporal pattern of divergence along Frontiers in Marine Science | www.frontiersin.org  the first principal coordinate (Figures 4A,B), and this divergence was further corroborated by Cluster, where the communities from November and April formed separate clusters (Figures 4C,D). PERMANOVA of pairwise dissimilarities between bacterioplankton communities indicated that seasonal compartmentalization comprises the largest source of variation ( Table 1). Although the latitude and offshore distances also showed significant effects (Table 1), the community structure did not exhibit a regular or obvious spatial pattern (Figures 4A,B).

Specific Functional Groups and Genera of Bacterioplankton Communities in Two Seasons
In view of the functional differences of the bacterioplankton community between seasons, we further explored the specific functional groups and functional genera for November and April. A model based on the random forest machine-learning algorithm was established to evaluate the importance of indicator  bacterial functions. Ultimately, five bacterial functions were defined as biomarkers when considering the minimum crossvalidation error of the model (Figure 5A). Of these, four bacterial functions showed higher relative abundance in April than November, and only one bacterial function showed higher relative abundance in November than April (one-way ANOVA, p < 0.001) ( Figure 5A). Furthermore, we observed that the dynamics of bacterial genera that drive the biomarker functions showed significant temporal divergence ( Figure 5B). For example, the genus Acinetobacter and Moraxella were abundant in the bacterioplankton communities collected in November, and the genus Rhodococcus and Stenotrophomonas were abundant in the bacterioplankton communities collected in April ( Figure 5B).

Ecological Processes Governing the Spatiotemporal Assembly of Bacterioplankton Communities
To explore the underlying driving forces of bacterioplankton community variation, we used the mean nearest taxon distance (MNTD) to quantify the relative contributions of major ecological processes that regulate the community assembly. The ecological process dominating the community assembly was dispersal limitation (accounting for 67.0-79.7%), which was relatively stable over the seasons (Figure 6). Second, 20.3-31.0% of the variation in community composition was primarily due to heterogeneous selection (Figure 6). In addition, there was still 0-3.7% of compositional variation unexplained by any ecological processes (Figure 6). These results indicate that stochastic rather than deterministic processes structured the bacterioplankton community assembly in the overwintering ground.

Distance-Decay of Similarities in Taxonomy and Function of Bacterioplankton Communities
The spatial dissimilarities in bacterioplankton taxonomic composition were significantly higher than that in functional composition ( Figure 7B). Furthermore, the taxonomic composition of bacterioplankton communities collected in both November (adjusted R 2 = 0.139, p < 0.001) and April (adjusted R 2 = 0.040, p < 0.001) exhibit a clear distance-decay pattern in similarities ( Figure 7A). This pattern was also FIGURE 6 | Summary of ecological processes governing the compositional variation of total bacterioplankton communities and the bacterioplankton communities collected in November and April. The percentages are given as the relative contribution of each process to the community variation as indicated by different colors. observed in the bacterial functional composition in November (adjusted R 2 = 0.142, p < 0.001) but not in April (adjusted R 2 = -0.001, p = 0.395) ( Figure 7A).

Drivers of Bacterioplankton Taxonomic and Functional Composition
To explore the deterministic mechanisms underlying community assembly, we used VPA for disentangling the relative role of environmental and spatial variables. Results of the variation partitioning showed that the drivers of bacterioplankton community composition were different between seasons and between taxonomic and functional levels (Figure 8). Specifically, the environment-related variables (pure environmental variables, environmental-trend covariation, environmental-PCNM covariation, covariation of all variables) dominated the variation of bacterioplankton taxonomic composition collected in all seasons (both November and April) and November (only November) and the variation of bacterioplankton functional composition collected in April (only April), while space-related variables (pure linear trend, pure PCNM, environmental-trend covariation, environmental-PCNM covariation, trend-PCNM covariation, covariation of all variables) dominated the variation of bacterioplankton taxonomic composition collected in April and the variation of bacterioplankton functional composition collected in all seasons and November (Figure 8). Remarkably, there was a large proportion (53.4-76.7%) of variation unexplained by the above variables. To further identify the geographic-independent environmental drivers, we correlated distance-corrected dissimilarities of taxonomic and functional community composition with those of environmental variables (Figure 9). Overall, nutrients were the strongest correlates of both taxonomic and functional composition in both seasons. Moreover, significant correlations were found between the temporal differences in community composition (between the seasons) and various environmental variables (Figure 9A), FIGURE 8 | Variation partitioning of the taxonomic and functional structure of total bacterioplankton communities (all) and the bacterioplankton communities collected in November (nov) and April (apr) by the seawater environmental variables (Env) and the spatial variables including linear trend (Trend) and PCNM variables (PCNM). and the spatial differences in community composition (within the season) were weakly driven by environmental variables (Figures 9B,C).

The Spatiotemporal Patterns and Characteristics of Surface Bacterioplankton Communities in an Overwintering Ground
The habitats of the large yellow croaker are being degraded due to the increased use of coastal areas for aquaculture operations and development (Kang et al., 2018). The failed recovery of the common snook (Centropomus undecimalis) indicates that habitat degradation is considered more significant than exploitation in explaining population declines (Brennan et al., 2008), which indicates the need for the study and restoration of croaker habitat. Planktonic bacteria are generally used as a bioindicator to systematically reflect or record the conditions of the marine environment or fish habitats due to their biogeochemical significance and direct relationship with fish health (Dowle et al., 2015;Dai et al., 2017). Over the past decade, dozens of studies have characterized patterns, associations, and drivers in the bacterioplankton community in various fish habitats and elucidated similarities and differences in the bacterioplankton communities between fish farms and natural areas (Tamminen et al., 2011;Xiong et al., 2015;Jing et al., 2019). However, comprehensive information on the spatial and temporal patterns of bacterioplankton communities in the croaker habitat is scarce.
Our study systematically investigated the surface bacterioplankton communities in an overwintering ground for large yellow croakers. We found that, during the overwintering period, the dominant bacterioplankton phyla and functions remained unchanged, but their relative abundance exhibited seasonality (Figures 2, 4 and Table 1). Specifically, the bacterioplankton community in autumn was characterized by animal parasites or symbionts (Figure 5A), and Acinetobacter was the representative genus ( Figure 5B). Acinetobacter baumannii, which was detected in our samples, is a virulent pathogen for many fish species that can trigger serious septicemia outbreaks (Brahmi et al., 2016;Li et al., 2017). More seriously, A. baumannii has been proven to be a serious human clinical pathogen involved in a large number of infections, including bacteremia, secondary meningitis, urinary tract infection, and ventilator-associated pneumonia (Mari-Almiralla et al., 2019). Therefore, the presence of Acinetobacter in overwintering grounds increases the latent risk of fish disease and the possibility of coinfection between fish and humans through handling or consumption (Xia et al., 2008). In addition, other functional genera found in autumn, such as Moraxella (Figure 5B), are opportunistic pathogens of fish (Walke et al., 2015). Once the environment deteriorates or suddenly changes, those bacteria might threaten the health of the croaker or other fish. In spring, many nutrient-related functional bacterioplankton were dominant (Figure 5). Rhodococcus, one of the dominant functional genera in April (Figure 5B), is generally regarded as the most versatile and efficient group of organisms that are suitable for the biodegradation of organic compounds (Larkin et al., 2005;Martínková et al., 2009). Given this characteristic, Rhodococcus has been found worldwide in FIGURE 9 | Environmental drivers of taxonomic and functional composition of total bacterioplankton communities (A) and the bacterioplankton communities collected in November (B) and April (C). Pairwise comparisons of environmental variables are shown with a color gradient denoting Spearman's correlation coefficients. The taxonomic and functional composition of bacterioplankton communities was related to each environment variable by partial (geographic distance-corrected) Mantel tests. Edge width corresponds to Mantel's r statistic for the corresponding distance correlations, and edge color denotes the statistical significance based on 9999 permutations. WT, water temperature; SAL, salinity; PO 4 , orthophosphate; NO 2 , nitrite; NO 3 , nitrate; NH 4 , ammonium; SiO 3 , silicate; TN, total nitrogen; TP, total phosphorus. eutrophic aquatic ecosystems (Wu et al., 2009;Lee et al., 2010), suggesting that it might be used as an organismal indicator of nutrient enrichment (Rowbotham and Cross, 1977). Moreover, Alcanivorax (Figure 5B), a hydrocarbon-degrading bacterium, could also be enriched in eutrophic environments and could be related to the presence of pollutants (Patel et al., 2014). Other functional genera in spring, such as Stenotrophomonas (Jia et al., 2019) and Neptunomonas (Lanfranconi et al., 2010), were also organic-degrading or nitrogen-related organisms that prefer organic-rich or eutrophic conditions. Therefore, the functional composition of the bacterioplankton community in April suggests that the nutrient level in the surface water of the habitat in spring was higher than that in the winter, which was verified by environmental data (Supplementary Table S1) and implies that this sea area has a high risk of developing eutrophication.

The Mechanism of Surface Bacterioplankton Community Assembly in the Overwintering Ground for Large Yellow Croaker
One key question in ocean microbial ecology is the extent to which environmental filtering and various biological interactions, on the one hand, and limited dispersal and historical contingency, on the other, are responsible for spatiotemporal patterns (Sunagawa et al., 2015). Nevertheless, the answer to this question in marine ecosystems is still controversial. For example, some studies have proposed that marine bacterioplankton communities are structured by deterministic factors (Pontarp et al., 2012;Sunagawa et al., 2015;Wang et al., 2019), and other studies have demonstrated that stochastic factors can potentially overwhelm the influences of selection on community structure (Beaton et al., 2016;Evans et al., 2017;Zhou and Ning, 2017). The latter theory is supported by our findings as dispersal limitation dominated the processes governing both the spatial and temporal scales of bacterioplankton community assembly (Figure 6). It has historically been considered that free-living bacteria are small, numerous, and easily transported and, hence, do not experience passive dispersal limitations (Cottenie, 2005;Shurin et al., 2009). However, numerous studies have recently documented that bacterial communities in water show spatial differences and reveal strong biogeographic patterns (Verreydt et al., 2012;Müller et al., 2014;Wu et al., 2018). The distance-similarity decay relationships for the bacterioplankton communities observed in our study (Figure 7) also suggest the impact of dispersal limitation on the generation of biogeographical patterns in marine bacteria (Bell, 2010;Adams et al., 2013). It is interesting to note that, even though the study area is located in the East China Sea and the geographic scale of the overwintering ground we investigated is comparatively broad, the distancedecay patterns in our study were weaker than those in other studies as the distance decay rates (the regression slopes) of our bacterioplankton communities were lower Dai et al., 2017). The strong distance-decay patterns in Wang et al. (2015) and Dai et al. (2017) were mainly attributed to spatially structured environmental gradients, which essentially suggests that the spatial heterogeneity of bacterioplankton communities in their studies was structured by the deterministic heterogeneous selection of abiotic and biotic environmental conditions (Dini-Andreote et al., 2015). However, stochastic dispersal exerted a dominant role in shaping our bacterioplankton communities (Figure 6), contributed to the relatively random distribution of the bacterial community compared with that under harsh environmental filters (Evans et al., 2017), and attenuated the heterogeneity of bacterial taxa in our study. Additionally, lower or no distance-similarity decay relationships were found in bacterial functions compared with the relationships found in bacterial taxa (Figure 7). This is a universal phenomenon because functional redundancy, which means that different populations share similar or the same functions, appears to be quite high in the microbial community, especially compared with the level of functional redundancy in plant or animal communities (Delgado-Baquerizo et al., 2016;Louca et al., 2018).
Comparatively, environmental stressors wielded less power than stochastic processes in creating temporal and spatial differences in the bacterioplankton communities' taxonomical and functional composition in the overwintering habitat for croaker. Our variation partitioning results showed that over 50% of the variation in the bacterioplankton community structure was unexplained by the environment-related variables (Figure 8). The lower power of the environmental variables could be explained by the restraining force of stochastic effects (Zhou and Ning, 2017); dispersal limitation can obscure the relationships between environmental variables and microbial community composition, resulting in apparently weak deterministic effects (Evans et al., 2017).
Among the purely environmental filters that were independent of spatial effects, nutrient-related variables, such as silicate, nitrate, and ammonium, could have exhibited some degree of heterogeneous selection on temporal and spatial variations in the bacterioplankton community composition in the overwintering grounds (Figure 9). The mechanism by which silicate influences the composition of the bacterioplankton community could be related to its influence on diatom growth; it may shape bacterioplankton community composition due to the tight coupling between diatoms and bacteria (Topper et al., 2010;Wang et al., 2015). The correlations of bacterioplankton community composition with nitrogenous nutrients could be attributed to the abundant biodegrading bacteria in the habitat (Figure 5), which are capable of decomposing organic substances into various inorganic nutrients (Manser et al., 2006;Li and Boyd, 2016). An equally important driver of variation in the bacterioplankton community was temperature (Figure 9); this finding is in accordance with those reported in freshwater ecosystems, oceans, and forests at the global scale (Penuelas et al., 2012;Sunagawa et al., 2015;Chen C. Y. et al., 2019). The importance of salinity in shaping the composition of the bacterioplankton community was also observed (Figure 9). Salinity may contribute to density gradients that physically separate water masses and their resident bacterial communities from near-shore to offshore areas (Lozupone and Knight, 2007;Fortunato et al., 2012) (Supplementary Table S1). Notably, more correlations were found between the temporal differences in community composition and the various environmental variables than between spatial differences and environmental variables (Figure 9); this difference may have resulted from the greater differences in environmental variables between seasons. Interestingly, the interactions among community composition and purely environmental filters were more complex and closer in spring than in autumn (Figures 9B,C). One possible explanation for this difference is that the high eutrophication level in April (Supplementary Table S1), as a long-term and continuous disturbance, could increase the importance of deterministic environmental filtering to bacterioplankton community assembly (Figure 6) (Dai et al., 2017).
Overall, we obtained a unique taxonomic and functional profile of the bacterioplankton community in the surface water of the overwintering ground and revealed the assembly mechanisms that structure it on the spatial-temporal scale. The systematic map of the surface bacterioplankton community has implied that current conditions are not conducive to the high primary productivity and the health of upper elements of the croakers' food chain. There are, however, potential limitations to our approach that merit further discussion. First, although the above results and discussion display the patterns and processes of bacterioplankton community assembly in the overwintering grounds for large yellow croaker, samples collected from the surface water can only somewhat reflect habitat conditions for the croaker who live at moderate and lower water depths (Lü et al., 2008). Therefore, we will further study the habitat microbiome of the water masses in which croaker is located. Additionally, we did not sample replicates at the study sites, which might have limited the analysis of the within-site variability and thereby might affect the accuracy of the results of the random forest classification and the distance decay analysis to a certain extent.

DATA AVAILABILITY STATEMENT
All the raw sequences have been uploaded to the NCBI Sequence Read Archive (SRA) database under the BioProject number PRJNA604077 and the accession number SRP247264.

AUTHOR CONTRIBUTIONS
WY and Z-MZ led manuscript writing and data analysis. S-ZZ, S-HZ, and BL contributed with sequence analysis, statistical analyses, and microbiological interpretations. LZ and RN collected and processed samples. J-YZ and C-HL were responsible for organizing field sampling and inter-institution communication. All authors contributed to manuscript revision, read, and approved the submitted version.