Coastal Eukaryotic Plankton Diversity of the Southern Adriatic as Revealed by Metabarcoding

Abstract

Plankton studies serve as a basis for marine ecosystem research, but knowledge of marine plankton is still incomplete due to its extreme taxonomic and functional complexity. The application of metabarcoding is very valuable for the characterisation of the plankton community. The plankton community of the Southern Adriatic is subject to strong environmental fluctuations and changes, which underlines the need for frequent, reliable and comprehensive characterisation of the plankton. The aim of this study was to determine the taxonomic composition and seasonal distribution of eukaryotic plankton in the Southern Adriatic. Plankton samples were collected monthly for one year at the coastal station of the Southern Adriatic and metabarcoding was used for taxonomic identification. The results showed a high taxonomic diversity and dynamic seasonal distribution patterns for both the protist and metazoan plankton communities. Metabarcoding revealed both the core, year-round plankton community and previously unrecorded plankton organisms in the Southern Adriatic. The results provide for the first time a comprehensive overview of the plankton community in this area by metabarcoding. The identified seasonal patterns of plankton genera and species in the Southern Adriatic will contribute to the understanding of plankton interactions and future changes in community diversity characterisation.

1. Introduction

Marine eukaryotic plankton comprises a huge diversity of organisms that have adapted to drift and are transported over long distances by ocean currents [1]. The diversity of marine plankton is reflected in taxonomic composition, body size and life history (meroplankton vs. holoplankton). Microscopic unicellular eukaryotes (protists or microeukaryotes) form an important part of the marine plankton community as primary producers (phytoplankton), consumers, parasites and symbionts [2,3]. Protists are often difficult to identify morphologically, as they are small and sometimes do not differ in their morphological characteristics from related taxa. Multicellular heterotrophic organisms (metazoans) participate in the plankton community primarily as consumers (zooplankton) and are often restricted to certain life stages (meroplankton). Although metazoans have a significantly larger body size than protists, many of them are still too small for microscopic identification, and there are also morphological similarities between related taxa (e.g., copepods and invertebrate larval stages). The planktonic larval stage plays a crucial role in the dispersal of many benthic invertebrates, but very little is known about the larval biology of most species. Often the taxonomy of larvae and adults is not yet unified, further complicating ecological analyses. Plankton functions as a complex ensemble of organisms that constantly interacts with environmental factors and organisms of other marine habitats (benthos, nekton) and, therefore, plays a fundamental role in marine food webs and the chemical and physical cycles of the oceans [4]. The interactions within the plankton community also determine the diversity and dynamics of the community. Because of this fundamental role, accurate information on the diversity and changes in the structure of the plankton community in space and time is needed. Plankton biodiversity is crucial for the functioning, sustainability and resilience of marine ecosystems. However, for most marine areas, the spatial and temporal patterns and taxonomic details of the plankton community are not defined.

Traditionally, plankton is identified using light microscopy. However, due to the small cell size, similar morphology and low abundance of many important components, an additional and complementary methodological approach such as metabarcoding reveals further details. Metabarcoding in combination with high-throughput sequencing (HTS) enables the simultaneous detection of multiple taxa over a large systematic range by using universal or more specific PCR primers to amplify taxonomically informative gene regions from environmental DNA [5]. Although metabarcoding (like any method) is subject to certain methodological biases (e.g., in sampling, pre-processing, DNA extraction, marker and primer selection, PCR amplification and sequencing), correct data analysis and interpretation that takes these biases into account leads to a highly reliable and comprehensive description of the plankton community. Metabarcoding and HTS have been used and described to study global plankton diversity and to recognise patterns of taxa occurrence. Metabarcoding studies have described spatial and temporal variation in complete eukaryote plankton assemblages [6,7,8], but more often studies have focused on selected plankton taxa [9,10,11,12,13,14,15]. As at the global scale, only those groups and taxa representative of the Adriatic plankton community have been analysed using metabarcoding [16,17,18,19,20,21,22,23,24], neglecting insight into the eukaryotic plankton community as a whole. For the Southern Adriatic, previous plankton research in this area also included mainly zooplankton and phytoplankton studies on certain representative groups. Zooplankton studies focused on copepods [25,26,27,28,29], tintinnids [30,31,32] and gelatinous invertebrates [33,34,35,36,37], while phytoplankton studies mainly characterised diatoms [38,39,40]. Moreover, the studies were predominantly based on traditional methods of taxa identification. Considering all this, there is a considerable lack of knowledge about the Southern Adriatic plankton diversity that has been recorded by metabarcoding. Our study aims to provide the first comprehensive description of the eukaryotic plankton community of this ecologically important marine region.

The Adriatic Sea is the semi-enclosed, northernmost part of the Mediterranean. Due to its geographical and oceanographic characteristics, the Adriatic Sea is an extremely dynamic but also vulnerable ecosystem. Biodiversity and climate changes affecting larger marine areas (Mediterranean and global oceans) are first reflected/recognised in the Adriatic, which underlines the importance of the Adriatic in monitoring the effects of climate change [41]. The Southern Adriatic is the deepest part of the Adriatic Sea and a key area for the regular exchange of water between the Adriatic and the Ionian Sea through the Strait of Otranto. The presence of a quasi-permanent cyclonic circulation [42] is characteristic of the Southern Adriatic, and these circulation events have a strong influence on the Southern Adriatic ecosystem and the plankton community [34]. Recent changes in the plankton composition of the Southern Adriatic have been linked to large-scale climatic changes in the area [34]. The Southern Adriatic is the entry point and the first capture area for potential new arrivals and invasive species in the Adriatic ecosystem [43,44]. As the Southern Adriatic is a dynamic environment, and the Adriatic Sea is the gateway to further marine influences from the Mediterranean and beyond, the characterisation of the SA plankton community is crucial for indicating changes and predicting the evolution of the community for the entire Adriatic Sea. The description of the Southern Adriatic plankton community contributes to the understanding of plankton interactions and changes in community diversity, especially in detecting new, non-indigenous members of the community over the course of future climate change scenarios.

