Diversity of the Summer Phytoplankton of 43 Waterbodies in Bulgaria and Its Potential for Water Quality Assessment

: The general awareness of the threats on biodiversity and water quality raised the number of studies that use phytoplankton in assessment procedures. Since most metrics require obtaining mean values, this paper presents data that may help speed up ﬁeld work and ﬁnd indicators for a rapid water quality assessment based on single samplings, allowing simultaneous work on many sites. The phytoplankton from 43 Bulgarian waterbodies collected during three summer campaigns (2018, 2019, 2021) at sites selected after drone observations was studied by conventional light microscopy (LM) and an HPLC analysis of marker pigments. Our results allowed us to recommend drones and the HPLC application as reliable methods in rapid water quality assessments. In total, 787 algae from seven phyla (53 alien, new for Bulgaria) were identiﬁed. Chlorophyta was the taxonomically richest group, but Cyanoprokaryota dominated the biomass in most sites. New PCR data obtained on anatoxin and microcystin producers conﬁrmed the genetic diversity of Cuspidothrix and Microcystis and provided three new species for the country’s toxic species, ﬁrst identiﬁed by LM. A statistical analysis revealed signiﬁcant correlations of certain algal phyla and classes with different environmental variables, and their species are considered promising for future search of bioindicators. This is especially valid for the class Eustigmatophyceae, which, as of yet, has been almost neglected in water assessment procedures and indices.


