1 INTRODUCTION

In recent years, there has been an increasing emphasis on the conservation of fish spawning habitats. Therefore, research on the protection of fish spawning habitat is crucial, and it has potential implications in fisheries resource restoration and environmental protection (Yan and Chen, 1998; Zhang and Yang, 2012; Chen et al., 2015). The large yellow croaker Larimichthys crocea is one of the most important commercial fish species in China. Therefore, exploring its spawning ground for the improvement of its ecological restoration is very important. However, there is an urgent need for quick monitoring and management initiatives to improve the spawning grounds for large yellow croaker in Sansha bay. Previous studies on the evaluation of restocking measures in the Sansha bay have mainly focused on the population dynamic of fish (Shen, 2011). However, studies on spatial-temporal dynamics of bacterioplankton communities in the inner-bay large yellow croaker habitat of Sansha Bay have not been sufficiently researched.

The nitrogen cycle, sulfur cycle, and carbon cycle in the marine ecosystems are inseparable from marine bacterioplankton. Most of marine bacterioplankton are decomposers, and a small part of them involves in the decomposition and transformation of organic compounds as producers. The degradation of toxic and harmful organics by bacterioplankton has greatly improved the self-purification ability of the marine ecosystem (Enns et al., 2012). Given the sensitivity of bacterioplankton to water quality (Sunagawa et al., 2015) and its special position in the marine geochemical circulation (Brown et al., 2014), there is a need to evaluate the health level of the marine ecosystem through the relationship between bacterioplankton and the marine water environment.

In addition, bacterioplankton also played an important role as potential health-beneficial bacteria and pathogenic bacteria of marine organisms (Xiong et al., 2015), and their role in the "microbia food web" (Fuhrman, 1992). Moreover, dissolved organic carbon (DOC) synthesized by heterotrophic bacteria enters a higher trophic level through the predation of zooplankton on bacterioplankton (Li et al., 2019), thus becoming an important part of the "microbial food web" of marine organisms. Currently, there have been studies on bacterioplankton as potential health indicators for a variety of marine organisms, including but not limited to fish, shrimp, and jellyfish (Qu et al., 2014; Tapia-Paniagua et al., 2014; Zhang et al., 2014).

Sansha Bay possesses abundant baits and is an appropriate place of reproduction and breeding for fish due to the mixing of sea and fresh water, among which the Guanjing Ocean is the only inner breeding and spawning ground of large yellow croaker in China (Xu and Xu, 2013). In the studies on the feeding habits of large yellow croakers at various growth stages, zooplankton, mainly copepods, was the main food source for juveniles (Zheng and Yang, 1965), and the abundance and diversity of zooplankton in Sansha Bay showed obvious seasonal changes. Therefore, given the relationship between zooplankton and bacterioplankton, studies on bacterioplankton community structure and seasonal succession in Sansha Bay might have great benefit for the health status of large yellow croaker habitats.

Although there are many studies on different sea areas of Sansha Bay (Lin et al., 2019; Xie et al., 2020), there are no enough researches on community structure and seasonal succession of bacterioplankton in the area. This article attempts to study the differences in the environmental parameters and bacterioplankton communities among the traditional breeding areas of large yellow croakers in Sansha Bay. Also, this study intends to explore the mechanism of seasonal succession of bacterioplankton to discuss the potential impact of bacterioplankton on large yellow croakers with the constant development of aquaculture industry and fisheries sector. The research aims to examine the following aspects: (1) the annual changes of bacterioplankton community structure in different areas of Sansha Bay; (2) drivers of change in bacterioplankton community structure in Sansha Bay across seasons and sea areas; and (3) the potential influence of bacterioplankton on the reproduction and growth of large yellow croaker in Sansha Bay.

2 MATERIAL AND METHOD 2.1 Study area and sampling strategy

