Phytoplankton Diversity of a Natural Karst Lake Combining Morphological and Molecular Approaches

Abstract

Phytoplankton are considered to be one of the most sensitive indicators of the ecological status of lakes. Nowadays, it is essential to recognize the prospects of the molecular approach (eDNA metabarcoding) in phytoplankton community assessments and combine them with the existing traditional microscopy-based morphological approach before its standardization. In this study, the aim was to characterize the phytoplankton community of a natural karstic lake by combining and comparing the morphological and molecular approach to check the applicability of eDNA metabarcoding as a biomonitoring tool. A total of 51 phytoplankton taxa were found using the morphological approach, whilst the molecular approach discovered 97 ASVs that corresponded to the algal community. The comparability of both approaches in describing phytoplankton communities is evident in the designation of centric diatoms, dinoflagellates and cryptophytes as descriptive taxa. Furthermore, both approaches proved reliable in detecting functional groups (Lo, C, X2, X3) with similar ecological demands. Moreover, the results have shown that euphotic zone samples can be reliably exchanged by composite samples to provide an accurate characterization of phytoplankton communities in the euphotic zone. It was confirmed that eDNA metabarcoding is an applicable tool for biodiversity monitoring of a natural karst lake and should be used as a feasible supplement to traditional microscopy in the phytoplankton community assessments, with regards to the drawbacks of each method.

1. Introduction

Natural lakes in Croatia are a phenomenon mainly associated with the karst landforms of the Dinaric ecoregion. The Dinaric ecoregion is a part of the Mediterranean Basin, well recognized as one of the Earth’s biodiversity hotspots [1]. Each karst lake is a unique freshwater ecosystem with its own geological, physical and chemical characteristics [2]. Looking from this perspective, biodiversity protection, conservation and sustainable management of freshwater karstic ecosystems require a set of applicable practices based on the best available technologies and establishment of relevant measures founded on up-to-date information and a quality knowledge database.

Phytoplankton play an essential role as the foundation of the food web and as primary producers ubiquitous in all aquatic ecosystems, being especially dominant in the pelagic zone of lakes [3]. As a key Biological Quality Element prescribed by the Water Framework Directive [4], phytoplankton are considered to be one of the most sensitive indicators of a lake’s ecological status because they respond rapidly to any environmental changes [5,6]. Phytoplankton biotic metrics used to assess the ecological status of surface water bodies are commonly constructed on traditional morphology-based microscopic identification of taxa, which is time-consuming, costly and has limitations in terms of reproducibility, comparability and applicability in biomonitoring programs [7,8]. Compared to the morphological approach, environmental DNA (eDNA) metabarcoding has been recognized as a tool with the potential to revolutionize biomonitoring and bioassessment, as it offers numerous advantages that include higher accuracy in species identification, increased detection of cryptic diversity and genetic variability, as well as high automation potential including high spatial and temporal resolution [7,9,10]. However, eDNA metabarcoding suffers from several setbacks, one being the incompleteness of the current reference database, which could potentially be compensated for by adding representative sequences of local species from the studied aquatic ecosystems [8]. This method may also be less suitable for estimating abundance, may not provide information on the age or size structure of a population [10] and needs to be standardized before it can be adequately applied in routine monitoring [11]. Nevertheless, the combination of microscopy and molecular methods can provide great improvements in phytoplankton community assessments [7,8,12]; therefore, it is important to recognize the potential of eDNA-based methods and harmonize them with the existing traditional morphological approach.

The aim of this study was to gain a better insight into phytoplankton diversity in the natural karst Lake Visovac (Krka River, southern Croatia) in order to achieve an adequate implementation of a reliable ecological status assessment. Specifically, the aims were to: (a) describe the horizontal and vertical distribution of phytoplankton in the Lake Visovac by applying the traditional morphological approach and the molecular eDNA metabarcoding approach using the amplicon sequencing of hypervariable region V9 of the 18S rRNA gene to investigate total eukaryotic phytoplankton diversity, (b) determine if composite samples could be used as a relevant substitute for discrete sampling in the characterization of phytoplankton community in the euphotic zone, (c) compare the results obtained by morphological and molecular approaches, and (d) establish the applicability of eDNA metabarcoding as a biomonitoring tool in Lake Visovac. The study results will enable establishing a functional system for monitoring changes in the environment and a contribution to the protection of unique karst ecosystems such as the Krka River.

2. Materials and Methods

2.1. Study Area

Krka River rises at the base of Dinara Mountain near the city of Knin in Croatia. It is a 72.5 km long karstic river situated in the central part of the eastern Adriatic coast. Its course is distinguished by alternating lotic and lentic parts as well as tufa deposits forming barrages and cascades. Lake Visovac has a volume of 103 × 10⁶ m³ and origins from the post-Würm period with the formation of the final and the largest tufa barrier in the Krka River hydrosystem, named Skradinski Buk [13]. Following the provisions of the national typology, Lake Visovac is classified as a medium-sized, medium-depth lowland lake on the carbonate substrate [14].

2.2. Sampling and Methods

Phytoplankton composition and biomass were investigated during August 2018 on 10 sampling stations (V1 to V10) along the limnetic and littoral zone of Lake Visovac (Figure 1). The study area extended from the northern part of the Lake close to Roški slap waterfall (V1) to the southern part immediately before Skradinski Buk waterfall (V8) and near the confluence of the tributary Čikola River (V9 and V10). Stations were divided into three groups according to their geographic position on the Lake: upper (V1, V2), central (V3, V4, V5 and V6) and lower stations (V7, V8, V9 and V10). Station V8 was excluded from the analyses due to shallow maximum depth (3 m).

The total number of samples was defined by the maximum depth of each station, Secchi depth, thermocline depth and mixing depth of the water column with respect to the temperature difference. Water column transparency (Z_SD) was determined with a Secchi disc and used for the calculation of euphotic zone depth (Z_EU) by multiplying with a standardized factor (2.5 × Secchi depth) for the Mediterranean geographical region [15]. The biological and chemical water samples were collected using the vertical sampler (Hydro-Bios Apparatebau Gmbh, Altenholz, Germany). For the morphological analysis of phytoplankton, discrete samples were taken at 5 m depth intervals (from the surface to the bottom) together with the composite samples taken from the euphotic zone on all stations. With respect to the euphotic zone depth, the discrete samples were divided into those from the euphotic and those from the aphotic zones, and means were calculated for each station, which were used in all further analyses. Grouping of discrete samples into euphotic zone and aphotic zone sets was also applied to spatially compare the phytoplankton community across sampling stations in Lake Visovac. For the DNA analysis of phytoplankton, composite and aphotic zone samples were taken at all stations.

Composite samples for chemical analysis of water were collected together with the phytoplankton samples and stored at −20 °C until laboratory processing. Chemical analysis included quantification of nitrate (NO₃⁻-N), nitrite (NO₂⁻-N), ammonium (NH₄⁺-N), total nitrogen (TN-N), orto-phosphate (PO₄³⁻-P), total silica (SiO₂), total inorganic carbon (TIC), dissolved inorganic carbon (DIC), total organic carbon (TOC), dissolved organic carbon (DOC) and bicarbonates (HCO₃⁻) using standardized methods [16].

Phytoplankton samples were placed into 250 mL volume plastic bottles, preserved with formaldehyde solution (2%) on the field and stored in the dark at 4 °C. Phytoplankton biomass was determined according to the Utermöhl method [17] using a Zeiss AxioVert inverted microscope equipped with an AxioCam MRc camera (Carl Zeiss AG, Jena, Germany). Taxa identification was performed using relevant literature [18,19,20,21,22,23,24] and names were assigned according to Algaebase [25]. Images of species were processed using the program AxioVision LE 4.8 (Carl Zeiss AG, Jena, Germany). The species were allocated into appropriate functional groups (FGs or coda) following the relevant literature [5,26]. All sampling and analytical procedures were performed according to following standards: HRN EN ISO 5667-3:2018, HRN EN 15204:2008, HRN EN 16695:2015 [27,28,29].

2.3. DNA Isolation

Samples for DNA extraction were filtered on Nucleopore track-etched polycarbonate membrane filters (47 mm diameter, 0.2 µm pore size; Whatman International Ltd., Maidstone, UK) in a volume of approximately 300 to 400 mL, depending on the suspended particles in the sampled water. After filtration, the filters were stored at −20 °C until further processing. Filters were cut into smaller pieces for DNA extraction using the DNeasy PowerWater Kit (Qiagen, Hilden, Germany). The manufacturer’s instructions were followed for isolation, with a minor change in the final step, where 60 µL of sterile DNA-free PCR-grade water was added instead of Qiagen’s C6 Solution. The quality of extracted DNA was assessed with a NanoDrop spectrophotometer (BioSpec—nano, Shimadzu Corporation, Kyoto, Japan).

2.4. PCR and Bioinformatic Processing