The main objective of this manuscript is to describe the diversity of the Southern Adriatic eukaryotic plankton community in its entirety, resulting in a detailed, reliable and comprehensive insight into the community structure and pointing to the current community interactions. The applied monthly sampling strategy, which enables the description of the plankton temporal pattern contributes to the general knowledge of plankton biodiversity and the further development of regional monitoring plans.

2. Materials and Methods

2.1. Sampling and DNA Isolation

Samples were collected at the Southern Adriatic coastal station (LOKRUM, depth 100 m, GPS coordinates: 42°37′21.0″ N, 18°06′00.0″ E) located in the vicinity of Lokrum Island, near Dubrovnik (Figure 1A). Vertical hauls with the plankton net (53 µm) equipped with the closing mechanism were performed at a monthly frequency for a one-year period (March 2021–February 2022). Plankton was sampled in two layers, namely surface (0–50 m) and subsurface (50–100 m). The sampling of both layers encompasses the euphotic zone (up to 100 m). Net samples were filtered on 1.2 µm cellulose filters and stored in an RNAlater at −20 °C until genomic DNA isolation. Isolations of the genomic DNA were performed using the NucleoSpin eDNA Water kit (Macherey-Nagel, Germany) according to manufacturer instructions. The quality and quantity of the isolated DNA were checked spectrophotometrically using nanodrop (Implen, Germany), and the samples were sent for metabarcoding library preparation and sequencing. DNA metabarcoding library preparation and sequencing were carried out by AllGenetics and Biology SL, Spain (www.allgenetics.eu, accessed on 4 May 2021).

2.2. DNA Metabarcoding Library Preparation and Sequencing

For library preparation, primers that amplify a region of around 430 bp of the V4 region of the 18S gene were used: TAReuk454FWD1 (5′ CCA GCA SCY GCG GTA ATT CC 3′) [45] and TAReukREV3 (5′ ACT TTC GTT CTT GAT YRA TGA 3′) [45,46]. These primers have the Illumina sequencing primer sequences attached at their 5′ ends. In the first amplification step, PCRs were carried out in a final volume of 12.5 μL, containing 1.25 μL of template DNA, 0.5 μM of the primers, 3.13 μL of Supreme NZYTaq 2x Green Master Mix (NZYTech, Portugal) and ultrapure water up to 12.5 μL. The reaction mixture was incubated as follows: an initial denaturation step at 95 °C for 5 min, followed by 35 cycles of 95 °C for 30 s, 48 °C for 45 s, 72 °C for 45 s and a final extension step at 72 °C for 7 min. The oligonucleotide indices that are required for multiplexing different libraries in the same sequencing pool were attached in a second amplification step with identical conditions but only for 5 cycles and at 60 °C annealing temperature [47]. A negative control that contained no DNA (BPCR) was included in every PCR round to check for contamination during library preparation. The library size was verified by running the libraries on 2% agarose gels stained with GreenSafe (NZYTech) and imaging them under UV light. Then, the libraries were purified using the Mag-Bind RXNPure Plus magnetic beads (Omega Bio–tek, Georgia, USA), following the instructions provided by the manufacturer. Finished libraries were pooled in equimolar amounts according to the results of a Qubit dsDNA HS Assay (Thermo Fisher Scientific Inc., Germany) quantification. These pools also contained a testimonial amount (1 μL) of the PCR negative controls. The pool was sequenced in a fraction of NovaSeq PE250 (Illumina, USA) aiming for a total output of 1 gigabase. Illumina paired-end raw data consist of forward (R1) and reverse (R2) reads stored in separate files, which are sorted by library and include the reads’ quality scores. We assessed the quality of the FASTQ files with the software FastQC, version 0.11.9 [48] and summarised the output using MultiQC, version 1.12 [49]. The raw FASTQ files were used for further bioinformatics steps.

2.3. Bioinformatics

Processing of raw reads for quality trimming and Amplicon Sequence Variants (ASVs) identification was conducted using R (R version 4.3.1, www.R-project.org, accessed on 4 July 2023) and the dada2 package (DADA version 1.26, [50]). After inspecting the read quality profiles (plotQualityProfile) of forward (R1) and reverse reads (R2), the quality-filtering step is completed with the filterAndTrim() function using the parameters maxEE = 2, truncQ = 5, maxN = 0 and trimLeft = 20. With the parameter truncLen, the reads were truncated at 240 nt (R1) and 210 nt (R2), maintaining overlap after truncation in order to merge the forward and reverse reads later. All the other arguments in filterAndTrim() were set to the default values. Additionally, within the dada2 quality-filtering step, the argument trimLeft was used to remove primers. The quality of the filtered files was also checked, and it was confirmed that truncation excluded average qualities Q-scores read areas <30. The error rates were learned using the errF and errR functions. The minimum number of total bases to use for error-rate learning was set to the default value (nbases = 1 × 10⁸). The maximum number of times to step through the self-consistency loop was also set to the default value (MAX_CONSIST = 10). The sample inference algorithm was applied to the filtered and trimmed sequence data, and ASVs were inferred using the Dada function. Paired reads were merged with a minimum overlap set to 15 nucleotides in the mergePairs function. Merged sequences outside the expected range (300–450 bp) were discarded. Chimeric sequences were removed using the removeBimeraDenovo function and consensus method. Taxonomy assignment is implemented using the naive Bayesian classifier method [51] and the PR2 database (PR2 version 4.14.0), a curated list containing only eukaryotic taxa [52]. The minimum bootstrap confidence for assigning a taxonomic level has been set to 90. A BLAST search [53] of unclassified ASVs (the taxonomy assigned to the level higher than the genus by DADA and the PR2 database) against the nonredundant nucleotide GenBank reference database [54] was conducted.