Introduction
Since mankind has existed, water has been one of the most important and precious resources of our planet. It is commonly recognized that the lifestyle, agriculture and industry of the modern society, experienced during the last century, led to climate changes and nutrient enrichment of waters, which, in turn, caused a considerable impact on the aquatic habitats. These changes provoked the interest of the scientific community with an increasing intensity of studies on all characteristics of water regarding its use, united by the term "water quality", and its assessment and management [1][2][3][4][5]. Since the end of the 19th century, they have been related with the inhabitants of aquatic biotopes and their potential role in bioindication (for historical details see [6]). Although today, different approaches serve to assess water quality, the use of primary producers with a short life cycle, such as phytoplankters, has a long and worldwide-known tradition [6][7][8][9][10][11]. The methodological tools applied involve certain indicator species or different functional groups, but also the total composition and indices based on diversity, sensitive to the number of species or to their quantitative role [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. Over the years, the research has also focused on so-called algal Table 1. Sampling sites in Bulgarian waterbodies and their environmental parameters during summer sampling campaigns in 2018, 2019 and 2021. Legend: WBN-name of the waterbody, IBWidentification number in the Inventory of Bulgarian Wetlands [37], Abbr-abbreviation of the name, Type-type of waterbody: M (small reservoir/"microreservoir", <100 ha), R (large reservoir, >100 ha) and L (natural lake), Alt-altitude above the sea level (m), WT-water temperature ( • C), CNconductivity (S m −1 ), TDS-total dissolved solids (µg L −1 ), DO-oxygen concentration (mg L −1 ), TP-total phosphorus (µg L −1 ), TN-total nitrogen (mg L −1 ). Waterbodies are presented according to their geographical location in counterclockwise order, starting from South-Western Bulgaria. Asterisks indicate the waterbodies which were sampled for the first time. WBN    The sampling was preceded by a drone sent to observe in real time the whole water area of each waterbody ( Figure 1) and to identify the sites with algal blooms [38][39][40][41][42][43][44][45]. In cases of visible water homogeneity, the sites from our previous studies were repeated, or new sites were selected in cases of waterbodies sampled for the first time. Two types of drones (each supplied by a photo camera) were used: DJI Mavic Pro, Model: M1P GL200A (SZ DJI Technology Co., LTD, Shenzhen, China) in 2018 and DJI Mavic 2 Enterprise Dual Pro (DJI Technology Co, LTD, Shenzhen, China) in 2019, 2021 because the latter had the ability to measure the surface water temperature [38][39][40][41][42][43][44][45].
The sampling was conducted from inflatable boats and by motorboats in the large reservoirs. Aquameter AM-200 and Aquaprobe AP-2000 from Aquaread's water-monitoring instruments, 2012 Aquaread Ltd., were applied for in situ measurement of the coordinates and the altitude of each site, as well as the water temperature, pH, water hardness (expressed by total dissolved solids), oxygen concentration, chlorophyll a and conductivity (Table 1). Total nitrogen (TN) and total phosphorus (TP) were measured ex situ with an Aqualytic AL410 photometer from AQUALYTIC ® , Dortmund, Germany- Table 1. The sampling was preceded by a drone sent to observe in real time the whole water area of each waterbody ( Figure 1) and to identify the sites with algal blooms [38][39][40][41][42][43][44][45]. In cases of visible water homogeneity, the sites from our previous studies were repeated, or new sites were selected in cases of waterbodies sampled for the first time. Two types of drones (each supplied by a photo camera) were used: DJI Mavic Pro, Model: M1P GL200A (SZ DJI Technology Co., LTD, Shenzhen, China) in 2018 and DJI Mavic 2 Enterprise Dual Pro (DJI Technology Co, LTD, Shenzhen, China) in 2019, 2021 because the latter had the ability to measure the surface water temperature [38][39][40][41][42][43][44][45].  [50,51]) with locations of the studied waterbodies and indication of their type. The waterbodies are represented by numbers that follow those in Table 1.
The sampling was conducted from inflatable boats and by motorboats in the large reservoirs. Aquameter AM-200 and Aquaprobe AP-2000 from Aquaread's water-monitoring instruments, 2012 Aquaread Ltd., were applied for in situ measurement of the coordinates and the altitude of each site, as well as the water temperature, pH, water hardness  [50,51]) with locations of the studied waterbodies and indication of their type. The waterbodies are represented by numbers that follow those in Table 1.

Algal Identification and Counting by Light Microscopy
At each site, a water sample was collected for algal determination and counting by light microscopy (LM). The samples were taken from the surface layer (0-50 cm) in a volume of 0.5 L in case of visible blooms and of 1-1.5 L in cases of bright color of the water. The samples were immediately fixed with 2-4% formalin and transported to the lab, where they were sedimented to 30 mL for at least 48 h [38][39][40][41][42][43].
The taxonomic LM work was performed twice for all samples: (i) almost immediately after the collection on a Motic BA microscope with a Moticam 2000 camera, supported by the Motic Images 2 Plus software program; (ii) some months later, all samples were processed in a repetitive and comparative way on a Motic B1 microscopes supplied by a Moticam 2.0 MP camera with the Motic Images 3 Plus software program. Here, we note that the identification and counting was done by the same person (MPSG), which ensured the consistency of the LM data.
Algae were counted on a Thoma blood-counting chamber, with a minimum of four iterations for each sample with the cell taken as the main counting unit and a further estimation of the biomass [10,[38][39][40][41][42]. The relative abundance of the species was expressed according to the following modification of the Starmach scale [59] in comparison with species' contribution to the biomass [60]: "rare species" were those seen as single specimens in the whole microscopic slide (<0.5% of the biomass), "occasional species" those represented by up to five specimens (<5% of the biomass), "common, or abundant species" those seen with 6 to 30 specimens in a slide (5-20% of the biomass), whereas dominants and subdominants were evaluated among the most numerous species which contributed to >20 and >25% of the biomass, respectfully.

Analysis of Phytoplankton Marker Pigments
For the estimation of the general phytoplankton composition and relative phytoplankton biomass, HPLC was applied for marker pigment analysis following the standard operational procedure SOP5 described by [61]. Phytoplankton samples in a volume of 0.5-1 L were filtered at the earliest possibility after collection through 0.45 cellulose filters Whatman NC45 ST/Sterile EO (Merck KGaA, Darmstadt, Germany). Pigments were extracted by two 15 min sonications in ice, separated by an overnight stay in darkness, at a temperature of 4 • C and the final application of 90% acetone. Afterwards, the samples were transported to the lab in plastic tubes in a box with dry ice. During this transportation, only 1 of 70 samples, the tube from the reservoir Sopot, was destroyed and, therefore, pigment data for this reservoir are not provided in the paper.
The pigment analysis was performed on a Waters HPLC system equipped with a photodiode array detector. Pigment concentrations were determined from calibration with chlorophyll and carotenoid standards (DHI, Denmark), and CHEMTAX was used for the calculation of the contribution of the main phytoplankton groups [38,[40][41][42][61][62][63][64]. The initial table of pigments, applied as a matrix, is provided in [38].
The chlorophyll a, measured by HPLC, was compared with its field measurement and used as an expression of total algal biomass for the assessment of the trophic status according to the OECD System [46] and of the ecological status according to the intercalibrations related with the WFD [47].

•
Molecular-genetic analysis for the identification of anatoxin producers Anatoxin-A (ATX) and its analogues, anatoxins (ATXs), are alkaloid neurotoxins released by more than 40 species of Cyanoprokaryota [65,66]. They are produced by eight ATX synthetase genes (ana genes) [67], among which anaB-anaG genes are common for different producing genera [68]. Therefore, the anaC gene was selected for the amplification by the set of the following primer sequences, F-ATGGTCAGAGGTTTTACAAG and R-CGACTCTTAATCATGCGATC [69], of the material extracted from the samples collected in 2021 in order to complete our data, obtained after the analysis of the samples from 2018 and 2019 [43].
DNA was extracted from the field samples through filtration performed on 0.45 cellulose filters Whatman NC45 ST/Sterile EO (Merck KGaA, Darmstadt, Germany). The extracted DNA was amplified following the procedure MyTaqHS Mix (Bioline), which included the 12.5 µL Tag Mix, 10 pmol (1 µL) primers (both straight and inverted) and 50 ng total DNA. A specified program was used for the incubation of the reaction mixtures in a QB-96 Thermal Cycler: 35 cycles of denaturation (each 10 s at 95 • C), annealing at 55 • C for 30 s, an extension for 30 s at 72 • C, followed by a final extension for 5 min at 72 • C.
GeneJET™ Thermo Scientific and Clone JET PCR kits (Thermo Fisher Scientific, Waltham, MA, USA) were used for the purification and cloning of the anaC PCR products, and the recombinant sequences were sent to Macrogen Inc. (Seoul, Republic of Korea) for Sanger sequencing with the same pJET primers. All resulting data were manually edited and initially analyzed using the Vector NTI 11.5 (Thermo Scientific) software package. The Mega 6.0. program [70], a BLAST [71] search in the National Centre for Biotechnology Information (NCBI) GenBank database [72] and the neighbor-joining method with 1000 bootstrap values were used for organizing the anaC sequences in a phylogenetic tree. The obtained sequences were deposited in the NCBI GenBank [72] under the accession numbers OQ3119995-OQ3200013 and OQ355032.

•
Molecular-genetic analysis for the identification of microcystin producers Microcystins are the best-known and most-studied cyanotoxins, produced by Cyanoprokaryota, considered as being the most widely spread toxins in freshwaters [25,73]. In this study, the amplification of the mcyA gene from the microcystin synthetase mcyA-J gene cluster [74] was applied to the samples from 2018 in order to complete our earlier investigations, in which the mcyB and mcyE genes were used [39,40,42]. The amplified region was 510 bp long, described from the toxic strains M. aeruginosa UWOCCPCC 7806 and M. aeruginosa UWOCCPCC 7820 [75]. The amplification was accomplished by the set of forward primer mcyA-102F-CGATGAACAAATCGGGCAATGGCA and reverse primer u-620R-TGCAAGTTTCGCACATCTCCAAGG following [76,77].
A specified manufacturer program was used for the incubation of the reaction mixtures in a QB-96 Thermal Cycler starting with the denaturation at 95 • C for 3 min, followed by 35 cycles of denaturation (each 30 s at 95 • C) and 30 s of annealing at 52 • C, an extension at 72 • C for 30 s with a final extension step lasting 5 min at 72 • C. The cloning and further steps coincided with those described above for anatoxin. The obtained sequences were deposited in the NCBI GenBank database [72] under the accession numbers OM525685-OM525722, and ON075818-ON075819.

Statistical Analysis
The statistical analysis was conducted using the records of all identified species organized by main taxonomic groups (phyla and classes) and their abundance (rare, occasional, abundant, subdominant and dominant species) in the studied waterbodies. All records were encoded for use by statistical software. Data processing was done by the crossdisciplinary tool SPSS version 19, developed by IBM [78] using descriptives, frequencies and crosstabs, which aimed to prove the relations between the taxonomic groups and environmental parameters.
In almost all waterbodies, chlorophytes were the main contributors to the biodiversity, followed by cyanoprokaryotes ( Figure 3). An exception was the phytoplankton of the small mountain reservoir Beglika, in which algae from these two phyla were not found by conventional LM. Cyanoprokaryotes were not found by LM in two other waterbodies-in the large reservoir Suedinenie and in the small reservoir Krapets ( Figure 3).    Most of the algal taxa (421, or 53%) were found in a single waterbody and the number of species found in more than five waterbodies was much lower-63, or 8%. The same trend was valid for the species from each of the recorded phyla ( Figure 4b). The most widely spread algae belonged to chlorophytes: Tetraedron minimum (25 sites), followed by Coelastrum astroideum (17), Nephrochlamys subsolitaria (14), Golenkinia radiata and Oocystis lacustris (each in 13 sites), Monactinus simplex and Tetradesmus lagerheimii (Syn. Scenedesmus acuminatus) (each in 12 sites). The most spread cyanoprokaryote was Planktolyngbya limnetica (14 sites), followed by Microcystis wesenbergii (12 sites), Microcystis aeruginosa and Raphidiopsis raciborskii (each in 11 sites), Aphanizomenon klebahnii and Coelomoron pusillum (each in 10 sites). The most widespread species from other taxonomic groups in descending order of findings were the streptophyte Cosmarium neodepressum var. planctonicum and the pyrrhophyte Parvodinium elpatiewskyi (each found in 12 sites), followed by the ochrophytes Lindavia comta (11 sites) and Aulacoseira granulata (10 sites), as well as the euglenophyte Trachelomonas volvocina (10 sites).
According to the available Bulgarian algological literature, out of all 787 species, at least 53 (7%), recorded for first time in the country, can be considered alien. Most of them were observed as rare species in a small number of waterbodies. The exceptions were: (i) the tropical cyanoprokaryotes Raphidiopsis acuminato-crispa and R. gangetica, which codominated in the small inland reservoir Mechka together with R. raciborskii, found earlier in the country [81][82][83][84][85][86][87][88][89][90][91]; (ii) the North-Asian cyanoprokaryote Aphanizomenon yezoense, described as being from Japan [92] but currently spread also in Northern and Central Europe [57], which dominated in the small reservoir Studena and was subdominant in the coastal natural lake Durankulak; (iii) the chlorophyte Tetrallantos lagerheimii, described as being from Sweden [93] but afterwards recorded on different continents except the Antarctic [57] and currently found as dominant in the small inland reservoir Hadzhidimovo ( Figure 5). Regarding the non-native, allochthonous species, we would like to note that during this study, the invasive R. raciborskii ( Figure 5) was found as abundant in 11 waterbodies, where in 9 of them, it was recorded for the first time (i.e., Byalata Prust, Kaynaka, Eleshnitsa, Malka Smolnitsa, Mechka, Preselka, Shabla, Tsonevo and Uzungeren).    (e) (f) Figure 5. Alien phytoplankters in Bulgarian waterbodies: (a) Raphidiopsis raciborskii (straight trichome, thick green arrow points to its akinete, thin green arrow points to the heterocyst) and Raphidiopsis gangetica (white arrow points to its coiled trichome) in the inland microreservoir
borskii raciborskii in Mechka (arrow indicates its pointed heterocyst); (d) Raphidiopsis acuminato-crispa in Mechka (thin green arrow points to the heterocyst, thick green arrow points to the akinete); (e) Aphanizomenon yezoense in the coastal lake Durankulak-aggregation of trichomes in a fascicle (long white arrow points to the akinete, short white arrow points to the long transparent apical cell, and green arrow points one of the heterocysts in the fascicle); (f) Tetrallantos lagerheimii-coenobium of bent cells from the small inland microreservoir Hadzhidimovo, black scale-5 µm, relevant to all figures.

Algal Blooms and Toxic Species
According to the drone observations, supported by conventional LM studies and the HPLC analysis of marker pigments, during the three summers of investigation, blooms of cyanoprokaryotes occurred in the microreservoirs Birgo, Duvanli, Izvornik 2, Malka Smolnitsa, Mechka, Mogila, Nikolovo, Plachidol 2, Poroy, Preselka, Sinyata Reka, and Studena, in the large reservoir Mandra, as well as in the coastal lakes Burgasko Ezero and Durankulak (Figures 3, 6 and 7 and details in [38][39][40][41][42][43][44]). A relatively high contribution of cyanoprokaryotes to the biomass was detected in the reservoirs Krapets and Dospat (Figure 6), for which chlorophyll a data clearly showed a lack of blooms and a low trophic state (Figure 7). Similar was the case of the large mesotrophic reservoir Shilkovtsi, in which blooms were not seen and the relatively high contribution of cyanoprokaryotes was explained by the identification of marker pigments of Synechococcus type T1, typical for picoplankters which cannot be detected by conventional LM [42].  Table 1. Asterisks indicate that the real value of chl a in Izvornik 2 was 765 µg L −1.

Algal Blooms and Toxic Species
According to the drone observations, supported by conventional LM studies and the HPLC analysis of marker pigments, during the three summers of investigation, blooms of cyanoprokaryotes occurred in the microreservoirs Birgo, Duvanli, Izvornik 2, Malka Smolnitsa, Mechka, Mogila, Nikolovo, Plachidol 2, Poroy, Preselka, Sinyata Reka, and Studena, in the large reservoir Mandra, as well as in the coastal lakes Burgasko Ezero and Durankulak (Figures 3, 6 and 7 and details in [38][39][40][41][42][43][44]). A relatively high contribution of cyanoprokaryotes to the biomass was detected in the reservoirs Krapets and Dospat (Figure 6), for which chlorophyll a data clearly showed a lack of blooms and a low trophic state (Figures 7). Similar was the case of the large mesotrophic reservoir Shilkovtsi, in which blooms were not seen and the relatively high contribution of cyanoprokaryotes was explained by the identification of marker pigments of Synechococcus type T1, typical for picoplankters which cannot be detected by conventional LM [42].
Most detected blooms, despite their different intensity, supported the development of microcystin-, anatoxin-and microviridin-producing species and of the different cyanotoxins (Table 3) with proved the natural water cytotoxicity [94] and demonstrated the effects of low cylindrospermopsin doses on the gastrointestinal human cells [95]. It has to be noted that toxic cyanoprokaryotes were also found in waterbodies without blooms at the moment of sampling, such as in Ezerets, Koprinka, Uzungeren and Zhrebchevo (Table  3 and Figures 6 and 7). Up to now, in the studied waterbodies, nodularins and their main producer, Nodularia, have not been found despite the conducted targeted microscopic, chemical and molecular-genetic analyses [39]. Although cylindrospermopsin was detected in Bulgarian waterbodies [38,96], the molecular-genetic studies also revealed that  Table 1. Asterisks indicate that the real value of chl a in Izvornik 2 was 765 µg L −1 .
Most detected blooms, despite their different intensity, supported the development of microcystin-, anatoxin-and microviridin-producing species and of the different cyanotoxins (Table 3) with proved the natural water cytotoxicity [94] and demonstrated the effects of low cylindrospermopsin doses on the gastrointestinal human cells [95]. It has to be noted that toxic cyanoprokaryotes were also found in waterbodies without blooms at the moment of sampling, such as in Ezerets, Koprinka, Uzungeren and Zhrebchevo (Table 3 and Figures 6 and 7). Up to now, in the studied waterbodies, nodularins and their main producer, Nodularia, have not been found despite the conducted targeted microscopic, chemical and molecular-genetic analyses [39]. Although cylindrospermopsin was detected in Bulgarian waterbodies [38,96], the molecular-genetic studies also revealed that the identified Raphidiopsis raciborskii, Raphidiopsis mediterranea and Chrysosporum bergii in our study did not contain the cyrJ gene responsible for its production [44]. Table 3. Toxins and toxin-producing cyanoprokaryotes in the considered Bulgarian waterbodies, sampled in 2018 and 2019. Legend: CPS-cylindrospermopsin; MC-microcystin, followed by the exact type (LR, RR or YR), MV-microviridin, followed by the letter indicating the specific type (A, B, C, etc.), SXT-saxitoxins; (?)-supposed toxicity based on a comparison of newly obtained genetic sequences with light microscopic data. Waterbodies are arranged by years in alphabetical order. Currently, by combining LM data and molecular-genetic studies based on anaC gene with 24 newly obtained sequences, the presence of anatoxin producing Cuspidothrix in the 2021 summer samples from Durankulak, Mechka, Nikolovo, Studena and Yunets was proved (Figure 8). A comparison of these results with our data from 2018 and 2019 ( Figure 8 and [43]) demonstrated the presence of toxic Cuspidothrix issatschenkoi in the samples from Durankulak, Nikolovo and Yunets, and suggested once more the potential toxicity of Cuspidothrix elenkinii (found in 2019 in Koprinka [43] and in 2021 in Yunets) and of Cuspidothrix tropicalis (found in 2018 in Burgasko Ezero, in 2019 in Sinyata Reka [43] and in 2021 in Studena). In Mechka, rarely, morphologically peculiar young nonheterocytous and sterile trihomes of Cuspidothrix were found. Due to a lack of reproductive and resting cells, akinetes, their morphological determination was unreliable. Molecular-genetic data separated the sequences from Mechka from all other identified Cuspidothrix strains. In this small reservoir, three different Raphidiopsis species (R. acuminato-crispa, R. gangetica, R. raciborskii) codominated and, considering the close phylogenetic position of both genera Cuspidothrix and Raphidiopsis (for details see [43]), a further analysis of more genes is needed for a clarification of the strains isolated from Mechka. The coincidence with sequences of Aphanizomenon sp. in the constructed phylogenetic tree was explained in detail in [43] as caused by the taxonomic separation of the genus Cuspidothrix and of its type species Cuspidothrix issatschenkoi, in particular, from the genus Aphanizomenon [52]. needed for a clarification of the strains isolated from Mechka. The coincidence with sequences of Aphanizomenon sp. in the constructed phylogenetic tree was explained in detail in [43] as caused by the taxonomic separation of the genus Cuspidothrix and of its type species Cuspidothrix issatschenkoi, in particular, from the genus Aphanizomenon [52].  During the PCR amplification of the mcyA gene, responsible for the microcystin synthesis, 47 sequences were obtained, 9 of which showed a 100% homology with strains in NCBI [72] and 38 had a 99% homology with them. Molecular-genetic studies based on mcyA gene outlined two clusters and four subclusters in the 2018 summer samples (Figure 9), which, in combination with the LM observations, confirmed the presence of microcystinproducing Microcystis as follows: Microcystis aeruginosa and Microcystis novacekii in Mandra During the PCR amplification of the mcyA gene, responsible for the microcystin synthesis, 47 sequences were obtained, 9 of which showed a 100% homology with strains in NCBI [72] and 38 had a 99% homology with them. Molecular-genetic studies based on mcyA gene outlined two clusters and four subclusters in the 2018 summer samples ( Figure  9), which, in combination with the LM observations, confirmed the presence of microcystin-producing Microcystis as follows: Microcystis aeruginosa and Microcystis novacekii in Mandra (cluster I), Microcytis botrys in Poroy (subcluster I of cluster II), Microcystis aeruginosa in Poroy, Mandra and Durankulak (subcluster II of cluster II), Microcystis aeruginosa in Poroy, Burgasko Ezero and Mandra (subcluster III of cluster II), Microcysis novacekii in Burgasko Ezero, Microcystis botrys in Durankulak, where Microcystis aeruginosa also occurred (subcluster IV of cluster II).

Figure 9.
Neighbor-joining phylogenetic tree based on nucleotides sequences from ten library samples and closest sequences retrieved after a BLAST [71] search in the NCBI database [72] with indication of their accession number. Bootstrap values are shown at branch points (percentage of 1000 trials). Legend for the abbreviations of the waterbodies: Man-Mandra; Dur-Durankulak; Vai-Burgasko Ezero (Vaya); Por-Poroy; Blu-Sinyata Reka (Blue River) with Arabic numerals after the abbreviation, indicating the exact site of sampling and number of the sequence. For the identical sequences obtained during this study, only one accession number received from NCBI [72] is provided in each cluster or subcluster: (i) OM525686 is representing the sequences from the reservoir Sinyata Reka and site 4 of lake Durankulak (Blu 1 and Dur (1) 4); (ii) OM525702 is relevant for the sequences obtained from sites 1 and 3 of the reservoir Mandra, Poroy (Man (1) 4, Man (3) 6 and Por 1); (iii) OM525709 is for sites 1-2 from Lake Burgasko Ezero (i.e., Vai (1-2) 4, and Vai (1-2) 9) and for two sequences from the reservoir Poroy (Por 2, Poroy 5); (iv) OM525721 represents the sequences from sites 3 and 4 of Lake Durankulak (Dur (3) 3 and Dur (4) 1).

Algal Groups and Environmental Variables-Results from Statistical Analysis
In the conducted statistical SPSS analysis [78], the species from each algal group found at certain environmental conditions were expressed as a percentage from all species of the relevant group. The first results from the data processing by the SPSS tool and the application of Cramer's V evaluation [80] of 1996 records of all algal taxa and their abundance showed different but insignificant correlations. Therefore, we decided to exclude all rare species, the presence of which in the waterbodies was considered as nonrepresentative due to their finding in single specimens and in single sites. The resulting correlations obtained in this way were of moderate significance (0.2 < effect size field) except those with pH, which showed low confidence. Most probably, the lack of strong significance in this case was due to the targeted sampling in mostly eutrophic and hypertrophic waters with an alkaline character. The results presented below concern only taxonomic groups that were significantly correlated with other environmental parameters (TN, TP, trophic status, water hardness, conductivity and altitude). They were obtained after conducting the SPSS analysis based on the common, abundant, subdominant and dominant species from all algal groups with the subsequent exclusion of classes and phyla that showed correlations of low confidence.
The significant negative correlations were found between four taxonomic groups and the exact chlorophyll a values, considered as a proxy of the trophic status; the occurrence of species from Euglenophyta, Streptophyta and Eustigmatophyceae increased with the rising trophic status, whereas Chrysophyceae demonstrated a clear preference for a lower trophicity ( Figure 10).  The occurrence of all main algal groups was correlated with the TN concentrations, and the SPSS analysis revealed the preference of most of the identified species for high water quality conditions with TN below 7 mg L −1 [4], with Eustigmatophyceae in particular concentrated in waters with a TN range of 0.4-7 mg L −1 , and only Cyanoprokaryota, Euglenophyta and Bacillariophyta were spread in waters with TN values ranging between 7 and 10 mg L −1 (Figure 11). Figure 10. Cumulative species number in main taxonomic groups, expressed as a percentage of the total within the relevant phyla and classes (SP%), in waterbodies with different chlorophyll a concentration (µg L −1 ) as an expression of the trophic state [37]. Figure 11. Cumulative species number in main taxonomic groups, expressed as a percentage of the total within the relevant phyla and classes (SP%), in waterbodies with different concentrations of total nitrogen [4,37]. Figure 11. Cumulative species number in main taxonomic groups, expressed as a percentage of the total within the relevant phyla and classes (SP%), in waterbodies with different concentrations of total nitrogen [4,37].
According to the SPSS analysis, the identified taxonomic groups were significantly reverse-correlated with different concentrations of the other important nutrient, TP, except Eustigmatophyceae ( Figure 12).  Although most species from all groups were found in lowland and plain waterbodies (0-500 m a.s.l.), the distribution of the following taxonomic groups was more specific according to the altitude location: xanthophyceans and eustigmatophyceaens were spread only in the lowland waterbodies (0-200 m a.s.l.), while pyrrhophytes and streptophytes occurred in all altitude groups but had a preference for lowland and plain waterbodies ( Figure 13).  Regarding water conductivity, we have to note that during the field studies we did not measure values below <10 µS cm −1 (typical for distilled water, Table 1), and the SPSS analysis conducted for the three other conductivity categories allowed us to reveal nine taxonomic groups that showed a significant reverse relationship with this parameter, among which Synurophyceae could be outlined as related with waters of lower conductivity ( Figure 14). In this way, only species of Bacillariophyceae and Eustigmatophyceae found in this study could be excluded from the search for potential indicators, as independent from the conductivity of the water. Regarding water conductivity, we have to note that during the field studies we did not measure values below <10 µS cm −1 (typical for distilled water, Table 1), and the SPSS analysis conducted for the three other conductivity categories allowed us to reveal nine taxonomic groups that showed a significant reverse relationship with this parameter, among which Synurophyceae could be outlined as related with waters of lower conductivity ( Figure 14). In this way, only species of Bacillariophyceae and Eustigmatophyceae found in this study could be excluded from the search for potential indicators, as independent from the conductivity of the water. Considering water hardness, the SPSS analysis revealed that the spread of the species of four phyla and three classes was significantly correlated with this variable (Figure 15). The number of species of most algal groups increased with the rise of water hardness, but only Cryptophyta showed a preference for very soft water and Eustigmatophyceae to very  Considering water hardness, the SPSS analysis revealed that the spread of the species of four phyla and three classes was significantly correlated with this variable (Figure 15). The number of species of most algal groups increased with the rise of water hardness, but only Cryptophyta showed a preference for very soft water and Eustigmatophyceae to very hard water (Figure 15). Considering water hardness, the SPSS analysis revealed that the spread of the species of four phyla and three classes was significantly correlated with this variable (Figure 15). The number of species of most algal groups increased with the rise of water hardness, but only Cryptophyta showed a preference for very soft water and Eustigmatophyceae to very hard water (Figure 15). Figure 15. Cumulative species number expressed as a percentage of the total within the phyla and classes (SP%), in different waterbodies according to the water hardness ( o dh), classified after [37]. Figure 15. Cumulative species number expressed as a percentage of the total within the phyla and classes (SP%), in different waterbodies according to the water hardness ( • dh), classified after [37].

Discussion
Results from the present study demonstrated a high phytoplankton diversity in the sampled waterbodies, which comprised 787 species from seven phyla with a clear predominance of the green algae with 330 species, or 42% from all identified taxa. The second taxonomically rich group was Cyanoprokaryota, represented by 160 species. All data obtained by the LM and HPLC studies indicated the generally high contribution of blue-green algae in the summer phytoplankton of the studied waterbodies, especially in eutrophic and hypertrophic ones. A comparison of the results from the LM observations on algal abundance and dominance with the HPLC data ( Table 2, Figures 6 and 7) once more demonstrated the reliability of the application of the HPLC analysis of marker pigments in rapid phytoplankton characterization for water quality assessment [38,42,61,62].
Although the use of dominants for indicative purposes has long been debated, focusing on them is supported by the fact that their dynamics is important for the community stability, and they enhance the evaluation of resources availability [34,[96][97][98]. In this study, blue-green algae dominated by 33 species in 60% of the sampled water bodies (Table 2). These data are consistent with the well-known summer dominance of cyanoprokaryotes in nutrient-rich waters (e.g., [11,18,99]). If such dominance in small, shallow, lowland and plain waterbodies can be taken as a normal seasonal event, finding the heterocytous cyanoprokaryote Dolichospermum planctonicum as a dominant in the highest (among the studied sites) large oligotrophic mountain reservoir Golyam Beglik (Table 2) can be considered as alarming for the potential decrease of its water quality. This finding is in accordance with previous observations on the enlarged spread of blue-green algae, and of their potentially toxic species in particular in our mountain reservoirs [100][101][102].
The phytoplankton quantitative structure revealed by the application of the HPLC marker pigment analysis combined with the use of chlorophyll a values as a proxy for trophic status showed that by contrast with the summer dominance of cyanoprokaryotes in nutrient-rich waters, green and most yellow-brown, pyrrhophyte, or euglenophyte algae dominated in the oligo-to mesotrophic waterbodies (Table 2, Figures 6 and 7). Since the water quality in such waters is traditionally considered as being better, we support the use of a lack of cyanoprokaryote dominants to rapidly indicate nonproblematic water quality in the case of single, snapshot samplings. In addition, we confirm our earlier opinion [10] that in water quality assessment and relevant ecological status of the waterbodies both autochthonous and allochthonous species have to be taken into account. This comes from our current results that 53, or 7% of the recorded species were alien, newly recorded in the country. Although most of them were rare, found in single specimens, a few occurred in dominant phytoplankton complexes: the green Tetrallantos lagerheimii and the cyanoprokaryotes Aphanizomenon yezoense, Raphidiopsis acuminato-crispa and R. gangetica. Since the last three species belong to well-known cyanotoxin-producing genera [66,73], their abundant development can be problematic, ensuring their future spread in the country, as it was earlier shown for the invasive Raphidiopsis raciborskii [89][90][91] and is supported by the newly obtained data from this study on the increases of its spread and abundance in the country.
The combined LM and molecular-genetic data provided here are in accordance with our previous results on the high genetic diversity of Microcystis in Bulgarian waterbodies [39,40,42]. They prove its toxicity, as suggested by us earlier for the species Microcystis novacekii in addition to the well-known toxicity of Microcystis aeruginosa [39,40,42]. With the current phylogenetic tree, based on the PCR amplification of the mcyA gene, we are the first to provide for Bulgaria genetic data on the presence of potentially toxic Microcystis botrys, identified also by LM in Durankulak and Poroy, and we genetically confirmed our earlier LM finding of Microcystis novacekii in Mandra and Burgasko Ezero [39]. The current PCR data, based on the anaC gene amplification from the 2021 summer phytoplankton samples confirmed the presence and relatively broad spread of three potentially toxic Cuspidothrix species in our waterbodies (mainly C. issatschenkoi, but also C. elenkinii and C. tropicalis) recorded in 2018 and 2019 [43]. They also indicated this finding in four more waterbodies (Durankulak, Nikolovo, Studena and Yunets) and revealed a yet unidentified Cuspidothrix sequence in the small reservoir Mechka. The diversity and wide spread of numerous toxigenic cyanoprokaryote strains has already been stressed as alarming for Bulgarian waterbodies and their water quality ( [37,38,49,85,[87][88][89]100,101], among others).
On one hand, the high phytoplankton biodiversity associated with the great variability from site to site (reaching 198 species in Durankulak) showed the phytoplankton sensitivity to water quality, but on the other hand, it complicated the identification of indicator species for its assessment. In order to try to identify taxa that reflected particular environmental parameters, we conducted an SPSS statistical analysis [78,80]. After obtaining the first results based on 1996 records of all taxa and their relative abundance, we had to exclude all rare species, which occurred in single specimens in a single waterbody. In this way, it was possible to demonstrate different responses of the algae from different groups to the environmental variables such as nutrients (TP, TN) and chlorophyll a as proxy of the trophic status, water hardness and conductivity, and altitude as well. After the exclusion of some groups whose correlations were statistically insignificant, we outlined that Chrysophyceae showed a preference for a lower trophic status, Bacillariophyceae were indifferent to the water conductivity and occurred in waters of high TN, Cryptophyta preferred more soft water, Eustigmatophyceae were indifferent to the water conductivity but were significantly correlated with the increased trophic status, TP, water hardness and lowland waterbodies, and Euglenophyta preferred waters of higher trophicity and TN concentration. These results may encourage further search for bioindicators from these taxonomic groups, and this is especially valid for Eustigmatophyceae, which showed significant correlations with most variables but up to now was almost neglected in water quality assessments. Most species of this group found in this study were recorded earlier by us as commonly occurring, with an increasing abundance in the summer periods in the coastal lake Durankulak during its ongoing eutrophication [103,104]. Although they never dominated, we believe that their increasing records and recent outlining by the SPSS analysis will sharpen the attention of phytoplanktonologists to this group.
Last, but not least, all study results strongly supported our earlier opinion about the successful application of remote vehicles in the studies of water quality based on phytoplankton diversity and its blooms in particular [38,45]. The usage of drones allowed us to quickly choose the representative sampling sites, and thus save time, efforts and fuel during the sampling process. Therefore, we strongly recommend the application of this method in future field studies related to rapid water quality assessments.