The hypervariable V9-region of the SSU rRNA gene (ca. 130 bp) was amplified using the universal eukaryotic primer pair according to the protocol of Stoeck et al. [30]. The primers used were 1391F (5′-GTACACACCGCCCGTC-3′) and EukB (5′-TGATCCTTCTGCAGGTTCACCTAC-3′), designed by Amaral-Zettler et al. [31]. The Polymerase chain reactions (PCR) program included the initial step at 98 °C for 30 s, 30 cycles of 94 °C for 30 s, 57 °C for 45 s and 72 °C for 30 s, with the final elongation step at 72 °C for 5 min [32]. PCR products were assessed by visualizing on a 1% agarose gel. Sequencing libraries were prepared using the NEB Next^® Ultra™ DNA Library Prep Kit for Illumina (NEB, USA). Libraries were sequenced on an Illumina MiSeq platform, generating 250-bp paired-end reads (SeqIT GmbH & Co. KG, Kaiserslautern, Germany).

The raw Illumina reads for the V9-region were demultiplexed using Cutadapt v3.0 [33], removing the barcodes in the 5′ to 3′ combination. Subsequently, the quality of the demultiplexed raw sequencing was checked using the FastQC tool [34]. After the initial steps, reads were processed using the QIIME2-2020.11 pipeline using the following steps: importing and demultiplexing raw sequencing data, then quality filtering and denoising using the DADA2 plugin [35] for correcting Illumina sequencing amplicon errors. Reads were trimmed at 5′ end for 20 bp (primer removal) and truncated to 185 nucleotides to remove the last, poor-quality nucleotides. Taxonomic assignment of the resulting amplicon sequencing variants (ASVs) was performed using the Naïve Bayes classifier. The Naïve Bayes classifier was pretrained on the Protist Ribosomal Reference (PR2) database v.4.14.0 [36] with a 99% OTU identity threshold. Metazoa sequences were filtered out from the dataset by using taxa filtering. A phylogenetic tree was created from the filtered taxa tables to support the phylogenetic diversity metrics used in the q2-diversity plugin. The raw sequence reads are deposited in the European Nucleotide Archive (ENA) under the project number PRJEB60049.

Ochrophyta, Dinoflagellata, Cryptophyta and Fungi accounted for the majority of eukaryote reads (Figure S1), and unassigned eukaryotes were excluded from graphical presentation. Only ASVs taxonomically assigned to the phytoplankton community were filtered from the main data and used for all further statistical analyses.

2.5. Statistical Analyses

All statistical multivariate analyses were carried out in PRIMER v7 for Windows (Primer-E Ltd., Plymouth, UK). Principal Component Analysis (PCA) was performed to outline and visualize the relationships between environmental variables. One-way analysis of similarity (ANOSIM) was used to test significant differences between composite and discrete samples and determine if composite samples can be used as a relevant substitute for discrete sampling in the characterization of phytoplankton communities in the euphotic zone. Non-metric Multidimensional Scaling (NMDS) was carried out to determine the spatial patterns in the phytoplankton community structure with respect to the sampling stations. Prior to all analyses, the data were normalized using log-transformation. Graphical charts were created in Microsoft Excel for Microsoft 365 (Microsoft Corporation, Redmond, WA, USA). The map of the study area was created using the Free and Open Source QGIS 3.16 software [37].

3. Results

3.1. Physical and Chemical Parameters

The environmental variables of water measured on nine sampling stations for composite samples are presented in Table 1. Secchi depth ranged from 2.5 m (station V2) to 5 m (station V9). Stations with the highest maximum depths were V7 and V9 (24 m and 23 m, respectively). The shallowest station was V10 (5 m) with the euphotic zone extending along the entire water column. The temperature of water ranged from 20.6 °C on station V7 to 26.2 °C on station V10. The lowest O₂ concentration was recorded on station V4 (9.03 mg L⁻¹) and the highest on station V1 (11.68 mg L⁻¹). Oxygen saturation was in the range from 102.6% (station V4) to 138% (station V1). The lowest pH was recorded on station V7 (7.75), whilst the highest was on station V1 (8). The electrical conductivity of water ranged from 494 µS cm⁻¹ on station V10 to 562 µS cm⁻¹ on station V6 (Table 1).

The environmental variables of water measured from discrete-depth samples in the euphotic and aphotic zone on nine sampling stations are presented in Table S1. The minimum value of water temperature in the euphotic zone was 22.4 °C, measured on station V2, and the maximum was 26.4 °C on station V10. The lowest concentration of O₂ was recorded on station V10 (10.13 mg L⁻¹), whilst the highest concentration and saturation of O₂ in the euphotic zone were measured on station V7 (12.28 mg L⁻¹ and 147.5%, respectively). The lowest saturation of O₂ was observed on station V5 (80.1%). The lowest pH (7.88) was recorded on station V2, whilst the highest (8.20) was on station V5. The electrical conductivity of water in the euphotic zone ranged from a minimum of 495 µS cm⁻¹ on station V10 to a maximum of 562 µS cm⁻¹ on station V2. In the aphotic zone, water temperature varied from 16.5 °C to 20.8 °C (stations V7 and V1, respectively). The lowest concentration and saturation of O₂ were recorded on station V9 (2.91 mg L⁻¹ and 27.4%, respectively) and the highest on station V1 (9.60 mg L⁻¹ and 108.2%, respectively). Value of pH in the aphotic zone ranged from 7.05 to 7.80 (stations V4 and V5, respectively), whilst electrical conductivity of water varied from 517 µS cm⁻¹ to 587 µS cm⁻¹ (stations V9 and V1, respectively).

Chemical parameters of water are presented in Table S2. The concentration of NO₃⁻-N was very low on most stations (<0.1 mg L⁻¹), with slightly higher values on stations V4 (0.7 mg L⁻¹) and V5 (0.3 mg L⁻¹). Very low concentrations of NO₂⁻-N and NH₄⁺-N were recorded on all stations (<0.001 mg L⁻¹ and <0.01 mg L⁻¹, respectively). The concentration of PO₄³⁻-P varied from the lowest measured on six stations in total (<0.01 mg L⁻¹) to the highest measured on station V9 (0.6 mg L⁻¹). The values of SiO₂ ranged from the lowest on station V9 (1.5 mg L⁻¹) to the highest on station V6 (4.5 mg L⁻¹). The highest concentration of TN was measured on stations V4 and V5 (1 mg L⁻¹), while it was very low on the other eight stations (<1 mg L⁻¹). The lowest concentrations of inorganic carbon compounds were detected on station V3 (TIC and DIC of 11.48 mg L⁻¹ and 10.32 mg L⁻¹, respectively), while the highest was present on station V7 (TIC and DIC of 13.98 mg L⁻¹ and 13.83 mg L⁻¹, respectively). TOC was in the range between 1.55 mg L⁻¹ (station V3) and 2.39 mg L⁻¹ (station V7). The values of DOC ranged from 0.80 mg L⁻¹ (V2) to 2.37 mg L⁻¹ (V1), whilst HCO₃⁻ ranged from 151 mg L⁻¹ (V4) to 212 mg L⁻¹ (V7).

Principal component analysis (PCA) performed for the 14 environmental variables explained 59% of the total variance on the first two PC axes (Table S3). NO₂⁻-N and NH₄⁺-N were excluded from PCA because their concentrations did not vary across the stations. The most important parameters for PCA axis 1 were temperature, DIC and electrical conductivity (intra-set correlations: 0.433, −0.361 and −0.338, respectively). Regarding axis 2, O₂ concentration, nitrate and TOC were the variables with the most weight for ordination (intra-set correlations: −0.464, 0.369 and −0.368, respectively). PCA arranged samples (Figure 2) into four groups: the first group consisted of samples from the central stations (V3, V4, V5, V6) and the upper station sample V2, whilst the second group included samples from the lower part of the Lake (V9 and V10). Sample V1 from the uppermost part and sample V7 taken from the lower part of the Lake were singled out.

3.2. Characterization of Phytoplankton Community According to Morphological Approach

Based on the morphological approach, a total of 51 phytoplankton taxa were found during the research period. Identified taxa belonged to eight major groups: Chlorophyta (23), Bacillariophyta (8), Ochrophyta (7), Cyanobacteria (4), Charophyta (4), Cryptophyta (2), Miozoa (2), Euglenozoa (1). In total, six taxa contributed more than 5% of the total biomass. The main descriptive phytoplankton species was dinoflagellate Ceratium hirundinella (O.F. Müller) Dujardin, followed by centric diatom Pantocsekiella ocellata (Pantocsek) K.T. Kiss & E. Ács, cryptophytes Cryptomonas sp. and Plagioselmis nannoplanctica (H. Skuja) G. Novarino, I.A.N. Lucas & S. Morrall, chlorophyte Tetraselmis cordiformis (N. Carter) Stein, and dinoflagellate Parvodinium inconspicuum (Lemmermann) Carty (Table S4).

