Sampling took place at two stations with long term sampling history and a continuous monitoring programme coordinated by the Centre for Marine Research (CIM), monthly during the years 2020 and 2021. Stations RV001 and RV004 (45°4'48''N, 13°36'36''E and 45°3'42.66''N, 13°32'56.976''E) are located 1 and 4 nautical miles off the Croatian coast and are used in routine monitoring for the collection of representative phytoplankton samples for the northern Adriatic Sea (Fig. 1). Samples for metabarcoding analysis were collected at each site using a phytoplankton net with a defined mesh size of 50 µm (Hydrobios, Germany). Samples were collected by vertical net hauls from ∼ 30 m depth to the surface to cover species from the entire water column. Concentrated samples were collected in 500 ml plastic bottles and filtered on 1.2 μm cellulose ester membranes (47 mm, Whatman, UK) until complete filter saturation (total volume ranging from 20 ml up to 150 ml). Filters were stored in small Petri dishes, at -80°C as suggested by Baricevic et al. (2022) until further lab processing.

Figure 1.  

Map of two sampling sites located near Rovinj with long sampling history used as reference points for sample collection.

Figure 2.  

Contribution of phytoplankton groups. Changes in contribution in genetic diversity of different phytoplankton groups depending on clustering method used (upper plot) and number of assigned species of each phytoplankton group depending on clustering method used.

Figure 3.  

Dendrogram of diatom genera. Dendrogram showing species as listed in Viličić et al. (2002) as recovered from all published and observed data until 2002, representing a potential diatom diversity observable in the eastern northern Adriatic. Genetic diversity acquired in this study is compared to genera found in the checklist (2002). Species recovered from metabarcoding data, but so far absent from the 2002 checklist, were added to the dendrogram: Filled purple dots at the tip points indicate diatom genera obtained with metabarcoding, while orange dots present metabarcoded genera that were not present in a checklist. The bubble plot shows the number of morphospecies from the checklist (column Species) and the number of different clustering units assigned per genera depending on the clustering method (columns ASVs, 99%-OTUs, 97%-OTUs) and removal of low abundant clustering units (columns ASVs*, 99%-OTUs*, 97%-OTUs*). The bar plot on the right showing the contribution of species/clusters for each genus/clustering method. Genera not found in metabarcoding data were left out. The bottom bar plot indicates the number of morphospecies (first column) and clustering units (other columns) generated by each clustering method for the diatom group.

Total DNA was extracted from stored filters using either Gentra Puregene DNA Isolation Kit (Qiagen, Germany) or NucleoSpin eDNA Water kit (Macherey-Nagel, Germany), both according to the manufacturer′s protocol. Isolated genomic DNA was eluted in 100 µl per sample and obtained concentrations were quantified by measuring DNA concentration/absorbance spectra using the Nanophotometer spectrophotometer (Implen, Germany). Samples were stored at 4°C until further processing. DNA metabarcoding library preparation and sequencing were carried out by AllGenetics & Biology SL (www.allgenetics.eu). DNA concentration was quantified for each extract using the Qubit High Sensitivity dsDNA Assay (Thermo Fisher Scientific). For DNA metabarcoding library preparation, eukaryote specific primers TAReuk454FWD1 (5′-CCAGCA(G/C)C(C/T)GCGGTAATTCC-3′) and TAReukREV3 (5′-ACTTTCGTTCTTGAT(C/T)(A/G)A-3′) (Stoeck et al. 2010, Piredda et al. 2017) were used, targeting around 430 bp of the V4 hypervariable region of the small subunit ribosomal RNA gene with the Illumina sequencing primer sequences attached at their 5’ ends. In the amplification step, PCRs were conducted in a final volume of 12.5 μl, containing 1.25 μl of template DNA optimised by diluting the starting template DNA, 0.5 μM of the primers, 3.13 μl of Supreme NZYTaq 2x Green Master Mix (NZYTech) 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 were attached in a second amplification step with identical conditions, but only 5 cycles and 60ºC as the annealing temperature (Vierna et al. 2017). 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 Biotek), 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) 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) aiming for a total output of 4 gigabases. Illumina paired-end raw data consisted of forward (R1) and reverse (R2) reads and were stored separately sorted by library, as fastq files.

