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

: 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, dinoﬂagellates 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 conﬁrmed 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.


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 Water 2023, 15, 1379 2 of 20 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.

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 6 m 3 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].

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.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 (ZSD) was determined with a Secchi disc and used for the calculation of euphotic zone depth (ZEU) 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 3− -P), total silica (SiO 2 ), total inorganic carbon (TIC), dissolved inorganic carbon (DIC), total organic carbon (TOC), dissolved organic carbon (DOC) and bicarbonates (HCO 3 − ) 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 appro-priate 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-29].

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 Pow-erWater 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).

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.11pipeline 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.

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].

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 2 concentration was recorded on station V4 (9.03 mg L −1 ) and the highest on station V1 (11.68 mg L −1 ).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 −1 on station V10 to 562 µS cm −1 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 2 was recorded on station V10 (10.13 mg L −1 ), whilst the highest concentration and saturation of O 2 in the euphotic zone were measured on station V7 (12.28 mg L −1 and 147.5%, respectively).The lowest saturation of O 2 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 −1 on station V10 to a maximum of 562 µS cm −1 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 2 were recorded on station V9 (2.91 mg L −1 and 27.4%, respectively) and the highest on station V1 (9.60 mg L −1 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 −1 to 587 µS cm −1 (stations V9 and V1, respectively).
Chemical parameters of water are presented in Table S2.The concentration of NO 3 − -N was very low on most stations (<0.1 mg L −1 ), with slightly higher values on stations V4 (0.7 mg L −1 ) and V5 (0.3 mg L −1 ).Very low concentrations of NO recorded on all stations (<0.001 mg L −1 and <0.01 mg L −1 , respectively).The concentration of PO 4 3− -P varied from the lowest measured on six stations in total (<0.01 mg L −1 ) to the highest measured on station V9 (0.6 mg L −1 ).The values of SiO 2 ranged from the lowest on station V9 (1.5 mg L −1 ) to the highest on station V6 (4.5 mg L −1 ).The highest concentration of TN was measured on stations V4 and V5 (1 mg L −1 ), while it was very low on the other eight stations (<1 mg L −1 ).The lowest concentrations of inorganic carbon compounds were detected on station V3 (TIC and DIC of 11.48 mg L −1 and 10.32 mg L −1 , respectively), while the highest was present on station V7 (TIC and DIC of 13.98 mg L −1 and 13.83 mg L −1 , respectively).TOC was in the range between 1.55 mg L −1 (station V3) and 2.39 mg L −1 (station V7).The values of DOC ranged from 0.80 mg L −1 (V2) to 2.37 mg L −1 (V1), whilst HCO 3 − ranged from 151 mg L −1 (V4) to 212 mg L −1 (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 2 − -N and NH 4 + -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 2 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.
Water 2023, 15, x FOR PEER REVIEW 6 o of PO4 3− -P varied from the lowest measured on six stations in total (<0.01 mg L −1 ) to highest measured on station V9 (0.6 mg L −1 ).The values of SiO2 ranged from the low on station V9 (1.5 mg L −1 ) to the highest on station V6 (4.5 mg L −1 ).The highest concen tion of TN was measured on stations V4 and V5 (1 mg L −1 ), while it was very low on other eight stations (<1 mg L −1 ).The lowest concentrations of inorganic carbon compou were detected on station V3 (TIC and DIC of 11.48 mg L −1 and 10.32 mg L −1 , respective while the highest was present on station V7 (TIC and DIC of 13.98 mg L −1 and 13.83 L −1 , respectively).TOC was in the range between 1.55 mg L −1 (station V3) and 2.39 mg (station V7).The values of DOC ranged from 0.80 mg L −1 (V2) to 2.37 mg L −1 (V1), wh HCO3 -ranged from 151 mg L −1 (V4) to 212 mg L −1 (V7).
Principal component analysis (PCA) performed for the 14 environmental variab explained 59% of the total variance on the first two PC axes (Table S3).NO2 --N and NH N were excluded from PCA because their concentrations did not vary across the statio The most important parameters for PCA axis 1 were temperature, DIC and electrical c ductivity (intra-set correlations: 0.433, −0.361 and −0.338, respectively).Regarding axi O2 concentration, nitrate and TOC were the variables with the most weight for ordinat (intra-set correlations: −0.464, 0.369 and −0.368, respectively).PCA arranged samples (F ure 2) into four groups: the first group consisted of samples from the central stations ( V4, V5, V6) and the upper station sample V2, whilst the second group included samp from the lower part of the Lake (V9 and V10).Sample V1 from the uppermost part a sample V7 taken from the lower part of the Lake were singled out.
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).
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).The descriptive Reynolds' functional groups on upper stations V1 and V2 w and X2 in both composite and euphotic zone samples (Figure 4).The composite a photic zone samples on station V3 were characterized by functional group X2 wi highest biomass share (50% and 37%, respectively), followed by associations C a (25% and 18%; 26% and 33%, respectively).Functional group Lo dominated the a blage in the composite sample on station V4, with 68% of the total phytoplankton bio whereas associations C and X2 were both contributing with 12% to the total phyto ton biomass.Meanwhile, the functional group X2 dominated the assemblage (51% coda C and Lo being subdominant (30% and 15%, respectively) in the euphotic zon ples on station V4.The functional group X2 was dominant, whilst codon C appea subdominant in both composite and euphotic zone samples on station V5.Group Lo became the most dominant in both samples on station V6, followed by coda X2 a Codon C prevailed in the composite sample on station V7, together with associat (43% and 39% of the total phytoplankton biomass, respectively), whilst in the eu zone samples, codon Lo gained dominance (50%).Codon Lo dominated the assem with 80% of the total phytoplankton biomass in the composite sample on station was also the case in the euphotic zone samples, although with a lower share (49% descriptive functional group in the composite sample on station V10 was X2 (57% total phytoplankton biomass), whereas in the euphotic zones sample, it co-dom with codon Lo.Species belonging to codon E (genus Dinobryon) were present in low biomass share on every station, except for the somewhat higher shares noted composite sample on station V10 and the euphotic zone sample on station V9 (11 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).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.The NMDS analysis of phytoplankton taxonomic composition according to phological approach (Figure 5) indicated segregation of almost all discrete-dept into two separate groups, the first one comprising the euphotic zone and the s including the aphotic zone.The samples from station V10 were all counted in photic zone.Considering the horizontal distribution, NMDS analysis of the composit showed grouping of samples from limnetic and littoral stations (Figure 6).Cry sp. and P. ocellata dominated on central limnetic stations V3 and V5.Stations littoral zone were clustered into two distinct groups as follows: the first group, d  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.

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).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.

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.
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).

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

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).
Water 2023, 15, x FOR PEER REVIEW 14 of 20 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 (Figures 7 and 8).

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 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 (Figures 7 and 8).

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 −1 ) 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].

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.

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

Figure 3 .
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 .
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 4 .
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 4 .
Figure 4. Relative biomass shares of Reynolds' functional groups (expressed in percenta composite (C), euphotic (EU) and aphotic (NONEU) zone samples from the sampling s to V10, excluding station V8) in Lake Visovac during August 2018 according to the mo approach.

Figure 5 .
Figure 5. Non-metric multidimensional scaling (NMDS) ordination based on the Bray-C larity distance in the taxonomic composition of phytoplankton community on sampling s to V10, excluding station V8) in Lake Visovac during August 2018 according to morpho proach.1-EU-euphotic zone samples, 2-NONEU-aphotic zone samples.

Figure 5 .
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 .
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 6 .
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 .
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 7 .
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 .
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 8 .
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.

20 Figure 9 .
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 9 .
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.

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.

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.