In terms of biomass share, the most dominant group in the composite samples on stations V1 and V2 (Figure 3) was Miozoa (54% and 48%, respectively), followed by Chlorophyta (24% and 29%, respectively) and Cryptophyta (14% and 12%, respectively). A corresponding situation regarding the phytoplankton assemblage was observed in the euphotic zone samples on the same stations, with a slight difference in shares of major taxonomic groups. Cryptophyta emerged as a dominant group (41%) in the composite sample on station V3, whilst Bacillariophyta and Miozoa were subdominant (25% and 18%, respectively). Conversely, the euphotic zone samples on station V3 were dominated by Miozoa (33%), followed by Bacillariophyta and Cryptophyta (26% and 18%, respectively). Phytoplankton community in the composite sample on station V4, situated close to the littoral zone in the central part of Lake Visovac, was characterized by the dominance of Miozoa with 68% of the total phytoplankton biomass, followed by subdominant groups Bacillariophyta (13%) and Chlorophyta (13%). The most dominant group in the euphotic zone samples on station V4 was Bacillariophyta (30%), accompanied by the subdominant groups Chlorophyta and Cryptophyta (28% and 23%, respectively). Cryptophyta and Bacillariophyta were the main groups characterizing the community in both composite and euphotic zone samples on station V5, positioned in the limnetic zone of the central part of the Lake (36% and 34%; 27% and 29%, respectively). Miozoa was the most dominant group on station V6 in both composite and euphotic zone samples (46% and 43%, respectively), with Bacillariophyta and Cryptophyta as subdominant groups. Bacillariophyta dominated in the composite sample on station V7, followed by Cryptophyta (44% and 30%, respectively), whilst Miozoa dominated the assemblage (50%) in the euphotic zone sample. Station V9, placed in the lower part of the Lake, was characterized by a complete dominance of Miozoa in the composite sample, reaching 80% of the total phytoplankton biomass, whilst in the euphotic zone sample, this share was 49%. Cryptophyta appeared as a dominant group in the composite sample on station V10 (46%). However, in the euphotic zone on the same station, Miozoa outnumbered Cryptophyta in their biomass share (36% and 30%, respectively).

Regarding the aphotic zone samples, stations V1 and V3 were co-dominated by Chlorophyta and Miozoa (39% and 34%; 41% and 43%, respectively). Chlorophyta was also a dominant group on station V2, followed by Miozoa (50% and 29%, respectively). In the aphotic zone samples on stations V4, V5 and V7, the domination of Miozoa was observed (61%, 57% and 64%, respectively). The aphotic zone on station V6 was characterized by Cryptophyta and Miozoa (30% and 28%, respectively), whilst on station V9, Chlorophyta took over domination (51%). Ochrophyta were present in all samples in the range from 1% (composite sample on station V6) to 11% (composite sample on station V10). There was no aphotic zone on station V10 due to the shallow maximum depth of 5 m (see Table 1).

The descriptive Reynolds’ functional groups on upper stations V1 and V2 were Lo and X2 in both composite and euphotic zone samples (Figure 4). The composite and euphotic zone samples on station V3 were characterized by functional group X2 with the highest biomass share (50% and 37%, respectively), followed by associations C and Lo (25% and 18%; 26% and 33%, respectively). Functional group Lo dominated the assemblage in the composite sample on station V4, with 68% of the total phytoplankton biomass, whereas associations C and X2 were both contributing with 12% to the total phytoplankton biomass. Meanwhile, the functional group X2 dominated the assemblage (51%), with coda C and Lo being subdominant (30% and 15%, respectively) in the euphotic zone samples on station V4. The functional group X2 was dominant, whilst codon C appeared as subdominant in both composite and euphotic zone samples on station V5. Group Lo again became the most dominant in both samples on station V6, followed by coda X2 and C. Codon C prevailed in the composite sample on station V7, together with association X2 (43% and 39% of the total phytoplankton biomass, respectively), whilst in the euphotic zone samples, codon Lo gained dominance (50%). Codon Lo dominated the assemblage with 80% of the total phytoplankton biomass in the composite sample on station V9, as was also the case in the euphotic zone samples, although with a lower share (49%). The descriptive functional group in the composite sample on station V10 was X2 (57% of the total phytoplankton biomass), whereas in the euphotic zones sample, it co-dominated with codon Lo. Species belonging to codon E (genus Dinobryon) were present in a very low biomass share on every station, except for the somewhat higher shares noted in the composite sample on station V10 and the euphotic zone sample on station V9 (11% and 9%, respectively).

In the aphotic zone samples, stations V1, V2 and V6 were dominated by codon X2 (51%, 64% and 45% of the total phytoplankton biomass, respectively), followed by codon Lo (34%, 29% and 28%, respectively). Functional groups X2 and Lo co-dominated the phytoplankton assemblage in the aphotic zone on station V3 (47% and 43%, respectively). On stations V4, V5 and V7 the descriptive codon was Lo (61%, 57% and 64%, respectively), followed by codon X2 (33%, 37% and 28%, respectively). Functional group X2 clearly dominated the assemblage on station V9 (77% of the total phytoplankton biomass).

The NMDS analysis of phytoplankton taxonomic composition according to the morphological approach (Figure 5) indicated segregation of almost all discrete-depth samples into two separate groups, the first one comprising the euphotic zone and the second one including the aphotic zone. The samples from station V10 were all counted into the euphotic zone.

Considering the horizontal distribution, NMDS analysis of the composite samples showed grouping of samples from limnetic and littoral stations (Figure 6). Cryptomonas sp. and P. ocellata dominated on central limnetic stations V3 and V5. Stations from the littoral zone were clustered into two distinct groups as follows: the first group, dominated by C. hirundinella (stations V1, V2, V4, V6, V9), and the second group, dominated by P. ocellata (station V7) and Cryptomonas sp. (station V10).

The results of one-way ANOSIM pairwise tests (Table 2) have demonstrated significant differences when comparing composite samples vs. aphotic zone samples (p = 0.042), and euphotic vs. aphotic zone samples (p = 0.006). On the other hand, the negative R statistic close to zero and the high significance level indicated very low differences between composite samples and euphotic zone samples.

3.3. Characterization of Phytoplankton Community According to Molecular Approach

ANOSIM pairwise testing (Table 2) confirmed that the composite samples can be used as representative of the phytoplankton community in Lake Visovac and further characterized by molecular analyses. Correspondingly, eDNA metabarcoding analyses were performed on composite and aphotic zone samples. A total of 1,010,539 quality reads were yielded in 19 samples for eukaryotes, featuring 7140 ASVs at the 99% similarity level. After the Metazoan sequences were filtered out, a total of 902,798 quality reads with 7002 ASVs were found. The mean frequency per sample was 47,515 (min. 9779; max. 108,716). Only 97 ASVs corresponded to algal community with a total of 150,005 quality reads (Table S5).

Based on the sum of the relative abundance of the family ranks determined using the molecular approach, the taxonomically unassigned phytoplankton ASVs were found to have the highest relative abundance in all recorded composite and aphotic zone samples, with approximately 50% recorded in composite sample V7 (Figure 7). The four most abundant family ranks were polar centric diatoms (Mediophyceae), cryptophyta (Cryptomonadales), chrysophytes (Chrysophyceae_X), and dinoflagellates (Peridiniales). Interestingly, the two most abundant family ranks had the highest relative abundance in aphotic zone samples: 38% in V6 for Mediophyceae and 58% in V9 for Cryptomonadales. Furthermore, two other most abundant family ranks had the highest relative abundance in composite samples in V2 (31%) for Chrysophyceae_X and in V10 (67%) for Peridiniales. The families with total relative abundance between 65% and 20% in all samples were: Chlorodendrales, Suessisales, Pyrenomonadales, Katablepharidales, Cryptophyceae_X, Dictyochophyceae_X and Pseudodendromonadales. Of these family ranks, four had the highest relative abundance in the composite samples as follows: Suessisales (V1), Dictyochophyceae_X (V4), Catablepharidales (V6), and Pseudodendromonadales (V9). In contrast, Cryptophyceae_X (V2), Chlorodendrales (V3), and Pyrenomonadales (V7) had the highest relative abundance in the aphotic zone samples, while their relative abundance in the composite samples was about or less than 1%. The total relative abundance of the two family ranks corresponding to the Sphaeropleales and Synurales was less than 20% in all samples, but their abundance was highest in the aphotic zone samples (V3 for Sphaeropleales and V4 for Synurales). Perkinsida_X, Chrysochromulinaceae and Bicoecales corresponded to families with a total relative abundance in all samples equaling less than 10%. Families with a proportion of less than 1% were Gymnodiniaceae, Dolichomastigales, Cyanidiales, Charophyceae_X, Chlorellales, Distigmidae, Euglenophyceae, Petalomonadida, Choanoflagellida_X, Nucleariida, Araphid-pennate, Bacillariophyta_X, Raphid-pennate, Eustigmatophyceae_X and Xanthophyceae_X.