Sansha Bay is a typical inner bay, including Yantian Harbor, Baima Harbor, Guanjing Ocean, and Dongwu Ocean (Shen, 2011), located in the northeast of Fujian Province of China. Water samples were collected in January, April, July, and October 2019 from the Sansha Bay at a depth of 0.5 m, respectively (Fig. 1). Except for the winter in January 2019, which only investigated the waters inside Sansha Bay and collected 21 samples, the other three seasonal further optimized the survey sites and investigated both the internal and external waters of Sansha Bay, collecting 28 samples each. A total of 105 samples were collected from four seasonal voyages, all of which were analyzed in this study. The sampling sites were divided into three groups based on historical geographical factors, namely, Yantian Harbor (YH), Breeding Area (BA), and Bay Margin (BM). YH includes 5 sites: S01, S02, S03, S04, and S05. BA included both Guanjing Ocean and Fuding Ocean, which have 16 sites: S06, S07, S08, S09, S10, S11, S12, S13, S14, S15, S19, S20, S21, S22, S23, and SF. BM included 7 sites: S24, S25, S26, S27, S28, S29, and S30.

2.2 Collection of bacterioplankton and determination of environmental parameters

An amount of 1 500 mL of water samples were collected and filtered through a 100-μm silk yarn on the sampling day, and then the bacterioplankton were filtered to three 0.22-μm polycarbonate membranes (47-mm diameter, Millipore, Boston, MA, USA), one of which was used for subsequent analysis, two for backup, and each membrane filtered 500 mL of water. All samples were stored in an ultra-cold refrigerator at -80 ℃ until further analysis. The DNA isolation was performed within a week of sampling.