The CIM Phytoplankton 18S ribosomal RNA reference database (CIMPhy18) was assembled using ARB software v.7.0 (Ludwig et al. 2004) by merging downloaded phytoplankton sequence records from the Silva (Quast et al. 2013) and Protist Ribosomal Reference (PR2) (Guillou et al. 2013) databases. Phytoplankton groups included Dinophyceae (dinoflagellates), Ochrophyta containing Bacillariophyceae (diatoms) and Chrysophyceae (golden algae), Cryptophyceae (cryptophytes), Haptophyta (mostly coccolithophorids) and Chlorophyceae (green micro algae). First, the Silva-aligned database was imported into ARB software with all additional data and NCBI accession numbers linked to phytoplankton taxa were compared between Silva and PR2 databases. PR2 accession number entries that were not present in the Silva database were downloaded from Entrez Molecular Sequence Database System (www.ncbi.nlm.nih.gov/Web/Search/entrezfs.html) in genbank (.gb) format and imported into ARB. Phytoplankton groups/phyla chosen for the reference database assembly were aligned in ARB Sequence Editor and a reference database phylogenetic tree was constructed. Relevant in-house sequences from the in-house CIM phytoplankton culture collection, representing barcoded sequences from northern Adriatic common phytoplankton taxa, were used to extend the database and to include species not included in the standard release of the databases. To acquire information about the taxonomic completeness of the constructed database, phytoplankton taxa entries were compared to: (1) the eastern Adriatic phytoplankton checklist (Vilicic et al. 2002) and (2) the list of northern Adriatic phytoplankton species from long-term monitoring data. Both lists were obtained from long-term studies with light microscopy. Before comparison, taxonomic data were normalised using the Species matching tool from GBIF (www.gbif.org) to ensure accurate classification and taxonomic uniformity. Sequences for missing taxa were searched in NCBI (www.ncbi.nlm.nih.gov), curated and, if suitable, added to the CIMPhy18 database along with accompanying metadata. Sequences for the reference database were extracted in fasta format and their accompanying taxonomy in a tax file, prepared to be imported into Mothur (Schloss et al. 2009).

Illumina paired-end reads for the small subunit region were processed using a custom script for the Mothur pipeline v.1.47.0 (Schloss et al. 2009), based on MiSeq standard operating procedure (https://mothur.org/wiki/miseq_sop/). Contigs were assembled from fastq files containing forward and reverse Illumina MiSeq reads and trimmed to the overlapping section, with no ambiguous bases allowed; the maximum homopolymer size was 8 bp and the maximum tolerated sequence length 450 bp. Sequences were screened for chimeras using the VSEARCH command (Rognes et al. 2016), removing sequences with chimeras from analysis from their group only (dereplicate=t). Clustered sequences are reported as amplicon sequence variants, ASVs (Callahan et al. 2017), with one nucleotide tolerance (cutoff = 1) eliminating erroneous sequences formed due to sequencing mistakes and operational taxonomic units, OTUs were assembled, based on fixed sequence 97%- and 99%-identity (99%-OTUs and 97%-OTUs) threshold. Both times column-format distance matrices were calculated by default option (one gap) and were assigned to OTUs using OptiClust algorithm (Westcott and Schloss 2017). The consensus classification for each ASV and OTU clustering units was performed using a naïve Bayesian classifier (Wang et al. 2007) trained using the CIMPhy18 constructed reference database, with a 90% bootstrap confidence threshold. Two ASV and OTU tables were built: one containing all informations, regarding read number present in a cluster (ASVs, 97%-OTUs or 99%-OTUs) and the second, including only clustering units represented with five of more reads per cluster (ASVs*, 97%-OTUs* or 99%-OTUs*). Species scientific names and corresponding authors were checked according to AlgaeBase (www.algaebase.org) prior to checklist formation. All downstream analyses were performed using R Statistical Software (v.4.2.3; R Core Team 2023) and data visualisations were made using ‘ggplot2’ (Wickham 2009) and ‘ggtree’ (Yu et al. 2017) Rpackages and assembled using ‘aplot’ (Yu 2022) Rpackage.