The descriptive functional groups in the eDNA composite sample of station V1 were Lo (47%) and C (30%) (Figure 8). In the composite eDNA sample of station V2, the dominant association was group C (with 34%), whereas coda X3, X2 and Lo were subdominant (with 29%, 18% and 13%, respectively). Functional groups C and Lo dominated the assemblage on station V3 (37% and 34%, respectively), whilst X3 and X2 were subdominant (14% and 10%, respectively). Codon X3 and X2 emerged as dominant on the station V4 (45% and 34%, respectively), whilst codon C appeared as a subdominant (11%). Functional groups C and X2 were the most dominant on stations V5 (34% and 32%, respectively) and V6 (33% and 32%, respectively), whilst codon X3, Lo and J were subdominant. Functional group X3 dominated the assemblage again on station V7 (53%), followed by functional groups X2, C and Lo (16%, 15% and 7%, respectively). Codon X2, X3 and Lo prevailed on station V9 (with 31%, 26% and 21%, respectively). At station V10, codon Lo showed clear dominance (73%).

Regarding the eDNA aphotic zone samples, functional group X2 prevailed with over 50% of total phytoplankton biomass on the majority of studied stations (V3 with 74%, V4 with 62%, V5 with 52%, V7 with 72%, and V9 with 68%). Station V1 was characterized by the dominance of codon C (36%) with coda Lo and X2 as subdominant (27% and 22%, respectively), whilst coda X2 and C had the highest biomass share on station V2 (38% and 30%, respectively). Aphotic zone samples from station V6 were characterized by functional groups C, X2 and X3 (38%, 25% and 23%, respectively).

NMDS analysis of phytoplankton eDNA metabarcoding has disclosed two groups, the first one including composite samples and the second one comprising aphotic zone samples (Figure 9).

3.4. Morphological and Molecular Diversity of Phytoplankton in Lake Visovac

The highest number of phytoplankton taxa in the composite samples based on the morphological approach (Figure 10a) was detected on station V10 (26) and the lowest on station V1 (16). The Margalef Richness Index in the composite samples ranged from 0.99 on station V1 to 1.79 on station V10 (Figure 10b). Pielou’s Evenness Index (Figure 10c) in the composite samples varied between 0.39 (station V4) and 0.78 (station V9). Concordantly, the lowest Shannon–Wiener Diversity Index value (Figure 10d) was present on station V4 (1.22), while the highest was on station V9 (2.38).

The highest number of phytoplankton ASVs provided by eDNA metabarcoding in the composite samples (Figure 10a) was detected on station V10 (49), whilst the lowest on station V6 (15). The lowest Margalef Richness Index in the eDNA composite samples (Figure 10b) was obtained on station V6 (2.04), whilst the highest taxa richness was recorded on station V10 (4.76). Results obtained by eDNA metabarcoding showed higher values for both indices compared to morphological approach data, with the exception of station V6. Pielou’s Evenness Index in the eDNA composite samples (Figure 10c) varied between 0.41 (station V10) and 0.82 (station V6). Compared to morphological composite samples, Pielou’s Index was higher in the eDNA composite samples on stations V2, V4, V6, V8 and V9 but lower on stations V1, V3, V7 and V10. The Shannon–Wiener Diversity Index (Figure 10d) for eDNA composite samples varied between 1.61 (station V10) and 2.86 (station V9), and provided higher values than morphological approach data, except for station V10.

The total number of taxonomically assigned ASVs detected by the molecular approach was two times higher than the number of taxa revealed by traditional morphological identification using microscopy. Considering the results on each station, the number of ASVs was higher than taxa number on every station except for station V6, which proved to be an outlier just as for morphological samples. The Margalef Richness Index was higher in all eDNA samples compared to morphological samples. Pielou’s Evenness index was lower in the eDNA samples on stations V1 and V3 (about 40% of not assigned ASVs and high share domination of Polar centric Mediophyceae), V7 (about 60% of unassigned ASVs) and V10 (clear domination of Peridiniaceae). Results obtained by the eDNA metabarcoding showed higher Shannon–Wiener Diversity Index values, except for station V10 due to the explicit prevalence of Peridiniaceae.

According to the morphological approach, the phytoplankton community was characterized by dinoflagellate taxa (Miozoa) belonging to functional group Lo, followed by cryptophyte and chlorophyte taxa (mainly Cryptomonas sp. and Tetraselmis cordiformis, respectively) from codon X2 (Figure 4, Table S4). A higher share of centric diatom Pantocsekiella ocellata (Bacillariophyta) assorted into codon C was apparent on central stations (V3–V6) and station V7. Mixotrophic genus Dinobryon (Ochrophyta) from codon E was present across all stations in a small share. According to eDNA metabarcoding approach, members of the family Cryptomonadales belonging to functional group X2 had the highest phytoplankton biomass share on central stations (V3-V6), whilst coda Lo (Peridiniales), C (Polar centric Mediophyceae) and X3 (Chrysophyceae_X) appeared either as dominant or with a high biomass share on other stations (Figure 7 and Figure 8).

4. Discussion

Comparison of identified taxa according to morphological approach indicated a high level of similarity between composite samples and euphotic zone samples in the karst Lake Visovac. A similar number of recorded species was found in both types of samples as well as similar descriptive taxa. Especially, a higher level of similarity in both euphotic zone samples and composite samples was found for the following taxa: C. hirundinella, Cryptomonas sp. and P. ocellata. The high level of similarity was further supported by the ANOSIM analysis, which indicated that euphotic zone samples can be reliably exchanged by composite samples to provide an accurate characterization of phytoplankton communities in the euphotic zone. Moreover, apart from some mismatches in samples from stations V4 and V7, in particular codon Lo, which consequently affected other relative biomass shares, the results on the phytoplankton community from composite samples and euphotic zone samples were mainly congruent. Gaps in shares of phytoplankton coda can be explained due to deviations in biomass when determining species with a large biovolume value, i.e., a difference of just one or two found cells in the observed sample can influence the results [38]. Nevertheless, the composite sampling of phytoplankton can be used as an optimal sampling approach in Lake Visovac during regular monitoring. Lake phytoplankton communities depend on factors such as nutrient availability, light and temperature, as well as hydrological conditions [39], all of which condition a predictably lower phytoplankton biomass in the aphotic zone of Lake Visovac than in the euphotic zone. Moreover, a significantly lower oxygen concentration (<5 mg L⁻¹) was measured on several stations (V3, V5 and V9) in the aphotic zone of Lake Visovac indicating hypoxic conditions. As for the horizontal distribution, the similarity within the phytoplankton community was more evident with respect to the sampling microlocation than to the sampling station, thus resulting in a strong separation between limnetic and littoral zone (the nearshore environment) samples in Lake Visovac. The relative predominance of centric diatoms in both the limnetic and littoral zones may be associated with their capability to resuspend from the bottom due to short retention time, fast water flow, and a relatively deep mixing depth, which can prevent centrics from sinking to the hypolimnion and allow them to dominate in the water column [40]. On the other hand, the flagella-bearing Cryptomonas sp. and Ceratium hirundinella can actively swim in the water column to obtain sufficient amounts of light and nutrients [41,42]. The segregation of littoral samples could be attributed to several biotic factors, such as the impact of macrophytic vegetation on development of phytoplankton, grazing by zooplankton, and interspecific competition, and abiotic factors such as variations in water chemistry and differences in nutrient availability [43,44].

The phytoplankton functional groups are based on ecological sensitivities and tolerances of taxa [5,26], so their main advantage over traditional taxonomic division is the possibility of generalizing the results [45]. Moreover, they can also be successfully integrated in ecological assessments where eDNA metabarcoding datasets are present [8]. The euphotic and aphotic zones were both dominated by a large mixotrophic swimming dinoflagellate C. hirundinella, known to build up higher biomass during the summer stratification period [46,47]. The flagellar motility enables cells to vertically migrate throughout the water column, thus facilitating effective exploitation of nutrients and boosting photosynthesis [26,48,49]. This conclusively allows them to dominate in thermally stratified mesotrophic lakes during the summer period [5]. Besides C. hirundinella, dinoflagellate Parvodinium inconspicuum was also a descriptive representative of Lo, a functional group characteristic of the summer epilimnion of mesotrophic lakes and tolerant of nutrient deficiency [5,26]. Dominance of the Ceratium species has risen significantly with warming caused by climate change [50]. The mixotrophic nutrition strategy is widespread and often dominant in freshwater ecosystems [51,52], enabling phytoplankton to be capable of bacterivory in nutrient-depleted conditions [53,54]. Different mixotrophic species may vary in their ability to regulate the shift between autotrophic and heterotrophic lifestyles, with regards to the changes in the environment, which might cause differences in their temperature response [55]. Changes in the functional role of mixotrophs from primary producers to consumers may cascade through food webs, altering species interactions, as well as the magnitude and direction of the carbon flux [55].