2.4. Community Diversity and Composition Analyses

Plankton community composition analyses, statistical analyses and generation of figures were conducted using R (R version 4.3.1, www.R-project.org) and the packages vegan [55], ggplot2 [56], phyloseq [57] and netcomi [58]. The plankton ASV dataset was reduced by removing ASVs classified as land-originated taxa. ASVs classified to taxonomic ranks higher than genus were also excluded from the dataset. Samples with very low read abundances (<300 reads; samples JUG3 and JUG22) were removed from the dataset. The number of reads in each sample was normalised using median sequencing depth. For taxa abundance visualisations (bar plots), sample counts of a taxa abundance matrix were transformed to relative abundances. The counts of each sample were transformed individually. Each sample plankton community’s alpha diversity was calculated as ASVs richness and diversity indices (Shannon and InvSimpson index) value. An analysis of variance (ANOVA) was used to test for differences in alpha diversity parameters among different seasons and sampling depths. Patterns of sample dissimilarity (beta diversity) (NMDS plot) were analysed based on ASVs and genera presence–absence (Jaccard distance), visualised using unconstrained ordinations of non-metric multidimensional scaling (NMDS), and statistical confidence was tested with multivariate PERMANOVA (permutations = 9999). Permutational multivariate analysis of variance (PERMANOVA) is a common method for assessing whether a variable is associated with an overall difference in sample community composition. ASVs were merged at the genus level (tax_glom function in phyloseq), and subsets of dataset groups (protist and metazoan) were used for heatmap and network composition analyses. For seasonal taxonomic composition analyses (hierarchical clustering and heatmap analyses), frequencies of occurrence for genera [59] were calculated as the number of times a genus occurred in a time frame (season) divided by the total number of samples in that time frame (season). To generate heatmaps, we applied agglomerative clustering specifically with the clustering method ward.D2, and correlation was used as a clustering distance measure with the R package pheatmap. The definition of seasons in the study followed the meteorological (climatological) seasons in the Northern Hemisphere (https://www.ncei.noaa.gov/news/meteorological-versus-astronomical-seasons, accessed on 4 July 2023) for spring (March, April, May), summer (June, July, August), autumn (September, October, November) and winter (December, January, February).

The plankton community association network (nodes are taxa/genera) was constructed from the phyloseq object using the compositionally aware correlation estimator SparCC [60] and the R package netcomi [58]. Read counts data were log-ratio transformed and the Bayesian approach was applied for zero handling. A threshold of 0.4 was used as a sparsification method so that only genera with an absolute correlation greater than or equal to 0.4 are connected. To detect clusters in the network, cluster optimal modularity optimisation was used, and unsigned transformation was set so that the edge weight between strongly correlated genera is high, no matter the sign (+ or − correlation). Betweenness and degree centrality measures were used for identifying hub nodes. Betweenness centrality measures the importance of a node in connecting other nodes in the network, thus capturing the role of a genus as an intermediary between other genera in the network. Degree centrality describes only the number of connections (edges) a genus has established with other genera in the network. Betweenness considers all geodesics between two nodes. R package igraph [61] was used for visualizing and analysing the network. Global network properties were defined for the whole network.

An online tool of The Interactive Tree Of Life (https://itol.embl.de, accessed on 4 July 2023) [62] was used to display and annotate common phylogenetic trees of the marine eukaryotic plankton community.

2.5. Data Availability

The data for this study have been deposited in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB75449.

3. Results and Discussion

3.1. Overview of the Diversity

High-throughput sequencing of the Southern Adriatic dataset (24 samples) resulted in 1.945,388 raw reads overall. After bioinformatics processing, 1,267,410 reads, corresponding to 932 ASVs, were used for the taxonomy assignment. A large part of the ASVs were assigned to Metazoa (711 ASVs and 1.251,960 reads), followed by protists (microeukaryotes) (182 ASVs and 13,789 reads) (Figure 1B). The remaining 39 ASVs assigned to land-originated taxa of Embryophyta (5 ASVs, 161 reads), Fungi (10 ASVs, 987 reads), seaweed Rhodophyta (2 ASVs, 25 reads) and unclassified Eukaryota (22 ASVs, 481 reads) (Figure 1B) were removed from further plankton community analyses. In the dataset, 475 ASVs (48% of all ASVs in the dataset; 62 protists and 413 metazoan ASVs), represented by 649,177 reads (51% of all dataset reads), were identified to a taxonomic rank higher than genus. Since our focus was to characterise plankton genus and/or species assemblages, a subset of ASVs reliably identified at the species or the genus level was produced. This finalised (genus/species level) Southern Adriatic plankton dataset contained 418 ASVs, with an average of 54 ASVs per sample. Rarefaction curves of all samples reached saturation, indicating that the sequencing depth was sufficient (Figure 1C). Unclassified ASVs (taxonomy assigned to levels higher than genus) were mainly assigned to unclassified Metazoa (copepods predominantly). A BLAST search found that 32% of these unclassified ASVs show high sequence identity (≥99%) to marine uncultured and environmental clones from different metagenomic surveys. This part of the dataset represents a reservoir of hidden diversity that is yet to be uncovered, as species isolation and barcoding will progress and enrich the reference databases used in metabarcoding. Future research hence will be able to mine the here-published dataset for the early presence of species for which, currently, no information or no sequence information is available.

In the finalised (genus/species level) Southern Adriatic plankton dataset, the highest ASV number was detected in winter, and the lowest was in summer samples (Figure 1D). Between seasons, significant differences were recorded for the detected number of ASV (ANOVA, p = 0.026 *, Tukey test = winter-summer, p adj = 0.0163399). The Shannon Index ranged from 3.2 in winter to 1.0 in the autumn sample and InvSimpson from 1.7 in spring to 13.6 in the winter sample (Figure 1D), but no significant differences (ANOVA, p < 0.05) in both diversity indices were present between the seasons. Sampling depth did not have a significant effect on any of the alpha diversity measures during the year. The average numbers were 56 and 5.9 for the surface (50–0 m) and 53 and 7.2 for the subsurface (100–50 m) layer for the ASV number and InvSimpson Index, respectively. The average number of the Shannon Index was 2.3 for both sampling layers (Figure 1D).

Community composition based on ASVs (Figure 2A) or genera (Figure 2B) differed significantly between seasons (PERMANOVA, p = 1 × 10⁻⁴), and 23% (ASVs) or 25% (genera) of the variation was explained by the season parameter. Separated clusters of plankton communities were revealed with NMDS analysis for the four seasons, with different seasonal similarities present for the ASV and genera-based NMDS (Figure 2). When the distance between samples was based on ASV presence/absence, the summer samples shared higher similarity with the spring and autumn samples, while the winter samples were grouped separately from other seasons. On the contrary, when the distance between samples was based on genera presence/absence, the similarity between summer and autumn, as well as the winter and spring samples, was indicated. A significant effect of the sampling depth on community composition was also present (PERMANOVA, p = 0.0191) for both ASVs and genera-based analyses but contributed only 7% of the variation in community distances (Figure 2).

3.2. Taxonomic Structure

The plankton community of the Southern Adriatic Sea consisted of 19 major taxonomic groups (pr2 database divisions), including metazoans (11 different groups) and 8 protist groups (Figure 3). A total of 298 ASV (93 different taxa; all at species and genus level) were assigned to metazoans, while 120 ASV (58 different taxa; all at species and genus level) were assigned to protist groups. Metazoan ASVs were uniformly abundant throughout the year. Among the protists, dinoflagellates dominated in abundance and were detected in all samples throughout the year and at all depths. Similarly, ciliates (Ciliophora) were detected in all seasons, although in lower abundance than dinoflagellates. Most of the other eukaryote groups showed distinct distributions across the seasons (Figure 3). The protist group Bacillariophyta was detected evenly in both sampling layers, but mainly in the spring and summer months and only in low abundance in winter. The protists of the Chlorophyta group were detected in all seasons except autumn, but in low abundances (<100 read/sample). A sister group of Chlorophyta, the Prasinodermophyta group, was only detected in the surface layer of the summer samples. Zooplankton protists of the Cercozoa and Radiolaria groups were present in all seasons except summer, and although they were generally recorded in low abundances (<100 read/sample) in the dataset, these groups contributed significantly to diversity in autumn. Among the other heterotrophic protists, a group of pico-sized ubiquitous marine stramenopiles (MASTs) was recorded throughout the water column in winter in the Southern Adriatic dataset (Figure 3).

About 20% of all recorded genera in the Southern Adriatic metabarcoding dataset were shared by all seasons; they were present throughout the year. These 27 shared genera (copepods: Acartia, Calocalanus, Candacia, Centropages, Clausocalanus, Corycaeus, Euterpina, Haloptilus, Isias, Lucicutia, Mecynocera, Oithona, Paracalanus, Pareucalanus, Sapphirina and Temora; shrimps: Meganyctiphanes; hydrozoans: Clytia, Liriope, Solmundella and Sulculeolaria; dinoflagellates: Blastodinium and Dino-Group-I-Clade-1; annelids: Polyophthalmus and Tomopteris; gastropod: Creseis; and tunicate: Doliolum) represent a year-round community in the southern coastal zone of the Adriatic Sea, which was detected by metabarcoding. As expected, most organisms of the year-round community spend their entire life cycle in the plankton (holoplankton); also, organisms of the meroplankton such as the annelids of the genus Polyophtalmus (Polyophthalmus pictus) formed the year-round plankton community of the coastal Southern Adriatic, suggesting that the eggs and/or larval stages of this genus/species in the plankton are not restricted to a specific season.

Together with the year-round community, autumn in the Southern Adriatic was characterised by the highest number of genera occurring only in a given season (24 genera, corresponding to 18% of all identified genera), followed by winter (16 genera, 12%), while spring had the lowest number of unique genera (6 genera, 4.5%). Between the seasons, winter and autumn shared the highest number of genera (12 genera, 9%), while spring and summer, as well as summer and winter, shared the lowest number of genera, only 3 genera (2%).

Among the most common marine eukaryotic plankton organisms [63] all major protist (Alveolata, Chlorophyta, Ochrophyta and Radiolaria) and metazoan (Annelida, Arthropoda, Bryozoa, Chordata, Cnidaria, Ctenophora, Echinodermata, Mollusca and Phoronida) groups (divisions) had genus/species representatives in the Southern Adriatic dataset (Figure 4), confirming a broad taxonomic range covered by 18S metabarcoding. In addition, the metabarcoding dataset enabled the detection of some marine eukaryotic plankton groups that are usually less frequently studied and/or taxonomically characterised, such as the metazoan phyla Nemertea and Platyhelminthes.

3.2.1. Protist (Microeukaryotes) Community Structure

The protist community, which was assigned to eight taxonomic groups, numbered 58 taxa (39 genera and 19 taxa from global environmental datasets) (Figure 5). For the coastal Southern Adriatic, a seasonal clustering of the protist community into two main groups (spring–summer and autumn–winter) was observed. The highest number of protist genera was shared between autumn and winter (Figure 5). Furthermore, autumn was characterised by the highest number of genera detected only in this season (autumn unique). Only two genera were found in all four seasons. The seasonal frequencies of occurrence of protist genera were generally low (Figure 5). There were more unique than shared genera in all seasons. Nevertheless, most protist groups showed seasonal differences in the composition of assigned genera (Figure 5). However, when characterising the seasonal patterns of protist groups’ characterisation, their low read abundance should be taken into account. In addition, the occurrence of the recorded genera should be confirmed and analysed over a period of several years.

The Chlorophyta were represented in the dataset by two taxa of prasinophytes (genus Chloropicon and prasinophyte clade IX), both of which showed a higher frequency of occurrence in winter (Figure 5) but were also recorded in other seasons except autumn. The picophytoplanktonic prasinophyte genus Chloropicon is the first record for the Adriatic. The metabarcoding record of another picophytoplankton genus, Prasinococcus (Prasinodermophyta), presented here is also the first record for the Adriatic. This genus of a non-flagellate coccoid alga was only detected in the summer (August 2021) during our one-year study. New discoveries of genera and species of coccoid and flagellate picoplankton in the Southern Adriatic indicate a large protist diversity in this area, which is still unexplored due to the morphological similarity of genera/species and limited morphological identification possibilities [44,64].

During the year, Cercozoa (genera Aulacantha and Coelodendrum) and Radiolaria (genera Sphaerozoum, Disolenia, Collozoum and Eucyrtidium) occur more frequently in the colder seasons. The annual occurrence of the radiolarian genera correlates with the occurrence of the dinoflagellate genus Brandtodinium, which is known to be a common symbiotic partner of these polycystine radiolarians in symbiotic interactions between pelagic hosts and microalgae [65]. Brandtodinium was only recently described as a new dinoflagellate genus [65], and this metabarcoding record represents a new occurrence for the Southern Adriatic. As Brandtodinium has only recently also been metabarcoded for the northern Adriatic [24], this may indicate a wider and more common occurrence in the Adriatic that was previously hidden due to microscopic identification limitations. While species of the genus Eucyrtidium have been described in previous studies as characteristic of the euphotic layer in the Southern Adriatic [66], the genera Sphaerozoum (Sphaerozoum trigenimum and Sphaerozoum fuscum), Disolenia (Disolenia quadrata) and Collozoum (Collozoum sp.) represent new records for the study area. All Cercozoa genera recorded in the metabarcoding dataset of the Southern Adriatic belong to the class Phaeodarea, a group of cosmopolitan protists traditionally categorised as radiolarians (Haeckel’s Radiolaria), but molecular phylogeny assigns them to the Cercozoa [67]. Phaeodarea (phaeodarians) are quantitatively the most important radiolarians in the Adriatic, and both Aulacantha scolymantha and Coelodendrum ramosissimum have already been recorded in the Southern Adriatic [66]. The diversity and seasonal occurrence (winter/autumn) of Cercozoa and Radiolaria recorded in this study using metabarcoding were consistent with previous studies of these two protist groups using the classical identification methodology [66].

Figure 5. Heatmap representing Southern Adriatic protist genera seasonal frequency of occurrence. Clustering of the dataset is indicated as dendrograms showing clustering of samples (merged according to season) with similar community structure and clustering of genera with similar seasonal occurrence patterns. Division affiliation of each genus is indicated on the heatmap side edges. (A). Upset plot of the shared protist genera between the seasons (B). There has never been a consensus with regard to whether or not Salpingacantha are morpho-types of Salpingella or distinct species [68], so the two taxa were, in our analysis, left as separate genera.

Figure 5. Heatmap representing Southern Adriatic protist genera seasonal frequency of occurrence. Clustering of the dataset is indicated as dendrograms showing clustering of samples (merged according to season) with similar community structure and clustering of genera with similar seasonal occurrence patterns. Division affiliation of each genus is indicated on the heatmap side edges. (A). Upset plot of the shared protist genera between the seasons (B). There has never been a consensus with regard to whether or not Salpingacantha are morpho-types of Salpingella or distinct species [68], so the two taxa were, in our analysis, left as separate genera.

Dinoflagellates dominated the protist dataset of the Southern Adriatic, with 48% of the total protist abundance, 76% of ASVs and 35% of all protist taxa. Among the Southern Adriatic dinoflagellates, parasitic (Blastodinium, Chytriodinium, Ellobiopsis, Hematodinium, Syndinium and environmental Dino-Group-Clades (Syndiniales)) and free-living (Goniodoma, Heterodinium, Polykrikos, Tripos and Warnowia) genera were recorded. In addition to Brandtodinium (Brandtodinium nutricula), Pelagodinium (Pelagodinium beii) as a foraminifera symbiont [65,69] and Ellobiopsis (Ellobiopsis chattonii) as a common parasite of pelagic marine copepods (Acartia, Calanus, Centropages) [70] were recorded for the first time in the Southern Adriatic. As Pelagodinium has only recently been metabarcoded for the Northern Adriatic [24], its detection in the Southern Adriatic may indicate a broader and more common occurrence in the Adriatic, which was previously hidden due to microscopic identification limitations. Blastodinium, Syndinium and Chytriodinium are cosmopolitan genera of parasitic dinoflagellates, also known as copepod parasites in the gut or lipid-rich copepod eggs [71,72]. As parasites have been shown to influence the mortality and fecundity of copepod populations [73] and copepods are the most abundant metazoans in the ocean and represent an important trophic link in pelagic food webs [74], the successful identification of copepod parasites in the Southern Adriatic dataset contributed significantly to the understanding of the diversity and health of the marine plankton community. The advantage of metabarcoding over classical microscopy in the identification of genera and species of parasitic dinoflagellates also lies in the fact that parasitic dinoflagellates often deviate from the typical dinoflagellate morphology (motile, bi-flagellated protist nucleus and dinokaryon as the nucleus) and change their morphology during their life cycle [71]. Parasitic dinoflagellates (Blastodinium, Hematodinium and Dino environmental clades (Syndiniales)) were the most abundant dinoflagellate taxa (>100 read/sample) in the Southern Adriatic dataset. The Syndiniales environmental sequences belonged to Groups I and II, further assigned to eleven clades. Sequences from global environmental datasets are also most frequently assigned to Groups I and II [75]. The dinoflagellate genus Blastodinium occurred in all seasons, as did the (Syndiniales) Dino-Group- I-Clade-1. Other Syndiniales clades were recorded interchangeably throughout the year (Figure 5A). Overall, dinoflagellate taxa were recorded in all seasons but with a higher taxa richness in autumn (15 taxa) and winter (12 taxa) compared to summer (9 taxa) and spring (8 taxa). The parasitic/symbiotic genera were clustered in autumn–winter, whereas the free-living genera occurred more frequently in a spring–summer cluster. A similar seasonal occurrence of dinoflagellates, with a high abundance and diversity of parasitic taxa in winter, has already been observed for other Mediterranean coastal study sites [13,46].

After dinoflagellates, diatoms (Bacillariophyta) were the second most abundant protist group in the dataset, accounting for 33% of the assigned protist reads. With the exception of Asterionellopsis as a winter-only diatom genus and Chaetoceros, which was detected throughout the year (except in autumn), the other diatom genera recorded (Rhizosolenia, Pseudo-nitzschia, Thalassiosira, Meuniera and Cerataulina) belonged to the spring–summer cluster (Figure 5). The more frequent presence and higher abundance in the warmer seasons (spring–summer) of Bacillariophyta were also found in previous Mediterranean plankton studies using metabarcoding [20,46] and also light microscopy [76,77] as taxa identification methods. The confirmation of Chaetoceros as an all-year present genus is also consistent with previous metabarcoding studies in the Mediterranean [21]. Interestingly, our metabarcoding dataset lacks autotrophs (phytoplankton), except for Bacillariophyta, Chlorophyta and some dinoflagellate taxa. Although these missing phytoplankton groups (autotrophs), such as Chrysophyceae, Dictyochophyceae, Haptophyta and Cryptophyta (Figure 4), usually contribute with lower read abundances (less than 10%) and taxa richnesses, their presence is usually successfully recorded with 18S barcodes in Mediterranean metabarcoding datasets [20,46]. The underrepresentation and non-detectability of autotrophic taxa in the coastal protist community of the Southern Adriatic should be further investigated, and possible methodological biases (sampling, sample preservation, DNA isolation and PCR amplification) should be considered. Although the methodology applied in our study provided a relevant and unified overview of the eukaryotic biodiversity of the Southern Adriatic Sea, due to the taxonomic complexity of the eukaryotic plankton [1], all methodological steps may still be unsatisfactory for investigating the composition of many plankton groups.

Among protists, ciliates were less abundant (7% of protist reads) but had similar ASV (Ciliophora 30, Bacillariophyta 21) and taxa richness (Ciliophora 13, Bacillariophyta 7) as diatoms. Ciliates were found in all seasons, and most taxa (61%) were assigned to the tintinnids (Tintinnida). The tintinnids (Xystonella, Salpingacantha, Salpingella, Codonellopsis, Dictyocysta and Rhabdonella) were characteristic of autumn–winter, but some genera were also found in summer (Amphorides and Eutintinnus). The tintinnids detected in our metabarcoding dataset are characteristic genera of the Adriatic surface zooplankton community (coastal and open waters) [30,31]. The ciliates from our Southern Adriatic dataset also included parasites of marine organisms (order Apostomatia; genera Chromidina and Synophrya) and environmental clades (PHYLL_4).

3.2.2. Metazoan Community Structure

To characterise the structure of the metazoan community, 298 metazoan ASVs (588,845 reads) associated with 11 divisions and counted for 73 different genera were used (Figure 6A). The metazoan taxa found characterise the meroplanktonic and holoplanktonic representatives of the zooplankton community in the Southern Adriatic (Figure 6A). For the coastal Southern Adriatic, a seasonal clustering of the metazoan community into two main groups (spring–winter and summer–autumn) was found, but a dominant clustering between all-year present and seasonally present genera was also evident. As for protists, the highest number of metazoan genera detected only in one season was found in autumn (13 unique genera) (Figure 6B). The phylum Arthropoda dominated the metazoan dataset in terms of abundance (505,433 reads; 86% of metazoan reads) and diversity (209 ASVs and 41 taxa). Most of the arthropod taxa belonged to the marine copepods (39 taxa, orders Calanoida, Cyclopoida and Harpacticoida). Cyclopoid and Harpacticoid copepods showed lower abundances and a lower taxa richness. Harpacticoid copepods were represented with only one genus, Euterpina, in the Southern Adriatic dataset. Among the Arthropoda, only two taxa were identified that are not copepods (shrimps of the genera Meganyctiphanes and Leptomysis) (Figure 6). Although Meganyctiphanes only occurs in low abundance, it was observed throughout the year. All recorded arthropod genera are part of the holoplankton. Most copepod genera (75%) were detected throughout the year and occurred with a high frequency of occurrence in all seasons (Figure 6A). The copepod genera Calocalanus, Centropages, Clausocalanus, Paracalanus and Oithona had the highest frequency of occurrence (100%) in all seasons, meaning that they were confirmed in all samples of the dataset (Figure 6A). In addition, among copepods, the genera Acartia, Temora and Mecynocera were also confirmed with high frequency (>80%) in all seasons. With 30 ASVs assigned to Calocalanus, the metabarcoding results were consistent with morphologically based copepod inventories describing a high species diversity of Calocalanus in the Adriatic (15 species, [27]). In contrast to the year-round genera, 36% of the copepod genera in the Southern Adriatic dataset were seasonal and were recorded only in autumn (Pontella, Pseudanthessius), winter (Pachos), spring/autumn (Neocalanus), summer/autumn (Diaixis and Copilia) and autumn/winter (Critomolgus). Season-specific genera were mostly recorded in lower read abundances. Pontella and Copilia were also recorded only sporadically and in lower abundances in previous Southern Adriatic zooplankton studies based on morphological identification [28,29]. Most of the copepod genera recorded by metabarcoding are known as common species of the Adriatic, living mainly in coastal and surface waters (Paracalanus, Clausocalanus, Calanus, Calocalanus, Acartia, Centropages, Oithona and Euterpina), but some were characterised as very rare copepods for the Adriatic (Pachos, Neocalanus) in previous studies [25,26,27]. In addition to free-living (pelagic or benthic) copepods, the Southern Adriatic dataset also enabled the identification of some marine invertebrate parasites/symbionts (Pseudanthessius) [78]. Our results also confirm that the Southern Adriatic dataset shows a similar copepod diversity as the dataset from the Tyrrhenian Sea, which was also analysed with 18S metabarcoding [15].

The division/phylum Cnidaria was represented in the dataset by the class Hydrozoa and Anthozoa (Figure 6A). While Anthozoa showed a low diversity with only two genera (Metridium and Tubastraea), the Hydrozoa comprised diverse orders; truly holoplanktonic (Siphonophorae and Trachymedusae) and hydrozoan medusae with obligate benthic life stagesthat form the meroplankton (Anthoathecata, Narcomedusae, Limnomedusae and Leptothecata). Together with copepods, hydrozoans were the most abundant group in the metazoan dataset. Hydrozoans accounted for 12% of all metazoan reads. In addition, most hydrozoan genera (Liriope (Limnomedusae), Clytia (Leptothecata), Sulculeolaria (Siphonophorae) and Solmundella (Narcomedusae)) were recorded throughout the year. The genus Aglaura (Trachymedusae) was only absent in summer. Obelia (Leptothecata) was characteristic for the winter–spring period, while Hydractinia (Anthoathecata) could only be detected in winter. All identified hydrozoans were previously recorded for the Southern Adriatic [34,36,37], with the exception of the two Siphonophorae genera, Nanomia (Siphonophorae), which was previously only identified for the Northern Adriatic [17,79], and Cordagalma (Siphonophorae), which represents a new record for the Adriatic but was previously recorded for the Mediterranean [80]. The results presented here, thus, clearly show that species from the genus Cordagalma enter the Adriatic. Nanomia (N. bijuga) was present most of the year, while Cordagalma (C. ordinatum) was only recorded in autumn (Figure 6). Meroplanktonic hydrozoans Liriope and Aglaura were also the most common and abundant hydromedusae in previous research for the area with juveniles present throughout the year [34].

All other metazoan groups (Annelida, Mollusca, Ctenophora, Phoronida, Bryozoa, Echinodermata, Nemertea, Platyhelminthes and Tunicata) together contributed 2% to the abundance of metazoans and were represented by a smaller number of genera. They also tended to occur sporadically and with lower seasonal abundance in the dataset (Figure 6A). Nevertheless, temporal patterns were identified for the metazoan community of the Southern Adriatic. In addition to copepods, the genera of polychaetes (Annelida), the meroplanktonic Polyophthalmus, the holoplanktonic Tomopteris, the holoplanktonic gastropod genus Creseis (Mollusca) and the pelagic tunicate genus Doliolum (Thaliacea) were recorded with high frequency of occurrence in all seasons (Figure 6A). The pelagic tunicate genus Oikopleura (Appendicularia) was absent from the Southern Adriatic metabarcoding dataset only in autumn. Similar to Tomopteris, another holoplanktonic polychaete genus (Typhloscolex) was recorded throughout the year, but Typhloscolex was absent in spring. The other polychaete genera in the Southern Adriatic dataset were mostly recorded only in winter–autumn, with the exception of Odontosyllis, a polychaete genus that was recorded in summer (Figure 6). These other polychaetes, which only occur at certain times of the year, make up the meroplankton community. Molluscs were generally present throughout the year, with holoplanktonic Creseis present year round and seasonal occurrences of meroplanktonic bivalve genera in winter (Barnea) and summer (Atrina and Hyotissa) or gastropods (Calliopaea and Haminoea) in summer/autumn.

The phyla Phoronida, Bryozoa, Platyhelminthes and Ctenophora were represented in the Southern Adriatic coastal dataset with the smallest number of genera (one or two) (Figure 6A), all of which were previously recorded for the studied area. The metabarcoding dataset record of Charistephane fugiens confirms the presence of this ctenophore in the Southern Adriatic, which was recently described as a new species for the Adriatic Sea [43]. This previous research characterised C. fugiens as being associated with the winter months, when there is a significant influx of Mediterranean water into the Southern Adriatic [43]. In our metabarcoding dataset, the annual occurrence of C. fugiens was also recorded only once, in spring (April), suggesting that a year-round population of this species has not yet been established/detected in the Southern Adriatic. However, further research, including a higher temporal and spatial resolution of the sampling strategies, will be necessary to clarify the distribution of the ctenophore C. fugiens distribution in the Southern Adriatic.

3.3. Network Analysis

In order to understand the complex interaction of the plankton community, an association network was created and analysed (Figure 7). Considering only associations with higher correlation (threshold 0.4), 72 genera composed the plankton network of the Southern Adriatic. Ten clusters (functional groups within the community) were identified within the network. The nodes (genera) in the clusters are strongly connected to each other but have few connections to genera outside their cluster. The clusters of the Southern Adriatic network counted between 12 and 2 nodes (genera) and were taxonomically diverse. The established connections between protists and metazoans confirmed the complexity of the interactions in the plankton community. The network clusters supported the seasonal patterns of co-occurrence of protist and metazoan genera (Figure 5 and Figure 6), and the identified hub nodes correspond to the genera present throughout the year. When considering betweenness centrality, the following seven genera (in descending order of betweenness centrality) were identified as hub nodes (keystone taxa) within the plankton community in the Southern Adriatic: Chaetoceros, Polyophthalmus, Zoothamnium, Meganyctiphanes, Metridium, Mecynocera and Temora (Figure 7). Considering the degree centrality measure, the following seven genera (in descending order of degree centrality measure) were identified as hub nodes: Lucicutia, Corycaeus, Chaetoceros, Meganyctiphanes, Metridium, Paracalanus and Acartia.

As the planktonic organisms share the same habitat and influence each other, the identified hub genera may comprise representative patterns of the plankton community [20], but the ecological role and importance of the associations established in the coastal Southern Adriatic network remain to be determined using a more comprehensive dataset. A comparison of the networks created for the future metabarcoding data will provide information on whether and how nodes (genera), clusters or the overall structure of the network change from year to year.

4. Conclusions

Our metabarcoding dataset provides the first comprehensive and united overview of the Southern Adriatic eukaryotic biodiversity patterns over the course of the year. Seasonal successions of the most important and highly diverse plankton groups were characterised, with the reliably identified genus/species level representatives characterised. In addition to a detailed description of the core, year-round eukaryotic plankton community, the metabarcoding approach allowed us to uncover a previously unknown diversity in the Southern Adriatic. This included taxa/genera ranging from pico-sized protist phytoplankton (Chloropicon and Prasinococcus), heterotrophic protist radiolarians (Sphaerozoum, Disolenia and Collozoum), parasitic dinoflagellates (Ellobiopsis) to metazoan Siphonophorae (Nanomia and Cordagalma). Taxa that occur consistently in a studied area represent a “core community” while taxa that occur only rarely form a “rare community”. For them, we must assume that our dataset may be unsuitable for characterizing seasonal patterns; characterisation and a sampling scheme with higher temporal resolution will address this issue. Nonetheless, the new findings suggest that there is great plankton diversity in this area that remains to be explored. In addition, the metabarcoding dataset of the Southern Adriatic contains a collection of ASVs that cannot currently be reliably identified to the species/genus level (due to a lack of reference sequences) but may make an important contribution to the ecology and biodiversity of the Southern Adriatic and Mediterranean. At the same time as characterising, the characterisation of taxonomic diversity is restricted by microscopy, and our metabarcoding dataset has successfully confirmed seasonal patterns and genus occurrences previously recorded by classical methods (microscopy) for the studied area. These results indicate a potential successful application of metabarcoding for marine monitoring, including the monitoring of invasive species and species not yet known. Importantly for monitoring applications, metabarcoding has the potential to reveal additional and previously unrecorded biodiversity. The results reported here show that a wide range of planktonic life can be recorded and analysed simultaneously, which is important for downstream analyses dealing with the intricate interactions in the planktonic network. The plankton association network provided a high-level insight into the global structure of the community, which needs to be further studied with more extensive spatiotemporal datasets. To summarise, our results highlight the added value of metabarcoding-based research for monitoring marine plankton. We have created the starting point to define the ecological context of the described plankton community of the Southern Adriatic and their interactions. In order to assess the consistency and ecological significance of the plankton diversity revealed by metabarcoding in this study, further observations are required in the future.