Temporal dynamics of zooplankton community in an oyster farming area of the Yellow Sea in China via metabarcoding

Introduction

Zooplankton functions as a bridge between producers and consumers. The changes in the zooplankton abundance affect primary and higher-order productivity in the oceans (Hernández-MIranda et al., 2022). The study of zooplankton composition is of great value to support research on ecosystem dynamics (Shi et al., 2020) and biofouling management. Many zooplankton species have short life cycles, and their richness and distribution are greatly affected by environmental factors. Community structure can serve as an indicator of local ecological conditions and environmental health and provide basic data for ecosystem detection (Pearson, 1978; Zhao et al., 2022). Oysters are the predominantly cultured species in China, potentially changing the composition of zooplanktonic and biofouling organisms around the oyster farm (Das et al., 2012). However, the structure and dynamic patterns of the zooplankton community in the oyster farm are poorly studied.

DNA barcoding technology has become increasingly important in ecological research (Valentini et al., 2009; Mugnai et al., 2021). Compared with traditional methods, it has the advantages of fast identification and high accuracy (Savolainen et al., 2005). The 18S rRNA gene includes conserved and highly variable regions, which can reflect genetic relationships among species and show the differences between species. It has been successfully applied in many studies: the discovery of new species (Kakui and Munakata, 2022; Kumar et al., 2022), the assessment of marine biological diversity (Liu et al., 2020), community characteristics and environmental correlation analysis (Liu et al., 2022). This study analyzed the species richness, diversity, and temporal dynamics of zooplankton in an oyster farm area using the third-generation sequencing-based DNA barcoding technique. The correlation analysis indicated significant interactions between zooplankton species and provided a basis for in-depth research on the interactions between human activities and marine ecosystems.

Materials and methods

Sampling sites and sample collection

The study was carried out from May 2020 to May 2021 in the Pacific oyster farm located in Longwan Bay, the northern part of the Yellow Sea (Figure 1A). A total of 33 monthly zooplankton samples were collected from three sites by a shallow water type- II plankton net (∼50 μm mesh size) at a depth of ~1 m where the Pacific oysters were cultured. Samples were fixed in 70% ethanol at 4°C for further processing. Water temperature, salinity, pH, and dissolved oxygen (DO) were measured in situ with a multi-parameter portable meter (Multi 3630 IDS, Germany).

Figure 1 Zooplankton richness and relative abundances. (A) Locations of 3 sampling sites (point A, point B and C) in the aquaculture area. (B, C) Operational taxonomic units (OTUs) richness and reads numbers of zooplankton phyla and classes. The pie chart represents richness. The bar chart represents abundance (the number of reads). (D) Reads number and the percentage distribution of reads per month changes for Top30 species (ranking based on reads number).

Genomic DNA extraction, PCR amplification and high-throughput sequencing

Samples were pretreated by centrifugation at 5000 ×g, and the pellets were washed twice with PBS. The total genomic DNA was then extracted using the Tissue DNA Kit (Omega, USA) according to the manufacturer’s protocol. DNA quality and quantity were examined by 1% agarose gel electrophoresis and NanoDrop 2000 (Thermo Fisher Scientific, USA). The 18S rRNA genes were amplified using universal primers: the forward primer 5’-AACCTGGTTGATCCTGCCAGT-3’, and the reverse primer 5 ‘-GATCCTTCTGCAGGTTCACCTAC-3’ (Bradley et al., 2016). The specific steps of PCR are as below: one cycle of denaturation at 95°C for 4 min, 32 cycles of amplification (95°C for 1 min, 65°C for 2 min, at 72°C for 2 min), and one cycle of extension at 72°C for 10 min. The products were confirmed by 1% agarose gel electrophoresis and purified by the Qiagen Gel Extraction Kit (Qiagen, Germany). The high-throughput sequencing of 18S rRNA gene was conducted on the PacBio Sequel platform by the BioMarker Technologies Company (Beijing, China).

Data analysis