Among other descriptive taxa, larger abundance of cryptophyte Cryptomonas sp. was observed on all stations, whilst centric diatom Pantocsekiella ocellata showed higher abundance on all central stations and lower station V7. Cryptomonas sp. belongs to functional group X2, together with another cryptophyte species Plagioselmis nannoplanctica and chlorophyte Tetraselmis cordiformis [5]. The representatives of codon X2 are acknowledged as meso-eutrophic indicators [5] and exhibit a wide range of tolerance to changes in ecological conditions in Lake Visovac [56]. P. ocellata is sorted into codon C, adaptated to high lake stability [57] and low light availability [5,58], and known to dominate in the mesotrophic ecosystems [59]. P. ocellata was confirmed to be one of the key descriptors of the phytoplankton community in previous investigations on Lake Visovac [8,46,47,56]. Chrysomonads of the mixotrophic genus Dinobryon from codon E were present in lower shares across all stations, and usually denote small, shallow, base-poor lakes or heterotrophic ponds [5,26].

The comparability of traditional light microscopy and eDNA metabarcoding approach is evident in the designation of centric diatoms, dinoflagellates and cryptophytes as descriptive taxa of Lake Visovac. Although dinoflagellate Ceratium hirundinella, diatom Pantocsekiella ocellata and cryptophytes Cryptomonas sp. and Plagioselmis nannoplanctica were identified to the species level by using traditional microscopy, the eDNA metabarcoding method designated the higher ranks of Polar centric Mediophyceae, Cryptomonadales and Peridiniales as the most dominant algal groups. The possible reason that the descriptive taxa were identified only by microsocpy and not by eDNA to the species level is that taxonomic assignments of short amplicon reads to the species level are still problematic because too many species are missing from the reference database [30,31]. For this reason, eDNA analyses of eukaryotic phytoplankton diversity were based at the family level, because metabarcoding at the V9-region of SSU rRNA genes allows correct identification from the genus to the higher taxonomic level, as recognized in previous studies [30,60]. Although there are several more specific primers for detecting diversity of diatoms, such as ribosomal (rbcL) genes [61,62], the small universal hypervariable V9-region of the 18S rRNA gene was chosen in this study because it provides a comprehensive overview of the community and has the ability to capture assemblages of photosynthetic organisms, especially when dealing with phytoplankton consisting of many different algal groups [8,30]. On the other hand, according to the relative abundance, there were a lot of taxonomically unassigned ASVs in the whole molecular dataset. First, the reason for this gap might be the choice of primers, as mentioned above, since they are crucial for species recognition [10]. Second, the method of sampling may also affect the results, as eDNA samples require different volumes of water to be filtered. In this case, water volumes were lower than usual to reduce potential PCR inhibitors from filtering larger volumes of water, but this could potentially limit detection of all taxa present in represented samples [10,63].

When comparing the results of the phytoplankton community defined by functional classification, both approaches proved reliable in detecting functional groups (Lo, C, X2, X3) with similar ecological demands and were congruent with previous studies [8]. Differences between the approaches can be attributed to the ability to distinguish indistinct morphological characteristics, as the existence of cryptic species, pico- and concealed phytoplankton can be difficult to ascertain from morphological analyses [7]. In addition, the most frequent coda were used in the assigning of class and family ranks into functional groups, which could also affect the results. The above mentioned could explain a higher share of functional group X3 to which the family rank Chrysophyceae_X was assigned in the molecular approach.

As for alpha diversity, in most cases, the eDNA metabarcoding results provided higher values than the morphological approach for all indices, except for Pielou’s Index which has shown contrasting results for all samples. This may be related to the occurrence of similar morphological features between microscopically recorded species, which may lead to difficulties in species discrimination [8,64]. In addition, certain small phytoplankton can be easily detected by eDNA metabarcoding, whereas they are usually missed by light microscopy, which may affect species diversity [65,66]. Several previous studies also reported that V9-region has potential for broader recognition spectrum when obtaining results for Shannon diversity [67,68]. On the other hand, there was a discrepancy in the Pielou index, which showed very low values mostly for eDNA samples, which was due to a clear dominance of Polar centric Mediophyceae and Peridiniales [67]. Although the total number of taxonomically assigned ASVs identified by the molecular approach was two times higher than the number of taxa identified by traditional morphological identification, these results should be viewed with caution, especially for those ASVs that could not be taxonomically assigned, as they do not reflect the same percentage or number of unidentified taxa [12]. There is still a problem in translating abundance from sequence data to biological abundance because the rDNA copy number varies among taxa. Therefore, caution should always be used when interpreting the most abundant taxa detected by amplicon sequences, because the sequences of Alveolata (dinoflagellates) show variation in rDNA copy numbers [69]. However, in this study the eDNA results were compared, and to some extent confirmed and verified with the results of the morphological approach.

Non-metric multidimensional scaling (NMDS) analysis demonstrated the corresponding grouping of samples in both morphological and molecular approaches. Similarly comparable results between molecular and morphological approaches using beta diversity were also confirmed in several recent studies [70,71,72]. The results of NMDS analysis indicated a significant dissimilarity between composite samples and aphotic zone samples in both approaches, thus confirming the applicability of eDNA metabarcoding in the routine biomonitoring of Lake Visovac. Significant differences in horizontal and vertical distribution of the phytoplankton in Lake Visovac were found previously by Ciglenečki-Jušić et al. [73].

5. Conclusions

The morphological and eDNA metabarcoding approaches offer comparable results in describing the phytoplankton community and are applicable tools for biodiversity monitoring in Lake Visovac. Moreover, eDNA metabarcoding should be used as a feasible supplement to traditional microscopy in the phytoplankton community assessments, with regards to the drawbacks of each method. It is important to emphasize the essential continual advancement of eDNA metabarcoding in providing more accurate results.

Author Contributions

Formal analysis, investigation, data curation, writing—Original draft preparation, visualization, M.Š.; formal analysis, investigation, writing—Review and editing, visualization, A.K.; formal analysis, investigation, validation, writing—Review and editing, visualization, P.Ž.; conceptualization, investigation, supervision, validation, writing—Review and editing M.G.U. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Acknowledgments

We thank Mirela Šušnjara for contributions in the laboratory and field work and Darija Damić for help with data visualization in QGIS.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Sampling stations on Lake Visovac during August 2018.

Figure 2. Principal Component Analysis (PCA) ordination diagram performed on the environmental variables and sampling stations (V1 to V10, excluding station V8) in Lake Visovac during August 2018. O2—oxygen concentration, O2%—oxygen saturation, T—temperature, SD—Secchi depth, EC—electrical conductivity, NO3−-N—nitrate, PO43−-P—ortophosphate, SiO2—total silica, TIC—total inorganic carbon, DIC—dissolved inorganic carbon, TOC—total organic carbon, DOC—dissolved organic carbon, HCO3−—bicarbonates.

Figure 3. Relative biomass of phytoplankton taxonomic groups (expressed in percentages) in the composite (C), euphotic (EU) and aphotic (NONEU) zone samples from the sampling stations (V1 to V10, excluding station V8) in Lake Visovac during August 2018 according to the morphological approach.

Figure 4. Relative biomass shares of Reynolds’ functional groups (expressed in percentages) in the composite (C), euphotic (EU) and aphotic (NONEU) zone samples from the sampling stations (V1 to V10, excluding station V8) in Lake Visovac during August 2018 according to the morphological approach.

Figure 5. Non-metric multidimensional scaling (NMDS) ordination based on the Bray–Curtis similarity distance in the taxonomic composition of phytoplankton community on sampling stations (V1 to V10, excluding station V8) in Lake Visovac during August 2018 according to morphological approach. 1-EU-euphotic zone samples, 2-NONEU- aphotic zone samples.

Figure 6. Non-metric multidimensional scaling (NMDS) ordination based on Bray–Curtis similarity distance in the taxonomic composition of phytoplankton community of composite samples on sampling stations (V1 to V10, excluding station V8) in Lake Visovac during August 2018 according to the morphological approach.

Figure 7. Relative abundance shares of phytoplankton family ranks (expressed in percentages) in the composite (C) and aphotic zone (NONEU) samples of sampling stations (V1 to V10, excluding station V8) in Lake Visovac during August 2018 according to the molecular (eDNA) approach.

Figure 8. Relative biomass shares of Reynolds’ functional groups (expressed in percentages) in the eDNA composite (C) and aphotic zone (NONEU) samples of sampling stations (V1 to V10, excluding station V8) in Lake Visovac during August 2018.

Figure 9. Non-metric multidimensional scaling (NMDS) ordination based on Bray–Curtis similarity distance in the taxonomic composition of phytoplankton community on sampling stations (V1 to V9, excluding station V8) in Lake Visovac during August 2018 according to the molecular eDNA approach. 1-C- composite samples, 2-NONEU-aphotic zone samples.

Figure 10. Phytoplankton diversity: (a) Total number of phytoplankton ASVs and taxa, (b) Margalef Richness Index (d), (c) Pielou’s Evenness Index (J’), and (d) Shannon–Wiener Diversity Index (H’) for eDNA (blue line) and morphological composite samples (red line) on sampling stations (V1 to V10, excluding station V8) in Lake Visovac during August 2018.

Table 1

Environmental variables on 9 sampling stations (V1 to V10, excluding station V8) for composite samples.