Multi 3630 IDS (WTW, Germany) was used to measure parameters at each station during sampling, including dissolved oxygen (DO), pH, and water temperature (T). Nitrate (NO₃^-), nitrite (NO₂^-), active phosphate (PO₄^3-), ammonia nitrogen (NH₄⁺), total nitrogen (TN), and total phosphorus (TP) was determined using an automated chemical analyzer (Smart-Chem 200 Discrete Analyzer, Westco Scientific Instruments, Brookfield, USA), all referred to China Marine Monitoring Code (General Administration of Quality Supervision, Inspection and Quarantine of the People's Republic of China and Standardization Administration, 2008). The values of the relevant environmental factors were available in Supplementary Table S1.

2.3 DNA extraction and sequencing

Bacterioplankton DNA was extracted using a DNA extraction kit (Minkgene Water DNA Kit), then concentrated, and measured DNA purity by NanoDrop One spectrophotometer (Thermo Fisher Scientific, MA, USA). The Ⅴ4 region of the 16S rRNA gene was amplified by universal primers 515F (5ʹ-GTGCCAGCMGCCGCGGTAA-3ʹ) and 806R (5ʹ-GGACTACNNGGGTATCTAAT-3ʹ) for bacterioplankton communities (Caporaso et al., 2011). Illumina HiSeq 2500 platform was used to generate a 250-bp double-terminal paired sequence (Guangdong Magigene Biotechnology Co., Ltd., Guangzhou, China).

2.4 Bioinformatics

Bioinformatics processing of sequencing data was performed by USEARCH (Ⅴ11.0.667_I8). The UNOISE3 algorithm was used to de-noise the combined sequence (unoise_alpha=2, minsize=4, according to the default settings), and the errors were corrected and chimerism was removed to generate zero-radius Operation Taxon Units (ZOTUs) and quantified ZOTUs (Edgar, 2016). The SILVA (Ⅴ 138) was used to assign 99% similarity representative sequences to each ZOTU using the RDP classifier.

Functional Annotation of Prokaryotic Taxa (FAPROTAX) was used to predict the function of bacterioplankton (Louca et al., 2016).

2.5 Statistics

All data pre-processing was done by Excel. All visualizations were done with the help of the package "ggplot2" in R 4.0.2, except for the sampling point map (Fig. 1), which was done with ArcGIS 10.6.

PERMANOVA was performed by the function adonis() in the package "vegan". analysis of variance (ANOVA) was performed by the function aov() in the package "vegan". analysis of similarity (ANOSIM) was performed by the function anosim() in the package "vegan". principal coordinate analysis (PCoA) was performed by the function cmdscale() in package "ape" and the function ddply() in package "plyr". The indicator ZOTUs among different seasons were identified using the function multipatt() in package "indicspecies" and likelihood ratio tests using the function DGEList() in the package "edgeR". The bioindicator functions were generated by the function randomForest(), and were identified by the most important season-discriminating features using the function rfcv() in package "randomForest". redundancy analysis (RDA) was performed by the function rda(), and variance inflation factors (VIFs) were calculated by the function vif. cca(), environmental factors were selected by a forward-selection procedure using the function ordiR2step() in the package "vegan".

3 RESULT 3.1 Variance of α-diversity of bacterioplankton

A total of 3 640 350 high-quality sequences were acquired after analysis and were further divided into 11 953 ZOTUs across all samples, and each sample contained 34 670 high-quality sequences. Chao1, Shannon, richness, and Simpson indices were selected as reference standards for the α-diversity of bacterioplankton community in Sansha Bay. All indices of BA were significantly lower than YH in winter and higher than YH in spring, except for the Simpson index. The indices of different sea areas in summer and autumn tended to be similar, except for the Simpson index, which was significantly lower in BA than that in the other two sea areas (Fig. 2).

3.2 Seasonal succession of bacterioplankton community composition

There were significant differences in the annual bacterioplankton community of different areas in Sansha Bay (PERMANOVA, P < 0.001). The differences were visualized by cluster analysis and PCoA based on the Bray-Curtis dissimilarity distance matrix (Fig. 3). The division of the bacterioplankton community was more obvious in winter and spring, and some loci overlapped in summer and autumn, but there were significant differences among different seasons (Supplementary Table S2). In the same season, there was a certain degree of difference in the community structure of bacterioplankton in different sampling areas (Supplementary Table S3). The difference of bacterioplankton community between BM and the other two areas was significant (PERMANOVA, P < 0.001), and there was also a significant difference between YH and BA (PERMANOVA, P < 0.01).

At the level of genera, there were differences in the relative abundance of the main bacterioplankton in Sansha Bay in different seasons (Fig. 4). Except in autumn, the relative abundance of Candidatus Actinomarina and Candidatus Aquiluna was relatively high in Sansha Bay. The relative abundance of Candidatus Aquiluna was low in autumn, while the relative abundance of Candidatus Actinomarina was significantly higher in that season than in other seasons (P < 0.001). Compared with the other three seasons, the relative abundance of AEGEAN-169_marine_group in spring was significantly lower (one-way ANOVA, P < 0.05), while the higher relative abundance of NS5_marine_group, Pseudoalteromonas, and Planktomarina occurred in that season (one-way ANOVA, P < 0.05).

3.3 Unique bacterioplankton community in different seasons

A total of ZOTUs with a relative abundance of > 0.01% were used to analyze the unique seasonal bacterioplankton community in different areas of Sansha Bay. A total of 669 ZOTUs with relative abundance greater than 0.01% were found in BA (Fig. 5a), of which 595 ZOTUs (88.94%) showed significant differences among different seasons (P < 0.001). Although differential ZOTUs screening was also carried out in the other two sea areas, the seasonal differences were not as obvious as those in BA. No seasonal differences were detected in 35.89% of ZOTUs in YH (Fig. 5b), while in BM this percentage reached 62.41% (Fig. 5c).

There was a seasonal succession of unique bacterioplankton community in different areas of Sansha Bay (Fig. 5d–f). For example, the relative abundance of Actinobacteria increased from January to October in BA and YH, while it showed the opposite trend from April to October in BM. Meanwhile, the relative abundance of Bacteroidetes in YH and BA was not significantly higher than in the three seasons of BM. It is worth noting that a significantly higher relative abundance of Cyanobacteria was found in YH in July and BM in October, while a non-significantly higher relative abundance was observed in each season of BA.

3.4 Bioindicator functions and corresponding genera

To explore the functional differences of bacterioplankton communities in different seasons and areas of Sansha Bay, five functional groups were selected as bioindicators based on a random forest machine-learning algorithm (Fig. 6a). Among them, the relative abundance of four functions was significantly higher in summer than in other seasons, and only the function aerobic_nitrite_ oxidation was significantly higher in autumn (one-way ANOVA, P < 0.05). Different bioindicator functions corresponded to different bacterioplankton (Table 1). In addition, temporal and spatial divergence was observed in the bacterioplankton genera corresponding to bioindicator functions, such as the relatively higher abundance of Rhodococcus in summer, especially in BM (Fig. 6b).

3.5 Drivers of bacterioplankton community and function

To specify the factors that drive the differences in the bacterial community functions of Sansha Bay, RDA and correlation analysis were used to visualize the relationship between environmental factors and bacterioplankton (Fig. 7). For the bacterioplankton with a relative abundance greater than 0.01%, most environmental factors were significantly correlated with them. For example, water temperature (T) was significantly positively correlated with summer and autumn, while DO did not follow a similar pattern (Fig. 7a). Furthermore, at the genus level, most environmental factors were also significantly correlated with bacterioplankton genera corresponding to different bioindicator functions (Fig. 7b).

4 DISCUSSION 4.1 Seasonal succession of bacterioplankton in Sansha Bay and its relationship with the environment

Numerous previous studies have shown that marine bacterioplankton community structure has significant interannual variability, with temperature and nutrient salinity being important factors influencing the structure of bacterioplankton communities in local waters (Fuhrman and Steele, 2008; Kataoka et al., 2009). Moreover, the study conducted by Yang et al. (2020) on the bacterioplankton in the overwintering area of large yellow croaker, indicated that there was a significant difference on the structure of bacterioplankton communities between autumn and spring. In this study, significant seasonal changes were observed in both environmental parameters and bacterioplankton community composition in Sansha Bay (Fig. 2; Supplementary Table S1), and there was some association between them.

Nutrient levels in Sansha Bay were high on an annual basis, especially in summer and autumn, which was related to the net-pen culture of large yellow croaker in Sansha Bay. Large yellow croaker tended to be harvested before the Chinese Lunar New Year (in February) and put into the next batch in spring, so the amount of feed used in spring, summer, and autumn was much greater than in winter (Xie et al., 2020), which explained the significantly higher nutrient concentration in summer and autumn. However, the nutrients in Sansha Bay did not show a significant increase in spring, which may be related to the lack of precipitation in Ningde spring. Thus, it can be seen that although aquaculture has negatively affected the water of Sansha Bay to some extent, the main factor contributing to the increase in nutrient salt concentration in Sansha Bay is still land-based pollution.

Moreover, at the same time, no large-scale algal bloom was observed in the all four seasons, which was relatively rare in the other high nutrient concentration water of Ningde (Xie, 2018), and large seaweed farms near the sampling area may explain this phenomenon. Large seaweeds provided a good buffer for nutrient salinity in Sansha Bay, reducing the ecological risk of algal blooms (Xie et al., 2020).

The results of a related study in Chesapeake Bay, USA, showed that there were significant differences in the structure of bacterioplankton communities caused by cold and warm seasons, with a large number of seasonal transitional taxa (Kan et al., 2006). In this case, the spatial and temporal differences in bacterioplankton communities in Sansha Bay may be originated from changes in water temperature. In Sansha Bay, where high temperatures prevail in summer and autumn and lower water temperatures occur in winter, Psychrobacter (Fig. 5d–e) was found in the main ZOTUs specific to both YH and BA in winter, which can effectively adapt to low-temperature environments and was found in previous studies to be of high abundance in hydrocarbon-contaminated waters (Prabagaran et al., 2007). Meanwhile, the genus Arcobacter, most of which were also found to be adapted to low temperatures (Canion et al., 2013), had a significantly higher relative abundance at both BA sites in winter (Fig. 3), and most of them were pathogenic to both humans and animals (Ho et al., 2006). Cyanobacteria (Fig. 5e–f) were found to have very significantly higher relative abundance in the major ZOTUs of YH in summer and BM in autumn, and were mostly dominated by Synechococcus. In previous studies, Synechococcus was found to be more adapted to warmer environments, and a certain degree of high nutrient salinity promoted its reproduction (Partensky et al., 1999). Our study also showed that the two groups with the highest relative abundance of Cyanobacteria were also the seasons with the highest water temperature and nutrient salinity in the same sea area (Fig. 5). However, in BM, where significantly higher relative Cyanobacterial abundance was observed in autumn, water temperature and nutrient salinity levels were still significantly lower than in the other two sea areas during the same period. Nevertheless, the nutrient salinity of BM was still higher in autumn compared to other seasons (Supplementary Table S1), and this was corroborated by the fact that wild large yellow croaker face migration to offshore wintering in autumn (Xu and Chen, 2011). The abundance of Cyanobacteria in BM autumn may provide a large source of carbon to provide abundant bait for the large yellow croaker population through the microbia food web.

On the other hand, the variability of eutrophication degree in Sansha Bay also interacted with the changes in the community of bacterioplankton. For example, the relative abundance of Bacteroidetes, which was relatively higher in both BA and YH, was significantly lower in BM (Fig. 5). In this study, the Bacteroidetes consisted mainly of NS5_marine_ group and Owenweeksia, both of which were found to have chemoheterotrophic ecological functions (Gómez-Pereira et al., 2010), and this was corroborated by the significantly lower nutrient salt levels in BM.

4.2 Potential effects of bacterioplankton community on large yellow croaker

Given the diversity of potential probiotic and pathogenic bacteria in marine bacterioplankton, different bacterioplankton community compositions could have an impact on the health of marine organisms (Zhang et al., 2014; Duarte et al., 2019). In the studies available so far, bacterioplankton involving marine organisms in breeding areas often differ significantly or to some extent in the distribution of potentially probiotic and pathogenic bacteria (Hu et al., 2015; Jing et al., 2019). For example, in a study of bacterioplankton in the cultured and non-cultured areas of Xiangshan Port, potentially pathogenic bacteria were found to be significantly higher in relative abundance in the non-cultured areas, while probiotic bacteria with antimicrobial activity were significantly higher in relative abundance in the cultured areas (Hu et al., 2015). In contrast, the difference of bacterioplankton in Sansha Bay in this study is mainly manifested in season, which may have an indirect effect on the reproduction and growth activities of large yellow croaker (Shen, 2011).

Spring is the peak of large yellow croaker spawning (Liu, 2004), and the large yellow croaker in Sansha Bay mainly carries out activities in BA dominated by the Guanjing Ocean. The relative abundance of Rhodobacteracae was significantly higher in spring BA compared to the other three seasons, and Rhodobacteracae contained many bacterioplankton with extensive antibacterial effects (Murray et al., 2011); while the relative abundance of Arcobacter, Pseudoalteromonas, and Vibro in the other three seasons was significantly higher than that in spring (Fig. 5), and they contain a large number of potentially pathogenic bacteria (Porsby et al., 2008; Prado et al., 2009; Soelberg et al., 2020). A relatively high abundance of Rickettsiales was observed in YH and BM in spring, and it has been extensively studied as a common pathogen (Fryer et al., 1990). Thus, at the level of bacterioplankton, BA in spring may provide a relatively suitable microecology environment for the spawning large yellow croaker. The speculated reason for this difference is the high abundance of Rhodobacteraceae which may inhibit the reproduction of a variety of potentially pathogenic bacteria, but the specific mechanisms of their interactions need to be further investigated.

FAPROTAX was used to further assess the ecological role of the bacterioplankton community in Sansha Bay on the growth and reproduction process of large yellow croaker in terms of elemental cycling. Similar to the results of previous studies in marine environments (Huang, 2018; Yang et al., 2020), chemoheterotrophy and aerobic_ chemoheterotrophy associated with C-cycling dominated the functional annotation results in Sansha Bay (Supplementary Fig.S1). As an organic-rich inner bay breeding area, Sansha Bay contained a large amount of aerobic heterotrophic bacterioplankton such as Owenweeksia and Pseudoalteromonas (Fig. 4), which performed those functions. Although the proportion of the above two types of functions was the highest in Sansha Bay, the ZOTUs (corresponding to these two types of functions) were relatively abundant in different seasons and sea areas (Supplementary Fig.S1).

The top five bioindicator functions from random forest machine learning included aerobic_nitrite_ oxidation (Fig. 6a), which is also involved in the N-cycle, in addition to compound degradation. Meanwhile, similar to many studies (Lücker et al., 2013; Levipan et al., 2014), the aerobic_nitrite_ oxidation-related bacterioplankton genus Nitrospina was the only genus among the corresponding genera that was positively correlated with NO₂^- (Fig. 7b), while it presented a higher relative abundance in the autumn when NO₂^- concentrations were higher (Supplementary Table S1). Among them, the relative abundance of these functions was always at a lower level in winter, which, according to the previous study (Zhou et al., 2021), maybe due to the generally lower nutrient salt concentration.

The bioindicator functions were not only found in significantly different relative abundances in different seasons (Fig. 6), but their corresponding bacterioplankton genera also differed significantly in different seas in the same season (Fig. 6b), indicating the different functions of different areas of Sansha Bay in the marine elemental cycle. It is noteworthy that multiple bioindicators functionally related bacterioplankton genera all showed the highest relative abundance in BM for that season (Fig. 6b), but there was no significantly higher concentration of nutrients in the BM compared to the other two marine areas (Supplementary Table S1). Considering that, BM is the only waterway for water exchanging between Sansha Bay and the offshore seas, the high presence of functionally related bacterioplankton in the BM may provide potential ecological supports for the growth and reproduction of large yellow croaker in the BA.

5 CONCLUSION

This study revealed the seasonal succession of the structure and function of the bacterioplankton community and its driving factors in different sea areas of Sansha Bay. The results showed that bacterioplankton community structure in the Breeding Area reflected a more suitable ecological environment in terms of composition than that in Yantian Harbor and Bay Margin. Meanwhile both the Yantian Harbor and Bay Margin functioned as baiting areas in summer and autumn, and many environmental factors were significantly associated with them. In addition, the results of functional prediction also differed among seasons and sea areas and were closely related to environmental changes. These findings suggest that exploring the seasonal succession of bacterioplankton community structure in different functional sea areas in Sansha Bay will help in understanding the ecological status of the traditional breeding area of large yellow croaker at the microbial level, and thus for better ecological resource restoration.

6 DATA AVAILABILITY STATEMENT

The raw paired-end sequences have been deposited in the NCBI Sequence Read Archive (SRA) database under BioProject number PRJNA807812 and accession number SRP360435.

7 ACKNOWLEDGMENT

The authors would like to thank the colleagues from Xiamen University, Second Institute of Oceanography, and other institutions for helping onboard. In addition, the authors thank Prof. Lingfeng HUANG and Prof. Jianyu HU from Xiamen University for their selfless sharing of environmental data.

Electronic supplementary material

Supplementary material (Supplementary Tables S1–S3 and Fig.S1) is available in the online version of this article at https://doi.org/10.1007/s00343-022-1431-8.