Raw sequences were firstly trimmed to remove the adapters and primers using Cutadapt v2.7 (Martin, 2011) and filtered by LImA v1.7.0 (https://lima.how/) and UCHIME v4.2 (Edgar, 2016) to eliminate most of the low-quality data. All low-quality sequences or unassembled sequences, including chimeric sequences, fuzzy sequences, short sequences (<1200 bp), long sequences (> 2000 bp), and non-specific amplifications, were eliminated for further analysis (Chen et al., 2022). The sequences were clustered into Operational Taxonomic Units (OTUs) based on 97% similarity by using the USEARCH (version 10.0) (Edgar, 2010; Edgar, 2013). Species annotation was mainly conducted in two steps: The feature sequences of OTUs gained by Blast in the QIIME2 software package (Bolyen et al., 2019) were aligned to the SILVA_138 reference database, and then supplementary annotations were obtained with Feature-classifier Classify-sklearn at the confidence threshold of 70% in the QIIME2 software package.

The richness and abundance of species were calculated to reflect the composition characterization of zooplankton communities. Spearman correlations between zooplankton at the phylum or class levels and environmental factors were analyzed by the R package Corrplot (Wei et al., 2021). Correlations between zooplankton OTUs and environmental factors were analyzed and visualized by the R packages Psych (Revelle, 2017) and ggplot2 (Wickham, 2016) respectively. Correlations between zooplankton from phyla to OTU levels were analyzed by R Package Psych and visualized using R Package Circlize (Gu et al., 2014) and ggplot2.

Results

Annual composition of the zooplankton community

A total of 369,875 circular consensus sequences (CCSs) were obtained through PacBio Sequel sequencing, of which 344,712 were retained after quality filtering. After cluster analysis by Usearch, 179 OTUs were obtained by filtering with a threshold of 0.005%. Twelve OTUs were annotated at the kingdom level, while others were annotated to the phylum level or below the phylum levels: Arthropoda (35 OTUs, 92.21%), Chordata (5 OTUs, 3.34%), Mollusca (16 OTUs, 2.22%), Cnidaria (8 OTUs, 0.77%), Ciliophora (6 OTUs, 0.72%), Annelida (14 OTUs, 0.31%), Bryozoa (2 OTUs, 0.12%), Chaetognatha (1 OTU, 0.093%), Echinodermata (1 OTU, 0.079%), Cercozoa (1 OTU, 0.067%), Endomyxa (1 OTU, 0.036%), Platyhelminthes (4 OTUs, 0.032%), Phoronida (1 OTU, 0.007%), Ctenophora (1 OTU, 0.003%), and Porifera (1 OTU, 0.002%) (Figure 1B). The top five zooplankton classes included Copepoda (29 OTUs, 74.76%), Thecostraca (3 OTUs, 17.41%), Appendicularia (1 OTU, 3.14%), Bivalvia (9 OTUs, 1.92%), and Hydrozoa (5 OTUs, 0.74%) (Figure 1C).

A total of 93 OTUs were annotated at the genera and species level, including 86 genera and 89 species. Most genera contain only one species, except for Centropages (3 species), Oithona (2 species) and Ciona (2 species). The Paracalanus parvus, Chthamalus stellatus, and Centropages abdominalis were the top 3 species with extremely high abundance (Figure 1D).

Variations of zooplankton abundance and richness

The species richness (species-level OTU numbers) and relative abundance (total reads of species-level OTUs) analysis of zooplankton samples from May 2020 to May 2021 revealed obvious temporal dynamics (Figure 2). The ANOSIM and NMDS analysis results indicated that samples from the same season were clustered together, and samples in autumn, winter, and spring had overlap, summer samples were independent of other seasons (Figure 2A). The highest richness value appeared in September, while the relative abundance showed the highest in March. Species richness in summer and autumn is much higher than that in spring and winter. Many species exhibit strong temporal dynamics with a seasonal preference, while others are evenly distributed throughout each month. For example, the arthropod C. stellatus mainly presented in May, while Paracalanus parvus was detected in each month with little difference in proportion (Figure 2B).

Figure 2 Temporal dynamics for zooplankton community. (A) NMDS and ADONIS analysis for the zooplankton community, demonstrated the differences in zooplankton composition in different seasons. (B) The number of reads in each month is used as a whole to show the proportion of each species in each month, the legend only shows the species with the largest proportion in each month. (C, D) Relative abundance and richness of top 5 phyla, two charts use the same legend. (E, F) Relative abundance and richness of top 10 class, two charts use the same legend.

With copepods being the dominant species, the species richness of Arthropoda is consistently the highest throughout the year. However, molluscan richness was mainly dominated by Bivalvia and Gastropoda, and showed significant seasonal patterns (high abundance in summer) (Figure 2C). Changes in the abundance of zooplankton communities also showed temporal variations. Among all detected OTUs, the relative abundance of Arthropoda was the highest throughout the year (>83%) (Figure 2D), with Copepoda being the most dominant Arthropoda taxon.

Impact of the environmental parameters and taxa interactions

The annual changes in the seawater temperatures showed an opposite pattern with the DO (Figure 3A), while the salinity (29~31.2) and pH (7.9~8.4) were relatively stable (Figure 3B). The highest temperature was recorded in August, with an average value of 28.2˚C. An abnormally high DO value was recorded during the August sampling, while a low salinity was also recorded. It was suspected to be caused by river flow or windy weather. The canonical correlation analysis (Kordas et al.) indicated that temperature, DO, and salinity significantly influenced the composition of zooplankton communities, in which, the temperature exerts a dominant influence (Figure 3C). Further correlation analysis with the spearman method on the class-level taxa and environmental factors indicated that the abundance of 48 taxa had significant correlations with environmental factors (50 negatives and 43 positives) (Figure 3D). The abundances of different taxa showed different relationships with the environmental parameters. For example, the Acantharea showed an extremely significant (p<0.01) correlation between salinity and temperature, while the Gastropoda abundance was significantly affected by dissolved oxygen and temperature (p<0.01). Gymnolaemata also showed an extremely significant correlation with salinity (p<0.01) besides a significant correlation with temperature (p<0.05) (Figure 3E).

Figure 3 Annual changes in environmental factors and correlation analysis with zooplankton. (A, B) Variation of the environmental factors (Tem, DO, pH and Sal) in the 13 months, take the numerical average of the three sample points. (C) Based on the scatterplot analyzed by CCA, A-M represents the sample collection of different months, showing the similarity between samples in each month and the influence of environmental factors. (D) The Edge Bundling plot based on spearman’s method shows the correlation between each environmental factor and phylum, which only shows a significant correlation (p<0.05), |r| indicates the strength of the correlation, represented by the thickness of the line segment, and Degree represents the number of environmental factors affecting OTU, represented by the size of the point. (E) Based on the sector chart of Spearman’s method, whether the correlation between different classes and environmental factors is significant, red represents negative correlation, blue represents positive correlation, ** represents p<0.01, * represents p<0.05.

Correlation analysis among 12 phyla indicated significant taxa interactions. The correlation between Mollusca and Annelida was extremely significant (p<0.01), while Ctenophora, Endomyxa and Phoronida had a significant positive correlation (p<0.05) (Figure 4A). Significant correlations were also identified among 14 classes. Except for the negative correlation between Spirotrichea and Copepoda, the other pairs were positively correlated (Figure 4B). Similarly, among the 898 OTU-level correlations, 868 were positively correlated, while only 30 were negatively correlated.

Figure 4 Interaction networks between zooplankton OTUs. (A, B) The Spearman method was used to analyze the correlation between zooplankton inside (phylum and class). (C) Pie chart based on the reads number showing the quantitative composition of each species in the dominant class. (D, E) Analysis of inter-species correlation based on OTU data shows the relationship between C. gigas and OTU, the most dominant species and OTU.

Dominant species in Longwan Bay and seasonal variations

Copepods were the year-round dominant species in this study except for May, when the barnacle C. stellatus boomed (Figure 2B). A total of 21 Crustacea species from 18 genera were identified, with the most abundant species from copepods Paracalanus and Centropages (Figure 4C). The maximum abundance of Centropages typicus was at the end of July, with the maximum richness of Centropages abdominalis was observed in January and April. Paracalanus parvus was the most abundant species in other months (Figure 2B). Spearman correlation analysis indicated that the abundances of Labidocera euchaeta, C. typicus, Canuella perplexa, and Acartia pacifica were significantly correlated with seawater temperature (p<0.01). The dominant genera Paracalanus, Centropages, Calanus showed significant correlations with 29 OTUs (Figure 4D). The abundance of C. gigas peaked in July, and showed a significantly positive correlation with 16 OTUs (Figure 4E)

Discussion

The zooplankton community diversity

A total of 89 zooplankton species were identified with the highly reliable 18S rRNA gene long reads sequencing, including the dominant species P. parvus (OTU1 + 172), Oithona similis (OTU7), C. abdominalis (OTU3), Ditrichocorycaeus anglicus (OTU20), and C. typicus (OTU124 + 136). The distribution of these species was consistent with their temperature preference (Hsieh and Chiu, 2004). The species number identified in this study was higher than some other similar studies on aquaculture areas (Guo et al., 2015), which may be a result of different local hydrological conditions (Domínguez et al., 2017) or species identification techniques (Gao et al., 2020). Compared to the traditional morphological identification, the metabarcoding technique has a great advantage in the case of zooplankton larvae or rare species, which leads to the discovery of a larger number of species and facilitates subsequent differential studies (Ayala et al., 2016; Yin et al., 2022). For example, Centropages hamatus and Phyllodoce groenlandica were rarely reported in the studied sea area by traditional morphological identification. However, Acartia bifilosa and Calanus sinicus, the dominant species in the Yellow Sea (Huo et al., 2012; Shi et al., 2020), were not identified in this study. This should be the consequence of the incomplete database which did not deposit the reference sequences of corresponding species (Djurhuus et al., 2018). With the improvement of the molecular marker database, the sequencing result can be re-annotated, and more comprehensive survey results are possibly revealed.

Copepoda is widespread (Cornils and Held, 2014) and one of the most abundant marine organisms in many marine ecosystems (Escribano and Hidalgo, 2000; Hwang et al., 2006). Copepods are thus often used as indicator species for water monitoring (Hsieh et al., 2005; Tseng et al., 2008; Shin et al., 2022). While early studies of copepods are mainly on large copepods (Shimode et al., 2012; Gong et al., 2013), the roles of small and medium-sized zooplankton in marine ecosystems were gradually realized (Middelbo et al., 2019; Tarrant, 2020). Many studies have been conducted on tiny zooplankton species composition, community characteristics, and life habits (Tiselius et al., 2013; Van and Park, 2016; Uttieri et al., 2021). Paracalanus (4 OTUs identified in this study) and Centropages (2 OTUs) are the most common small and medium-sized copepods and are the main contributors to the richness of copepods. They are dominant species in many ocean areas, such as Helgoland, Germany and Maryland Coastal Bays, USA (Oghenekaro et al., 2018; Hirche et al., 2019; Prusova and Galagovets, 2022), as well as the Bohai, Yellow Sea, and East China Sea (Wang et al., 2001; Wang et al., 2002). Previous studies have shown that copepods’ abundance change dramatically impacts the density and diversity of marine zooplankton (Michel and Herring, 1984; Zakaria et al., 2016). Centropages are omnivorous animals, usually preying on eggs and other copepods (Turner et al., 1985; Slater and Hopcroft, 2005). They have been widely identified in multiple sea areas (Dagg and Grill, 1980; Liang et al., 1994). The aquaculture area provides extraordinary biomass and may provide appropriate conditions for the growth and development of Centropages.

Temporal dynamics of zooplankton community

Tiny marine organisms are highly susceptible to environmental factors and show significant seasonal variation (Zhang et al., 2019). The abundance of P. parvus exhibited seasonal variations, consistent with the observations from other areas of the China Seas (Ke et al., 2001; Rong et al., 2002; Zheng et al., 2022). However, the maximum abundance season of P. parvus varied among different reports. The optimum temperature range for the growth of P. parvus is 13-24°C (Zhang et al., 2006), which should be responsible for the different peak occurrences for different sea areas. The abundance of C. typicus in the Longwan Bay reached the highest in summer, differing from the April boom in the southern Mediterranean Sea (Halsband-Lenk et al., 2001). It should also be caused by the local temperature of the corresponding seasons. At the same time, the growth of C. abdominalis starts at 5°C, and the optimal temperature is between 10-15°C (Slater and Hopcroft, 2005). It is also consistent with the seasonal distribution of its abundance observed in this study.

Barnacles (Thecostraca) are one of the most common marine biofouling crustaceans, with planktonic larvae and sessile adults (Pérez-Losada et al., 2012). The abundance of Thecostraca was mainly contributed by C. stellatus in this study. The interaction between oysters and barnacles may impact their abundance and distribution. They are all habited in the intertidal zone, forming partial competition on the limited substrates, which will then lead to a decrease in the population. However, barnacles and oysters showed significant niche differences in this area. The barnacle larvae occurrence (May) was earlier than oyster larvae (July).

Furthermore, barnacles adapt to longer dry exposure environments and mainly distribute around high tide lines. The high abundance of C. stellatus larvae in Longwan Bay indicated that oyster farming did not adversely affect its living. The relative abundance of other Mollusa peaks in March, May, and July, mostly were Mytilus edulis, Musculus cupreus, and Hiatella arctica. The molluscan abundance peaks in March and May were mainly contributed by the larvae of mussel M. edulis, which is consistent with its common spawning season (Nordsieck, 2006; Tyler-Walters, 2008). At the same time, the July molluscan occurrence peak was mainly contributed by M. cupreus larvae. The study results indicated a relatively low proportion of oyster larvae in the collected zooplankton samples, even within the oyster farm region. Oyster culture in Longwan Bay should have a limited impact on the marine zooplankton community. It should also benefit from the well-controlling aquaculture density in Longwan Bay.

Effects of environmental factors and species interactions

The diversity of species distribution in marine ecosystems is affected by physicochemical and biological factors (Tittensor et al., 2010), where temperature is the main controlling factor (Gillooly, 2000; Lewandowska et al., 2014; Yasuhara and Danovaro, 2016). Temperature determines the general dissolved oxygen level and affects the composition and species interactions of the zooplankton community (Portner and Farrell, 2008; Kordas et al., 2011). Temperature is also a fundamental determinant of organisms’ life cycle. Too high or too low temperatures can be associated with a high mortality rate (Heinle, 1969). Studies have suggested that increasing temperature accelerates zooplankton metabolism and sexual maturation and reduces zooplankton lifespan (Huang, 1985).

Adult survival and larval occurrence were also influenced by changing ocean temperatures. The Pacific oyster C. gigas is one of the main mariculture species in Longwan Bay, whose optimal water temperature for growth is 15-25°C. According to our field investigation, this bay is very close (~100 km) to the southern distribution limit (around Lianyungang, Jiangsu Province) of C. gigas. Severe summer mortality occurred in recent years, possibly because of environmental changes caused by global warming. The high temperature observed in August (~28°C) has exceeded the optimal temperature range of the Pacific oyster and could cause heat stress. Correspondingly, the development and maturity of oyster gonads may also be influenced by environmental change. The swimming larvae of C. gigas were mainly sampled at the end of July, indicating a dominant mass spawning during July. The newborn oyster spats will then face heat stress during August. If global warming persists, the wild Pacific oyster population and distribution limits may change in this area.

The composition of tiny zooplankton communities is also driven by their biological activities, which are sensitive to temperature and dissolved oxygen (Bērziņš and Pejler, 1989; Zhao et al., 2020). In marine ecosystems, many animals (e.g. ctenophores) feed on zooplankton and have a wide range of adaptations to DO. However, the movement frequency of copepods is reduced in the relatively low DO environment, which improves the predation rate and thus affects the composition of the community (Mary Beth et al., 2004). A previous study indicated that both the survival rate and the reproduction of the copepod Acartia tona and Oithona colcarva were inhibited at a low DO state (Roman et al., 1993).

Indirect and direct species interactions are widely present in marine communities and play an important role in the formation and maintenance of biodiversity (Sun, 2011). Zooplankton is crucial for the marine ecosystem as a hub in species interaction networks (Greve, 1977), where Spirotrichea is significantly associated with Copepoda. Spirotrichea is a microplankton which feeds on photosynthetic algae and bacteria and then as the primary food for copepods and young fish (Calbet and Saiz, 2005; Lima-Mendez et al., 2015; Santoferrara et al., 2017). OTUs in Ascetosporea, the parasitic protozoans infecting aquatic invertebrates (Bass et al., 2019), were significantly correlated with OTUs in Tentaculata. The species interactions revealed by this study indicated that food chain and parasite-host relationships also play important roles in stabilizing plankton communities.

Conclusion

The special and complex environment of the ocean is an important driving force for the formation and evolution of marine biodiversity. Human activities have created new challenges for marine zooplanktons, which were driven to constantly evolve many adaptive traits: morphological, behavioural, physiological traits and reproductive patterns (Shen and Shi, 2002). Species living in the same environmental pressure also exert selective pressure on each other to form a relatively stable biome. The metabarcoding analyses of the zooplankton community in this study not only lay the foundation for the spatiotemporal dynamic mechanism of zooplankton in the aquaculture area but also provide supporting evidence for further understanding of the impact of human activities and global climate change on biome composition. This study provides highly reliable 18s full-length sequence data and reveals major planktonic organisms’ structure and seasonal dynamics in Longwan Bay. The results will be helpful for further monitoring the planktonic community change and improving biofouling management efficiency.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: BioProject, PRJNA947076.

Author contributions

FX conceived and designed the study. XM and BQ collected and processed samples. BQ conducted data analysis. BQ, XM and FX wrote the paper. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the Key Research and Development Program of Shandong (2022LZGC015), the Key Deployment Project of Centre for Ocean Mega-Research of Science, Chinese Academy of Science (No. COMS2019Q11), and the Taishan Scholars Program.

Acknowledgments

We acknowledge Chonglin Li, Liangheng Li and Shaoxi Deng for their assistance with the sample collection. Most of the calculations were supported by the Oceanographic Data Center, IOCAS.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References