Station	Max Depth (m)	SD *(m)	ZEU *(m)	T *(°C)	O2 *(mg L−1)	O2(%)*	pH	EC *(µS cm−1)
V1	15.2	3.0	7.5	23.8	11.68	138.0	8.00	545
V2	18.0	2.5	6.25	23.3	11.03	129.8	7.99	547
V3	18.0	3.0	7.5	22.0	10.00	115.0	7.84	556
V4	15.5	4.0	10.0	21.5	9.03	102.6	7.79	559
V5	18.0	3.5	8.75	22.1	9.79	113.0	7.94	559
V6	15.0	3.5	8.75	21.7	9.09	103.3	7.77	562
V7	24.0	3.5	8.75	20.6	11.14	124.5	7.75	545
V9	23.0	5.0	12.5	25.9	10.18	125.6	7.84	505
V10	5.0	4.5	5.0	26.2	9.97	123.5	7.98	494

Note: * SD—Secchi depth, Z_EU—euphotic zone depth, T—temperature, O₂—oxygen concentration, O₂ (%)—oxygen saturation, EC—electrical conductivity.

Table 2

One-way analysis of similarity (ANOSIM) between the composite (C), euphotic (EU) and aphotic (NONEU) zone samples in Lake Visovac. Statistically significant values are marked in bold.

Pairwise Tests
Comparison of Samples	R Statistic	Significance Level (p)	Possible Permutations	Actual Permutations	Number >= Observed
C vs. EU	−0.023	0.548	24310	999	547
C vs. NONEU	0.168	0.042	24310	999	41
EU vs. NONEU	0.391	0.006	24310	999	5

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/w15071379/s1. Supplement file: Table S1. Physical and chemical parameters in the euphotic and aphotic zone samples on the sampling stations (V1 to V10, excluding station V8).; Table S2. Chemical parameters measured in the laboratory for composite samples, sampling stations (V1 to V10, excluding station V8).; Table S3. Relative variance explained and factor coordinates of the variables for the first two principal components (PC1 and PC2) of the Principal Component Analysis (PCA).; Figure S1. Shares of all Eukaryota class groups on the sampling stations (V1 to V10, excluding station V8) in Lake Visovac during August 2018 according to analysis of V9-region.; Table S4. Phytoplankton taxa biomass (mg L⁻¹) on the sampling stations (V1 to V10, excluding station V8) in Lake Visovac during August 2018.; Table S5. Assigned taxonomic amplicon sequencing variants for phytoplankton community on the sampling stations (V1 to V10, excluding station V8) in Lake Visovac during August 2018.

References

1. Myers, N.; Mittermeier, R.A.; Mittermeier, C.G.; da Fonseca, G.A.B.; Kent, J. Biodiversity Hotspots for Conservation Priorities. Nature; 2000; 403, pp. 853-858. [DOI: https://dx.doi.org/10.1038/35002501] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10706275]

2. Mikac, I.; Fiket, Ž.; Terzić, S.; Barešić, J.; Mikac, N.; Ahel, M. Chemical Indicators of Anthropogenic Impacts in Sediments of the Pristine Karst Lakes. Chemosphere; 2011; 84, pp. 1140-1149. [DOI: https://dx.doi.org/10.1016/j.chemosphere.2011.04.027] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21531436]

3. Gaedke, U. Trophic Dynamics in Aquatic Ecosystems. Encyclopedia of Inland Waters; Elsevier: New York, NY, USA, 2009; pp. 499-504. ISBN 978-0-12-370626-3

4. European Commision Directive of the European Parliament and of the Council. 2000/60/EC. Establishing a Framework for Community Action in the Field of Water Policy (Water Framework Directive); The Publications Office of the European Union (L’Office des publications de l’Union européenne: Luxembourg, 2000.

5. Reynolds, C.S.; Huszar, V.; Kruk, C.; Naselli-Flores, L.; Melo, S. Towards a Functional Classification of the Freshwater Phytoplankton. J. Plankton Res.; 2002; 24, pp. 417-428. [DOI: https://dx.doi.org/10.1093/plankt/24.5.417]

6. Bellinger, E.G.; Sigee, D.C. Freshwater Algae: Identification, Enumeration and Use as Bioindicators; 1st ed. Wiley-Blackwell: New York, NY, USA, 2015; ISBN 978-1-118-91716-9

7. Huo, S.; Li, X.; Xi, B.; Zhang, H.; Ma, C.; He, Z. Combining Morphological and Metabarcoding Approaches Reveals the Freshwater Eukaryotic Phytoplankton Community. Environ. Sci. Eur.; 2020; 32, 37. [DOI: https://dx.doi.org/10.1186/s12302-020-00321-w]

8. Hanžek, N.; Gligora Udovič, M.; Kajan, K.; Borics, G.; Várbíró, G.; Stoeck, T.; Žutinić, P.; Orlić, S.; Stanković, I. Assessing Ecological Status in Karstic Lakes through the Integration of Phytoplankton Functional Groups, Morphological Approach and Environmental DNA Metabarcoding. Ecol. Indic.; 2021; 131, 108166. [DOI: https://dx.doi.org/10.1016/j.ecolind.2021.108166]

9. Borics, G.; Abonyi, A.; Salmaso, N.; Ptacnik, R. Freshwater Phytoplankton Diversity: Models, Drivers and Implications for Ecosystem Properties. Hydrobiologia; 2021; 848, pp. 53-75. [DOI: https://dx.doi.org/10.1007/s10750-020-04332-9]

10. Pawlowski, J.; Apothéloz-Perret-Gentil, L.; Mächler, E.; Altermatt, F. Environmental DNA Applications in Biomonitoring and Bioassessment of Aquatic Ecosystems; Environmental Studies. no 2010 Federal Office for the Environment: Bern, Switzerland, 2020; 71. [DOI: https://dx.doi.org/10.5167/uzh-187800]

11. Hering, D.; Borja, A.; Jones, J.I.; Pont, D.; Boets, P.; Bouchez, A.; Bruce, K.; Drakare, S.; Hänfling, B.; Kahlert, M. et al. Implementation Options for DNA-Based Identification into Ecological Status Assessment under the European Water Framework Directive. Water Res.; 2018; 138, pp. 192-205. [DOI: https://dx.doi.org/10.1016/j.watres.2018.03.003]

12. Kulaš, A.; Gligora Udovič, M.; Tapolczai, K.; Žutinić, P.; Orlić, S.; Levkov, Z. Diatom EDNA Metabarcoding and Morphological Methods for Bioassessment of Karstic River. Sci. Total Environ.; 2022; 829, 154536. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2022.154536]

13. Bonacci, O.; Andrić, I.; Roje-Bonacci, T. Hydrological Analysis of Skradinski Buk Tufa Waterfall (Krka River, Dinaric Karst, Croatia). Env. Earth Sci.; 2017; 76, 669. [DOI: https://dx.doi.org/10.1007/s12665-017-7023-9]

14. Official Gazette of the Republic of Croatia No. 96/19 Regulation on Water Quality Standards. Zagreb, Croatia, 2019. Available online: https://leap.unep.org/countries/hr/national-legislation/regulation-water-quality-standard (accessed on 6 March 2023).

15. Poikane, S. Water Framework Directive Intercalibration Technical Report; Part 2: Lakes Office for Official Publications of the European Communities (OPOCE): Luxembourg, 2009.

16. APHA Standard Methods for the Examination of Water and Wastewater. American Public Health Association: Washington, DC, USA, 2017. Available online: https://www.scirp.org/(S(351jmbntvnsjt1aadkozje))/reference/referencespapers.aspx?referenceid=2459667 (accessed on 6 March 2023).

17. Utermöhl, H. Zur Vervollkommnung Der Quantitativen Phytoplankton-Methodik. SIL Commun.; 1958; 9, pp. 1-38. [DOI: https://dx.doi.org/10.1080/05384680.1958.11904091]

18. Krammer, K.; Lange-Bertalot, H. Bacillariophyceae 3. Teil: Centrales, Fragilariaceae, Eunotiaceae. Süßwasserflora von Mitteleuropa Band 2/3; Ettl, H.; Gerloff, J.; Heynig, H.; Mollenhauer, D. Spektrum Akademischer Verlag: Heidelberg, Germany, 2000.

19. Komárek, J.; Anagnostidis, K. Cyanoprokaryota 1. Teil: Chroococcales. Süsswasserflora von Mitteleuropa 19/1; Ettl, H.; Gärtner, G.; Heynig, H.; Mollenhauer, D. Gustav Fischer: Jena, Germany, Stuttgart, Germany, Lübeck, Germany, Ulm, Germany, 1998.

20. Komárek, J.; Anagnostidis, K. Cyanoprokaryota 2. Teil/ 2nd Part: Oscillatoriales. Süsswasserflora von Mitteleuropa 19/2; Büdel, B.; Krienitz, L.; Gärtner, G.; Schagerl, M. Elsevier/Spektrum: Heidelberg, Germany, 2005.

21. Krammer, K. Bacillariophyceae; Spektrum Akademischer: Heidelberg, Germany, 2000.

22. Kristiansen, J.; Preisig, H.R. Süßwasserflora von Mitteleuropa 1. Chrysophyte and Haptophyte Algae: Part 2: Synurophyceae; 2nd ed. Spektrum Akademisher Verlag: Berlin, Germany, 2007.

23. Hofmann, G.; Werum, M.; Lange-Bertalot, H. Diatomeen Im Süßwasser—Benthos von Mitteleuropa. Bestimmungsflora Kieselalgen Für Die Ökologische Praxis; Über 700 Der Häufigsten Arten Und Ihre Ökologie Gantner, A.R.G. Koeltz Scientific Books: Ruggell, Liechtenstein, 2011.

24. John, D.M.; Whitton, B.A.; Brook, A.J. The Freshwater Algal Flora of the British Isles: An Identification Guide to Freshwater and Terrestrial Algae; Cambridge University Press: Cambridge, UK, 2011.

25. Guiry, M.D.; Guiry, G.M. AlgaeBase. World-Wide Electronic Publication; National University of Ireland: Galway, Ireland, 2022; Available online: https://www.algaebase.org (accessed on 3 August 2022).

26. Padisák, J.; Crossetti, L.O.; Naselli-Flores, L. Use and Misuse in the Application of the Phytoplankton Functional Classification: A Critical Review with Updates. Hydrobiologia; 2009; 621, pp. 1-19. [DOI: https://dx.doi.org/10.1007/s10750-008-9645-0]

27. HRN EN ISO 5667-3:2018 Water Quality—Sampling—Part 3: Preservation and Handling of Water Samples. Available online: https://repozitorij.hzn.hr/norm/HRN+EN+ISO+5667-3%3A2018 (accessed on 16 January 2019).

28. HRN EN 15204:2008; Water Quality—Guidance Standard on the Enumeration of Phytoplankton Using Inverted Microscopy (Utermöhl Technique). Available online: https://repozitorij.hzn.hr/norm/HRN+EN+15204%3A2008 (accessed on 11 June 2020).

29. HRN EN 16695:2015; Water Quality—Guidance on the Estimation of Phytoplankton Biovolume. Available online: https://repozitorij.hzn.hr/norm/HRN+EN+16695%3A2015 (accessed on 11 June 2020).

30. Stoeck, T.; Bass, D.; Nebel, M.; Christen, R.; Jones, M.D.M.; Breiner, H.-W.; Richards, T.A. Multiple Marker Parallel Tag Environmental DNA Sequencing Reveals a Highly Complex Eukaryotic Community in Marine Anoxic Water. Mol. Ecol.; 2010; 19, pp. 21-31. [DOI: https://dx.doi.org/10.1111/j.1365-294X.2009.04480.x]

31. Amaral-Zettler, L.A.; McCliment, E.A.; Ducklow, H.W.; Huse, S.M. A Method for Studying Protistan Diversity Using Massively Parallel Sequencing of V9 Hypervariable Regions of Small-Subunit Ribosomal RNA Genes. PLoS ONE; 2009; 4, e6372. [DOI: https://dx.doi.org/10.1371/annotation/50c43133-0df5-4b8b-8975-8cc37d4f2f26]

32. Stoeck, T.; Kochems, R.; Forster, D.; Lejzerowicz, F.; Pawlowski, J. Metabarcoding of Benthic Ciliate Communities Shows High Potential for Environmental Monitoring in Salmon Aquaculture. Ecol. Indic.; 2018; 85, pp. 153-164. [DOI: https://dx.doi.org/10.1016/j.ecolind.2017.10.041]

33. Martin, M. Cutadapt Removes Adapter Sequences from High-Throughput Sequencing Reads. EMBnet. J.; 2011; 17, 10. [DOI: https://dx.doi.org/10.14806/ej.17.1.200]

34. Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data; Babraham Bioinformatics Babraham Institute: Cambridge, UK, 2010; Available online: http://www.bioinformatics.babraham.ac.uk/projects/fastqc (accessed on 3 May 2019).

35. Callahan, B.J.; McMurdie, P.J.; Rosen, M.J.; Han, A.W.; Johnson, A.J.A.; Holmes, S.P. DADA2: High-Resolution Sample Inference from Illumina Amplicon Data. Nat. Methods; 2016; 13, pp. 581-583. [DOI: https://dx.doi.org/10.1038/nmeth.3869]

36. Guillou, L.; Bachar, D.; Audic, S.; Bass, D.; Berney, C.; Bittner, L.; Boutte, C.; Burgaud, G.; de Vargas, C.; Decelle, J. et al. The Protist Ribosomal Reference Database (PR2): A Catalog of Unicellular Eukaryote Small Sub-Unit RRNA Sequences with Curated Taxonomy. Nucleic Acids Res.; 2012; 41, pp. D597-D604. [DOI: https://dx.doi.org/10.1093/nar/gks1160]

37. QGIS A Free and Open Source Geographic Information System. Open Source. Geospatial Foundation Project. 2009; Available online: https://qgis.org/en/site/ (accessed on 21 March 2023).

38. Zarauz, L.; Irigoien, X.; Fernandes, J.A. Changes in Plankton Size Structure and Composition, during the Generation of a Phytoplankton Bloom, in the Central Cantabrian Sea. J. Plankton Res.; 2008; 31, pp. 193-207. [DOI: https://dx.doi.org/10.1093/plankt/fbn107]

39. Ptacnik, R.; Lepistö, L.; Willén, E.; Brettum, P.; Andersen, T.; Rekolainen, S.; Lyche Solheim, A.; Carvalho, L. Quantitative Responses of Lake Phytoplankton to Eutrophication in Northern Europe. Aquat. Ecol.; 2008; 42, pp. 227-236. [DOI: https://dx.doi.org/10.1007/s10452-008-9181-z]

40. Žutinić, P.; Gligora Udovič, M.; Kralj Borojević, K.; Plenković-Moraj, A.; Padisák, J. Morpho-Functional Classifications of Phytoplankton Assemblages of Two Deep Karstic Lakes. Hydrobiologia; 2014; 740, pp. 147-166. [DOI: https://dx.doi.org/10.1007/s10750-014-2115-y]

41. Clegg, M.R.; Maberly, S.C.; Jones, R.I. Behavioral Response as a Predictor of Seasonal Depth Distribution and Vertical Niche Separation in Freshwater Phytoplanktonic Flagellates. Limnol. Oceanogr.; 2007; 52, pp. 441-455. [DOI: https://dx.doi.org/10.4319/lo.2007.52.1.0441]

42. Jamal, A.; Yusoff, F.M.; Banerjee, S.; Shariff, M. Littoral and Limnetic Phytoplankton Distribution and Biodiversity in a Tropical Man-Made Lake, Malaysia. Adv. Stud. Biol.; 2014; 6, pp. 149-168. [DOI: https://dx.doi.org/10.12988/asb.2014.4631]

43. Lemly, A.D.; Dimmick, J.F. Phytoplankton Communities in the Littoral Zone of Lakes: Observations on Structure and Dynamics in Oligotrophic and Eutrophic Systems. Oecologia; 1982; 54, pp. 359-369. [DOI: https://dx.doi.org/10.1007/BF00380005] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28309960]

44. Mulderij, G.; Van Nes, E.H.; Van Donk, E. Macrophyte–Phytoplankton Interactions: The Relative Importance of Allelopathy versus Other Factors. Ecol. Model.; 2007; 204, pp. 85-92. [DOI: https://dx.doi.org/10.1016/j.ecolmodel.2006.12.020]

45. Blaum, N.; Mosner, E.; Schwager, M.; Jeltsch, F. How Functional Is Functional? Ecological Groupings in Terrestrial Animal Ecology: Towards an Animal Functional Type Approach. Biodivers. Conserv.; 2011; 20, pp. 2333-2345. [DOI: https://dx.doi.org/10.1007/s10531-011-9995-1]

46. Gligora Udovič, M.; Kralj Borojević, K.; Žutinić, P.; Šipoš, L.; Plenković-Moraj, A. Net-Phytoplankton Species Dominance in a Travertine Riverine Lake Visovac, NP Krka. Nat. Croat.; 2011; 20, pp. 411-424.

47. Gligora Udovič, M.; Žutinić, P.; Kralj Borojević, K.; Plenković-Moraj, A. Co-Occurrence of Functional Groups in Phytoplankton Assemblages Dominated by Diatoms, Chrysophytes and Dinoflagellates. Fundam. Appl. Limnol.; 2015; 187, pp. 101-111. [DOI: https://dx.doi.org/10.1127/fal/2015/0759]

48. Sommer, U. Growth and Survival Strategies of Plankton Diatoms. Growth and Reproductive Strategies of Freshwater Phytoplankton; Sandgren, C.D. Cambridge University Press: Cambridge, UK, 1988; pp. 227-260. ISBN 0-521-32722-9

49. Grigorszky, I.; Padisák, J.; Borics, G.; Schitchen, C.; Borbély, G. Deep Chlorophyll Maximum by Ceratium hirundinella (O. F. Müller) Bergh in a Shallow Oxbow in Hungary. Hydrobiologia; 2003; 506, pp. 209-212. [DOI: https://dx.doi.org/10.1023/B:HYDR.0000008632.57769.19]

50. Olrik, K.; Cronberg, G.; Annadotter, H. Lake Phytoplankton Responses to Global Climate Changes. Climatic Change and Global Warming of Inland Waters: Impacts and Mitigation for Ecosystems and Societies; Goldman, C.R.; Kumagai, M.; Robarts, R.D. Wiley: Chichester, UK, 2012; pp. 173-199.

51. Stoecker, D.K.; Johnson, M.D.; deVargas, C.; Not, F. Acquired Phototrophy in Aquatic Protists. Aquat. Microb. Ecol.; 2009; 57, pp. 279-310. [DOI: https://dx.doi.org/10.3354/ame01340]

52. Sanders, R.W. Alternative Nutritional Strategies in Protists: Symposium Introduction and a Review of Freshwater Protists That Combine Photosynthesis and Heterotrophy1. J. Eukaryot. Microbiol.; 2011; 58, pp. 181-184. [DOI: https://dx.doi.org/10.1111/j.1550-7408.2011.00543.x]

53. Kamjunke, N.; Henrichs, T.; Gaedke, U. Phosphorus Gain by Bacterivory Promotes the Mixotrophic Flagellate Dinobryon Spp. during Re-Oligotrophication. J. Plankton Res.; 2007; 29, pp. 39-46. [DOI: https://dx.doi.org/10.1093/plankt/fbl054]

54. Hansen, P.J. The Role of Photosynthesis and Food Uptake for the Growth of Marine Mixotrophic Dinoflagellates1. J. Eukaryot. Microbiol.; 2011; 58, pp. 203-214. [DOI: https://dx.doi.org/10.1111/j.1550-7408.2011.00537.x]

55. Wilken, S.; Huisman, J.; Naus-Wiezer, S.; Van Donk, E. Mixotrophic Organisms Become More Heterotrophic with Rising Temperature. Ecol. Lett.; 2013; 16, pp. 225-233. [DOI: https://dx.doi.org/10.1111/ele.12033]

56. Šimunović, M.; Kulaš, A.; Žutinić, P.; Goreta, G.; Gligora Udovič, M. Phytoplankton Metrics for Trophic and Ecological Status Assessment of a Natural Karstic Lake. Acta Bot. Croat.; 2022; 81, pp. 185-196. [DOI: https://dx.doi.org/10.37427/botcro-2022-016]

57. Beamud, S.G.; León, J.G.; Kruk, C.; Pedrozo, F.; Diaz, M. Using Trait-Based Approaches to Study Phytoplankton Seasonal Succession in a Subtropical Reservoir in Arid Central Western Argentina. Environ. Monit. Assess.; 2015; 187, 271. [DOI: https://dx.doi.org/10.1007/s10661-015-4519-1]

58. Reynolds, C.S. Vegetation Processes in the Pelagic: A Model for Ecosystem Theory. Xxvii, 371 p. Oldendorf/Luhe, Germany: Ecology Institute, 1997. (Excellence in Ecology No. 9). Price DM 68.00. J. Mar. Biol. Ass.; 1997; 77, 919. [DOI: https://dx.doi.org/10.1017/S0025315400036432]

59. Hu, R.; Han, B.; Naselli-Flores, L. Comparing Biological Classifications of Freshwater Phytoplankton: A Case Study from South China. Hydrobiologia; 2013; 701, pp. 219-233. [DOI: https://dx.doi.org/10.1007/s10750-012-1277-8]

60. Choi, J.; Park, J.S. Comparative Analyses of the V4 and V9 Regions of 18S RDNA for the Extant Eukaryotic Community Using the Illumina Platform. Sci. Rep.; 2020; 10, 6519. [DOI: https://dx.doi.org/10.1038/s41598-020-63561-z]

61. Rimet, F.; Vasselon, V.; A.-Keszte, B.; Bouchez, A. Do We Similarly Assess Diversity with Microscopy and High-Throughput Sequencing? Case of Microalgae in Lakes. Org. Divers. Evol.; 2018; 18, pp. 51-62. [DOI: https://dx.doi.org/10.1007/s13127-018-0359-5]

62. Vasselon, V.; Domaizon, I.; Rimet, F.; Kahlert, M.; Bouchez, A. Application of High-Throughput Sequencing (HTS) Metabarcoding to Diatom Biomonitoring: Do DNA Extraction Methods Matter?. Freshw. Sci.; 2017; 36, pp. 162-177. [DOI: https://dx.doi.org/10.1086/690649]

63. Bruce, K.; Blackman, R.; Bourlat, S.J.; Hellström, A.M.; Bakker, J.; Bista, I.; Bohmann, K.; Bouchez, A.; Brys, R.; Clark, K. et al. A Practical Guide to DNA-Based Methods for Biodiversity Assessment; Pensoft Publishers: Sofia, Bulgaria, 2021; ISBN 978-619-248-053-0

64. Wilmotte, A.; Dail Laughinghouse IV, H.; Capelli, C.; Rippka, R.; Salmaso, N. Taxonomic Identification of Cyanobacteria by a Polyphasic Approach. Molecular Tools for the Detection and Quantification of Toxigenic Cyanobacteria; John Wiley & Sons, Ltd.: Chichester, UK, 2017; pp. 79-134.

65. Not, F.; Zapata, M.; Pazos, Y.; Campaña, E.; Doval, M.; Rodríguez, F. Size-Fractionated Phytoplankton Diversity in the NW Iberian Coast: A Combination of Microscopic, Pigment and Molecular Analyses. Aquat. Microb. Ecol.; 2007; 49, pp. 255-265. [DOI: https://dx.doi.org/10.3354/ame01144]

66. Xiao, X.; Sogge, H.; Lagesen, K.; Tooming-Klunderud, A.; Jakobsen, K.S.; Rohrlack, T. Use of High Throughput Sequencing and Light Microscopy Show Contrasting Results in a Study of Phytoplankton Occurrence in a Freshwater Environment. PLoS ONE; 2014; 9, e106510. [DOI: https://dx.doi.org/10.1371/journal.pone.0106510] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25171164]

67. Maritz, J.M.; Rogers, K.H.; Rock, T.M.; Liu, N.; Joseph, S.; Land, K.M.; Carlton, J.M. An 18S RRNA Workflow for Characterizing Protists in Sewage, with a Focus on Zoonotic Trichomonads. Microb. Ecol.; 2017; 74, pp. 923-936. [DOI: https://dx.doi.org/10.1007/s00248-017-0996-9] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28540488]

68. Tragin, M.; Zingone, A.; Vaulot, D. Comparison of Coastal Phytoplankton Composition Estimated from the V4 and V9 Regions of the 18S RRNA Gene with a Focus on Photosynthetic Groups and Especially Chlorophyta. Environ. Microbiol.; 2018; 20, pp. 506-520. [DOI: https://dx.doi.org/10.1111/1462-2920.13952]

69. Medinger, R.; Nolte, V.; Pandey, R.V.; Jost, S.; Ottenwalder, B.; Schlotterer, C.; Boenigk, J. Diversity in a Hidden World: Potential and Limitation of next-Generation Sequencing for Surveys of Molecular Diversity of Eukaryotic Microorganisms. Mol. Ecol.; 2010; 19, pp. 32-40. [DOI: https://dx.doi.org/10.1111/j.1365-294X.2009.04478.x]

70. Lanzén, A.; Lekang, K.; Jonassen, I.; Thompson, E.M.; Troedsson, C. High-Throughput Metabarcoding of Eukaryotic Diversity for Environmental Monitoring of Offshore Oil-Drilling Activities. Mol. Ecol.; 2016; 25, pp. 4392-4406. [DOI: https://dx.doi.org/10.1111/mec.13761]

71. Maitland, V.C.; Robinson, C.V.; Porter, T.M.; Hajibabaei, M. Freshwater Diatom Biomonitoring through Benthic Kick-Net Metabarcoding. PLoS ONE; 2020; 15, e0242143. [DOI: https://dx.doi.org/10.1371/journal.pone.0242143]

72. Kahlert, M.; Bailet, B.; Chonova, T.; Karjalainen, S.M.; Schneider, S.C.; Tapolczai, K. Same Same, but Different: The Response of Diatoms to Environmental Gradients in Fennoscandian Streams and Lakes—Barcodes, Traits and Microscope Data Compared. Ecol. Indic.; 2021; 130, 108088. [DOI: https://dx.doi.org/10.1016/j.ecolind.2021.108088]

73. Ciglenečki-Jušić, I.; Ahel, M.; Mikac, N.; Omanović, D.; Vdović, N. Investigation of Natural Characteristics and Assesment of Antropogenic Influences on Water Quality of Visovac Lake; Ruđer Bošković Institute: Zagreb, Croatia, 2013; Available online: https://www.bib.irb.hr/662477 (accessed on 21 June 2020).