Long-Term Eutrophication and Dynamics of Bloom-Forming Microbial Communities during Summer HAB in Large Arctic Lake

: Harmful algal blooms (HABs) in arctic lakes are recent phenomena. In our study, we performed a long-term analysis (1990–2017) of the eutrophication of Lake Imandra, a large subarctic lake, and explored the biodiversity of bloom-forming microorganisms of a 2017 summer HAB. We performed a 16Sr rRNA metabarcoding study of microbial communities, analysed the associations between N, P, C, and chlorophyll concentrations in the lake water, and developed models for the prediction of HABs based on total P concentration. We have demonstrated that blooms in Lake Imandra occur outside of optimal Redﬁeld ratios and have a nonlinear association with P concentrations. We found that recent summer HABs in a lake occur as simultaneous blooms of a diatom Aulacoseira sp. and cyanobacteria Dolichospermum sp. We have studied the temporal dynamics of microbial communities during the bloom and performed an analysis of the publicly available Dolichospermum genomes to outline potential genetic mechanisms beneath simultaneous blooming. We found genetic traits requisite for diatom-diazotroph associations, which may lay beneath the simultaneous blooming of Aulacoseira sp. and Dolichospermum sp. in Lake Imandra. Both groups of organisms have the ability to store nutrients and form a dormant stage. All of these factors will ensure the further development of the HABs in Lake Imandra and the dispersal of these bloom-forming species to neighboring lakes.


Introduction
Lake Imandra is the largest freshwater lake located in the southwestern part of the Kola Peninsula, Russia. The lake has a complex natural shape and is subdivided into three major parts connected by narrow straits with strong currents (Figure 1). The lake has a surface area of 876 km 2 , an average depth of 16 m, and a maximum depth of 67 m and is used for hydropower generation. The lake catchment area (12,342 km 2 ) is crowded, with various mining sites dumping their wastewater into the lake and storing tailings nearby. The area is populated by approximately 80,000 people, and along with industrial effluents, the lake receives treated municipal wastewater from three settlements [1]. The apatite processing industry is a major source of phosphorous for Lake Imandra, while the mining industry (MI) and municipal wastewater (MWs) are sources of nitrogen. The average mineralization of the lake water has increased 2.3 times from 27.7 mg/L in 1950 to 64.2 mg/L in 2017 due to the excessive input of various soluble organic compounds as well as soluble inorganic compounds, including mineral macronutrients. The load distribution is not equivalent among the major parts of the lake [2]. According to our estimations, the most eutrophic part of the lake, Bolshaya Imandra (average mineralisation: 73.2 mg/L), receives 1233 tonnes of anthropogenic nitrogen (N) and 300 tonnes of phosphorous (P) annually. The two remaining parts of the lake, Jokostrovskaya Imandra (average mineralisation: 68.1 mg/L) The seasonal dynamics of the planktonic communities vary between different parts of the lake and correlate in general with the concentrations of macronutrients. For the least eutrophic areas (southern part of Jokostrovskaya Imandra and Babinskaya Imandra), the biomass of the phytoplankton reaches the summer maximum in July and is normally dominated by diatoms and/or Chlorophyta spp. [3]. Heavily eutrophic Bolshaya Imandra and the northern part of Jokostrovskaya Imandra have multiple planktonic biomass peaks formed by blooming diatoms, dinoflagellates, and Chrysophyta, Chlorophyta, and Cyanobacteria species. Summer blooms in Bolshaya Imandra and the northern part of The seasonal dynamics of the planktonic communities vary between different parts of the lake and correlate in general with the concentrations of macronutrients. For the least eutrophic areas (southern part of Jokostrovskaya Imandra and Babinskaya Imandra), the biomass of the phytoplankton reaches the summer maximum in July and is normally dominated by diatoms and/or Chlorophyta spp. [3]. Heavily eutrophic Bolshaya Imandra and the northern part of Jokostrovskaya Imandra have multiple planktonic biomass peaks formed by blooming diatoms, dinoflagellates, and Chrysophyta, Chlorophyta, and Cyanobacteria species. Summer blooms in Bolshaya Imandra and the northern part of Jokostrovskaya Imandra start with the succession of diatoms and is usually dominated by Aulacoseira islandica [1][2][3]. The proliferation of diatoms may be accompanied by an increase in the biomass of mixotrophic Dinobryon spp. (Chrysophyta), described as early spring bloom formers in oligo-and mesotrophic aquatic ecosystems [6][7][8]. Planktonic biomass during diatom blooms exceeds 20 g/m 3 and gradually decreases during the summer. Based on the results of manual taxonomic analysis, the earlier described blooms occurring at the end of July and the beginning of August are formed by dinoflagellates (Ceratium hirundinella, Peridinium aciculiferum, and Peridinium goslaviense), Chlorophyta (Mucidosphaerium spp., Micractinium spp., Paulschulzia spp., Eudorina spp., Hindakia spp., and Pseudopediastrum spp.), or cyanobacteria species (Dolichospermum spp., Anabaena spp., and Coelosphaerium spp.) [1][2][3]9].
Cyanobacteria in Lake Imandra demonstrate typical generalist species characteristics, as their blooming may suddenly occur in any part of the lake. Various cyanobacterial species are recorded in planktonic samples from May to October and normally yield between 0.25 and 4.3 g/m 3 planktonic biomass. The most intensive blooming was observed in Bolshaya Imandra and the northern part of Jokostrovskaya Imandra, with cyanobacterial biomass reaching 25-31 g/m 3 at the peak of blooming [3]. According to our earlier published manual taxonomic identifications, the most intensive blooms are formed by Dolichospermum lemmermannii, creating thick green films on the surface of the lake, and are accompanied by high fish mortality ( Figure 2). D. lemmermannii blooms may develop as soon as the water temperature exceeds 16 • C and predominantly begin in semi-closed bays with little or no current [1,3,9]. Blooming over large open areas occurs during wind-free periods and is relatively rare. September blooms take place at water temperatures of approximately 7-8 • C and may be formed by Anabaena contorta, Anabaena ellipsoides, Anabaena subcylindrica, Dolichospermum circinale, Dolichospermum planctonicum, and Dolichospermum spiroides, as identified by manual taxonomic analysis [1,9]. Jokostrovskaya Imandra start with the succession of diatoms and is usually dominated by Aulacoseira islandica [1][2][3]. The proliferation of diatoms may be accompanied by an increase in the biomass of mixotrophic Dinobryon spp. (Chrysophyta), described as early spring bloom formers in oligo-and mesotrophic aquatic ecosystems [6][7][8]. Planktonic biomass during diatom blooms exceeds 20 g/m 3 and gradually decreases during the summer. Based on the results of manual taxonomic analysis, the earlier described blooms occurring at the end of July and the beginning of August are formed by dinoflagellates (Ceratium hirundinella, Peridinium aciculiferum, and Peridinium goslaviense), Chlorophyta (Mucidosphaerium spp., Micractinium spp., Paulschulzia spp., Eudorina spp., Hindakia spp., and Pseudopediastrum spp.), or cyanobacteria species (Dolichospermum spp., Anabaena spp., and Coelosphaerium spp.) [1][2][3]9]. Cyanobacteria in Lake Imandra demonstrate typical generalist species characteristics, as their blooming may suddenly occur in any part of the lake. Various cyanobacterial species are recorded in planktonic samples from May to October and normally yield between 0.25 and 4.3 g/m 3 planktonic biomass. The most intensive blooming was observed in Bolshaya Imandra and the northern part of Jokostrovskaya Imandra, with cyanobacterial biomass reaching 25-31 g/m 3 at the peak of blooming [3]. According to our earlier published manual taxonomic identifications, the most intensive blooms are formed by Dolichospermum lemmermannii, creating thick green films on the surface of the lake, and are accompanied by high fish mortality ( Figure 2). D. lemmermannii blooms may develop as soon as the water temperature exceeds 16 °C and predominantly begin in semi-closed bays with little or no current [1,3,9]. Blooming over large open areas occurs during windfree periods and is relatively rare. September blooms take place at water temperatures of approximately 7-8 °C and may be formed by Anabaena contorta, Anabaena ellipsoides, Anabaena subcylindrica, Dolichospermum circinale, Dolichospermum planctonicum, and Dolichospermum spiroides, as identified by manual taxonomic analysis [1,9]. All prior studies of the phytoplankton in Lake Imandra were performed using traditional manual taxonomic identification techniques and thus provided limited data on the All prior studies of the phytoplankton in Lake Imandra were performed using traditional manual taxonomic identification techniques and thus provided limited data on the entire biodiversity and co-occurrence of the bloom-forming species. In this paper, we present the results of a long-term analysis of the eutrophication of the lake, associations of the planktonic biomass with environmental parameters and metabarcoding analysis of the typical recent summer bloom in Lake Imandra. The sampling station was placed in the vicinity of Jokostrov Island, which is located at the intersection of the two most eutrophic parts of the lake (Figure 1). We also report nonlinear regression models allowing the prediction of the HABs for Lake Imandra based on the P concentrations in surface water samples. Our study is the first 16S rRNA metabarcoding study of a summer HAB conducted for a heavily eutrophic subarctic lake in Fennoscandia, and for the first time, we suggest nonlinear regression models as tools for the prediction of the HABs based in arctic lakes based on P concentrations in surface water samples.

Hydrochemical Analysis
Water samples were collected between 1990 and 2017. The sampling was performed according to the standardised national and international guidelines for the collection of water samples for environmental monitoring. Analysis of the hydrochemical parameters was performed at the ISO17025:2009 compliant laboratory certified by the Russian Federal Agency on Technical Regulation and Metrology (Certificate РОCC RU.0001.517126). The methods used in the current study are described in detail in respective ISO standards. The concentration of total organic carbon [10] was determined as described in ISO 8467:1993 (2019) by heating a sample in a boiling water bath with a known amount of potassium permanganate and sulfuric acid. The amount of consumed permanganate was determined by the addition of oxalate solution followed by titration with potassium permanganate. Concentrations of total phosphorus [11] and PO 4 3− ions were measured by the ammonium molybdate spectrometric method after potassium persulfate decomposition [11] and digestion with sulfuric acid (PO 4 3− ) as described in ISO 6878:2004(2019). The concentration of total nitrogen [12] was measured spectrometrically in the presence of Griess reagent after treatment of the samples with peroxodisulfate and subsequent cadmium reduction to NO 2 − (ISO 11905-1:1997(2019)). The concentration of NO 3 − was measured with a conductivity detector after ion exchange chromatography separation (ISO 10304-1:2007 (2016)). The concentration of NH 4 + was determined as described in ISO 7150-1:1984 (2017). The concentration of chlorophyll-a (Chl-a) was determined spectrometrically after acetone extraction (ISO 10260:1992 (2017)).
Calculation of the lake critical phosphorus loading (CLP) level was performed according to the method suggested by Vollenweider [13]. The actual annual phosphorous loads (LP) were calculated according to the method proposed by Kirchner and Dillon [14] by using the average annual concentration of P in water for the respective year.

Prediction of Chlorophyll-a Concentrations and Visualization of TC/TN/TP/Chl-a Ratios
LOESS regressions [15] were built based on chlorophyll-a concentrations measured for sets of lake water samples collected in 2012. Modelling was performed in the R environment [16], and results were visualized using the ggplot2 package [17]. The LOESS (LR) models were used to predict Chl-a levels based on total P concentrations observed in 1990-2011 and 2013-2017. Both the predicted and measured Chl-a concentrations were used for the visualization of TC/TN/TP/Chl-a ratios in the ternary graphical system suggested by Smith et al. [18]. Linear regression analysis and PCA as well as visualization of the TC/TN/TP/Chl were performed in R [16] using integrated functions/ggplot2 visualization [17] and the ggtern package, respectively [19].

DNA Preparation and Processing
The surface water samples were collected from a boat using sterile one litre bottles (67.595249 N; 33.003494 E). The bottles were sealed and were transported to the shore, where in the laboratory, the 1 L water samples were filtered through standard line bottle top filters (VWR) with 0.22 µm membranes. Membranes were used for DNA extraction with a DNeasy PowerWater Kit (QIAGEN). The blooming intensity and the respective time points for sample collection were determined based on the total chlorophyll-a concentrations. The extracted DNA was quantified, accessed for purity, and used to generate 16S amplicons with primers specific to the V1-V3 hypervariable region (27F-AGAGTTTGATCCTGGCTCAG and 534R-ATTACCGCGGCTGCTGG). The library was prepared using the Nextera XT DNA Library Prep Kit (Illumina). The amplicons were analysed using an Illumina MiSeq sequencer in 2 × 300 bp paired-end read sequencing mode and a 600-cycle v3 MiSeq Reagent Kit. The raw sequencing reads are available from the NCBI through BioProject accession number PRJNA526000.

Analysis of Bacterial Communities
Analysis of microbial communities was performed using the Quantitative Insights into Microbial Ecology pipeline [20]. The raw sequence reads were trimmed and filtered by quality score to retain sequences with a base call > Q20. The taxonomic analysis was performed using the Greengenes 13.8 16S rRNA database [21]. The naive Bayes classifier [22] was trained using the V1-V3 primer pair 27F (AGAGTTTGATCCTGGCTCAG) and 534R (ATTACCGCGGCTGCTGG), and the same primers were used for the generation of amplicons. For general differential abundance analysis, an operational taxonomic unit (OTU) table was created using the DADA2 pipeline [23]. Differential abundance analysis of the filtered table was performed in QIIME2 using the Gneiss algorithm [10]. The analysis of correlations between the groups of microorganisms and environmental factors was performed in R [16], and the results were visualized with ggplot2 [17].

Analysis of the Publically Available Dolichospermum Genomes
Analysis of the NCBI available Dolichospermum genomes (NZ_AP018316.1 and NZ_ CP043056.1) was performed locally using the BLAST+ algorithm and a key procaryotic B 1 and B 12 vitamins biosynthesis enzyme sequence as a query. Statistically significant matches were curated manually and were verified against the UniProt database to confirm the predicted biosynthetic functions.

Eutrophication Dynamics
Lake Imandra has a complex snow/groundwater/river water feeding regime. The contribution of each factor greatly depends on the season. In the winter, the lake is mainly fed by groundwater and rivers, while in spring, snowmelt has a greater impact on its hydrological regime. The only outlet from the lake is used for hydropower generation, and as a result, the lake has a relatively slow water exchange rate (Table 1). It takes 2.38 years to completely exchange the water in the lake, which, in turn, makes the lake highly sensitive to the excessive input of macronutrients.
The general decline in industrial activities in 1986-1993 had positive effects on the eutrophication of the lake. The study site used for metabarcoding analysis was located at the intersection of the two most eutrophic parts of the lake. Between 1990 and 2017, the annual critical load of anthropogenic P was not exceeded in 1990-1991. In 1992-2017, the annual critical phosphorous load exceeded the critical load by 1.1-3.7 times ( Figure 3). The annual average concentrations of the main macronutrients in the lake water in 1990-2017 exceeded the initial parameters of the lake water recorded in 1930, as it contained 0-35 µg/L NO 3 − and 0-8 µg/L PO 4 3− [5]. The main physicochemical parameters of the surface water vary in different areas of Lake Imandra and correlate with the localisation of the eutrophication sources, river inlets, and the directions of the currents. The NH 4 /NO 3 balance is typically shifted towards NH 4 ions in Jokostrovskaya Imandra and reflects the decomposition of the organic matter in the sediment. NO 3 ions are dominant in water samples collected from Bolshaya Imandra, which receives the mining wastewater enriched with both NO 3 and P. N and P concentrations both have seasonal and annual variations and depend on the intensity of the mining and enrichment activities of API. In the summer months, the concentrations of the main macronutrients are greatly affected by proliferation of the photosynthetic microorganisms actively consuming excessive N and P from the lake water. The average physicochemical parameters of the surface water of different parts of the Lake ( Table 2) remained relatively stable over the years, as no new API processing plants or mines were opened in the area in 1990-2017. Table 1. The morphometric characteristics and critical and actual annual phosphorus loadings of Lake Imandra and its parts. H max -the maximum lake depth; Z-the average lake depth; S lake -the lake surface area; V-the lake volume; S cathc -the catchment area; R-the water drainage rate; t-the time of water discharge; CLP-critical phosphorus loading; Q (P)max -the maximum permissible phosphorus admission to the lake; M (P) -the maximum permissible drainage module of phosphorus; L (P) -annual phosphorous load in 2017. Year) The annual average concentrations of the main macronutrients in the lake water in 1990-2017 exceeded the initial parameters of the lake water recorded in 1930, as it contained 0-35 µ g/L NO3 − and 0-8 µ g/L PO4 3− [5]. The main physicochemical parameters of the surface water vary in different areas of Lake Imandra and correlate with the localisation of the eutrophication sources, river inlets, and the directions of the currents. The NH4/NO3 balance is typically shifted towards NH4 ions in Jokostrovskaya Imandra and reflects the decomposition of the organic matter in the sediment. NO3 ions are dominant in water samples collected from Bolshaya Imandra, which receives the mining wastewater enriched with both NO3 and P. N and P concentrations both have seasonal and annual variations and depend on the intensity of the mining and enrichment activities of API. In the summer months, the concentrations of the main macronutrients are greatly affected by proliferation of the photosynthetic microorganisms actively consuming excessive N and P from the lake water. The average physicochemical parameters of the surface water of different parts of the Lake ( Table 2)

Nutrient Stoichiometry and Chlorophyll-a
The training datasets for the generation of the LR regression models were based on total phosphorous concentrations and concentrations of chlorophyll-a ( Figure 4) measured for the same sampling points during July 2012. The LR curves describing the relationships between the TP and Chl-a concentrations have dual Chl-a concentration peaks that are specific to each part of the lake and represent the blooming events that took place in July 2012 (Figure 5a). The sampling site was placed at the intersection of the two most eutrophic parts of the lake to cover the entire range of chlorophyll and TP concentrations. The PCA (Figure 5b) also demonstrated three distinct TP and Chl-a concentration distribution patterns for different parts of the lake, thus supporting the results of the LOESS regression modelling.

Nutrient Stoichiometry and Chlorophyll-a
The training datasets for the generation of the LR regression models were based on total phosphorous concentrations and concentrations of chlorophyll-a ( Figure 4) measured for the same sampling points during July 2012. The LR curves describing the relationships between the TP and Chl-a concentrations have dual Chl-a concentration peaks that are specific to each part of the lake and represent the blooming events that took place in July 2012 (Figure 5a). The sampling site was placed at the intersection of the two most eutrophic parts of the lake to cover the entire range of chlorophyll and TP concentrations. The PCA (Figure 5b) also demonstrated three distinct TP and Chl-a concentration distribution patterns for different parts of the lake, thus supporting the results of the LOESS regression modelling.  Regular cyanobacterial blooming in Lake Imandra started in 2000, while irregular blooming had been observed prior to 2000, since 1993. Plotting TOC, TP, and TN alongside Chl-a concentrations for the period between 1993 and 2017 demonstrated that the Chl-a concentrations followed the TP and NT concentrations (Figure 6b). Correlation analysis of associations between the two hydrochemical parameters and Chl-a concentrations supported this observation, as TP/Chl-a had a strong positive correlation (r = 0.86, p < 0.05), and TN/Chl-a had a medium positive correlation (r = 0.66, p < 0.05). The associations between TOC and Chl-a were not significant.    In the ternary diagram, the parameters of all of the studied samples are located away from the optimal Redfield value, and thus, a balanced C/N/P ratio is not needed for the development of cyanobacterial blooms in Lake Imandra. All points are clustered in the In the ternary diagram, the parameters of all of the studied samples are located away from the optimal Redfield value, and thus, a balanced C/N/P ratio is not needed for the development of cyanobacterial blooms in Lake Imandra. All points are clustered in the same area of the ternary diagram and to a great extent only differ in chlorophyll  (Figure 6b).

Microbial Communities
In the current study, we used primers in the V1-V3 16S rRNA region, allowing for the amplification of bacterial and chloroplast 16S rRNA sequences. The amplicons obtained from eukaryotic photosynthetic microorganisms (EPMs) were included in the analysis, and the relative abundances of the OTUs were calculated considering the relative abundances of the OTUs derived from both the prokaryotes and chloroplasts of EPMs.
We identified 37 bacterial phyla among the whole set of the studied samples. The prokaryotic communities were dominated by Proteobacteria, Bacteroidetes, Actinobacteria, Verrucomicrobia, Planctomycetes, and Cyanobacteria followed by members of other less abundant taxa. A total of 18 of 37 prokaryotic phyla comprised unclassified microorganisms having no significant 16S rRNA similarities with the cultured and characterised bacterial species. An analysis of the alpha and beta diversity performed for both prokaryotes and EPMs showed statistically significant changes (p < 0.05) in the biodiversity and relative abundances of the species at the peak, middle, and terminal blooming phases identified based on the total chlorophyll-a concentrations. Shannon's H index dropped three times between the peak and the end of blooming ( Figure 7a) through a dramatic decrease in the relative abundances of the EPM species. same area of the ternary diagram and to a great extent only differ in chlorophyll concentrations. The blooming in 1993-2017 occurred under N-and P-limited conditions for plankton if judged based on the Redfield optimal ratio (Figure 6b).

Microbial Communities
In the current study, we used primers in the V1-V3 16S rRNA region, allowing for the amplification of bacterial and chloroplast 16S rRNA sequences. The amplicons obtained from eukaryotic photosynthetic microorganisms (EPMs) were included in the analysis, and the relative abundances of the OTUs were calculated considering the relative abundances of the OTUs derived from both the prokaryotes and chloroplasts of EPMs.
We identified 37 bacterial phyla among the whole set of the studied samples. The prokaryotic communities were dominated by Proteobacteria, Bacteroidetes, Actinobacteria, Verrucomicrobia, Planctomycetes, and Cyanobacteria followed by members of other less abundant taxa. A total of 18 of 37 prokaryotic phyla comprised unclassified microorganisms having no significant 16S rRNA similarities with the cultured and characterised bacterial species. An analysis of the alpha and beta diversity performed for both prokaryotes and EPMs showed statistically significant changes (p < 0.05) in the biodiversity and relative abundances of the species at the peak, middle, and terminal blooming phases identified based on the total chlorophyll-a concentrations. Shannon's H index dropped three times between the peak and the end of blooming ( Figure 7a) through a dramatic decrease in the relative abundances of the EPM species. Amplicons derived from photosynthetic microorganisms on average comprised 34.0% of the total microbial population during the blooming peak, 18.6% during the middle of the blooming period, and only 1.69% at the end of the blooming period. The photosynthetic communities at the peak and middle phases of the blooming period were dominated by EPMs that made up 24.92% and 16.98% of the total number of amplicons, respectively. At the end of the blooming period, the EPMs were found in only 3 of 10 studied samples on average, comprising 0.15% of the total number of amplicons. The death of the EPMs was accompanied by 18.62%, 0.76%, 0.12%, and 25.56% increases in the mean relative abundances of Proteobacteria, Firmicutes, Fusobacteria, and unclassified prokaryotic organisms, respectively, leading to significant changes in beta-diversity metrics ( Figure  7b). Chlorobi, Chloroflexi, and Planctomycetes detected in our samples had statistically significant associations with several cyanobacterial genera and EPMs (Figure 8). Their mean relative abundances at the peak of blooming, however, ranged between 0.07% and 0.8% of the total microbial population. Amplicons derived from photosynthetic microorganisms on average comprised 34.0% of the total microbial population during the blooming peak, 18.6% during the middle of the blooming period, and only 1.69% at the end of the blooming period. The photosynthetic communities at the peak and middle phases of the blooming period were dominated by EPMs that made up 24.92% and 16.98% of the total number of amplicons, respectively. At the end of the blooming period, the EPMs were found in only 3 of 10 studied samples on average, comprising 0.15% of the total number of amplicons. The death of the EPMs was accompanied by 18.62%, 0.76%, 0.12%, and 25.56% increases in the mean relative abundances of Proteobacteria, Firmicutes, Fusobacteria, and unclassified prokaryotic organisms, respectively, leading to significant changes in beta-diversity metrics (Figure 7b). Chlorobi, Chloroflexi, and Planctomycetes detected in our samples had statistically significant associations with several cyanobacterial genera and EPMs (Figure 8). Their mean relative abundances at the peak of blooming, however, ranged between 0.07% and 0.8% of the total microbial population. PCA showed that changes in the relative abundances of EPMs and eight Cyanobacterial genera (Dolichospermum, Pseudanabaena, Nostoc, Chloroidium, Paulinella, Synechococcus, Planktothrix, and Phormidium) made the largest contribution to the distribution of the data in two-dimensional coordinates (Figure 9). The pool of EPM amplicons comprised OTUs assigned to Bacillariophyta, Chlorophyta, Euglenozoa, Haptophyceae, Stramenopiles, and Streptophyta. EPM OTUs at the peak of blooming were dominated (24.1%) by a single amplicon derived from a diatom Aulacoseira sp. (Bacillariophyta). All of the remaining OTUs thus comprised 0.82% of all of the amplified sequences and therefore likely had no major impact on the environmental consequences of the bloom. In the middle of the bloom, the mean relative abundance of the dominant OTU dropped to 16.34%, and at the end of the bloom, this amplicon was not detected in any of the studied samples ( Figure  10). PCA showed that changes in the relative abundances of EPMs and eight Cyanobacterial genera (Dolichospermum, Pseudanabaena, Nostoc, Chloroidium, Paulinella, Synechococcus, Planktothrix, and Phormidium) made the largest contribution to the distribution of the data in two-dimensional coordinates (Figure 9). The pool of EPM amplicons comprised OTUs assigned to Bacillariophyta, Chlorophyta, Euglenozoa, Haptophyceae, Stramenopiles, and Streptophyta. EPM OTUs at the peak of blooming were dominated (24.1%) by a single amplicon derived from a diatom Aulacoseira sp. (Bacillariophyta). All of the remaining OTUs thus comprised 0.82% of all of the amplified sequences and therefore likely had no major impact on the environmental consequences of the bloom. In the middle of the bloom, the mean relative abundance of the dominant OTU dropped to 16.34%, and at the end of the bloom, this amplicon was not detected in any of the studied samples ( Figure 10).

Putative Genetic Triats Beneath Simultaneous Blooming
At present, the NCBI nucleic acid archive contains two fully sequenced genomes of Dolichospermum sp., both containing genes essential for the synthesis of B1 and B12 vitamins. The prokaryotic B1 biosynthesis pathway requires at least four enzymes. Searches

Putative Genetic Triats Beneath Simultaneous Blooming
At present, the NCBI nucleic acid archive contains two fully sequenced genomes of Dolichospermum sp., both containing genes essential for the synthesis of B1 and B12 vitamins. The prokaryotic B1 biosynthesis pathway requires at least four enzymes. Searches Figure 10. The relative abundance of pro-and eukaryotic photosynthetic species during various phases of HAB observed in Imandra Lake in July 2017.

Putative Genetic Triats Beneath Simultaneous Blooming
At present, the NCBI nucleic acid archive contains two fully sequenced genomes of Dolichospermum sp., both containing genes essential for the synthesis of B 1 and B 12 vitamins. The prokaryotic B 1 biosynthesis pathway requires at least four enzymes. Searches of the genomes of Dolichospermum compactum (NZ_AP018316.1) and Dolichospermum sp. (NZ_CP043056.1) showed that both species contain all four genes required for thiamine biosynthesis (ThiC > ThiD > ThiN > ThiL). The thiamine phosphate synthase ThiC is a key enzyme that is responsible for the biosynthesis of the aminopyrimidine ring of the thiamine molecule, and ThiL is responsible for the final stage of the biosynthetic process or formation of thiamine diphosphate from thiamine monophosphate. De novo biosynthesis of vitamin B 12 or cobalamin requires approximately 30 biosynthetic genes and may occur through the aerobic or anaerobic pathway. Several cobalamin biosynthesis genes are shared by two pathways, but both have enzymes specific to each pathway. Our analysis has shown the presence of homologs of at least five key cobalamin biosynthesis enzymes in both available genomes. The CobE protein participates in the conversion of cobalt-precorrin 5 into cobalt-precorrin 6; CobD is involved in cobalamin biosynthesis and is likely responsible for conversion of adenosylcobyric acid to adenosylcobinamide or adenosylcobinamide phosphate; CobW is involved in the storage of cobalt ions prior to their use in cobalamin biosynthesis; the CobQ protein catalase amidations at positions B, D, E, and G in adenosylcobyrinic A,C-diamide intermediates; and the CbiM protein is a part of the transporter complex CbiMNOQ, which is involved in cobalt import.

Selection of the Primers
The current version of Illumina dye sequencing technology allows the amplification of 2 × 300-350 bp DNA fragments; thus, the precision of microbiome analysis greatly relies on the selected primers. The microbial 16S rRNA gene has approximately 1600 base pairs and consists of nine hypervariable regions. Amplification of conserved regions the allows reliable identification of higher taxa, while analysis based on quickly evolving areas may provide better discriminating power at the genus or even species level. A number of studies have shown that the selection of the area for amplification influences the outcome of taxonomic identification [24,25]. At the same time, the selection of primer pairs did not affect the ecological interpretation of the results for planktonic communities, as the identified taxa had similar correlation patterns for all primer pairs [26]. A study by Parulekar et al. comparing the performance of the V1-V3 and V3-V4 primer pairs for eutrophic lakes in southern Norway showed that the two primer pairs equally identified Cyanobacteria, Proteobacteria, Actinobacteria, Acidobacteria, and Parcubacteria, while Bacteroidetes, Verrucomicrobia, Armatimonadetes, Firmicutes, and Planctomycetes were better represented in the samples amplified with the V3-V4 pairs [27]. Chloroplasts are of cyanobacterial origin, and thus, prokaryotic primers adequately amplifying cyanobacterial 16S rRNA sequences provide good insight into the biodiversity of eukaryotic species containing chloroplasts [28]. Due to the variations of the chloroplast copy numbers [11] the approaches based on 16S rRNA are not suitable for the absolute quantification of the OTUs; at the same time, the use of relative abundances of the OTUs allow the reliable and statistically significant comparison of the various phases of the blooming events. Our study was focused on the blooming dynamics of photosynthetic organisms in general (both proand eukaryotic), and thus, we selected V1-V3 primer pairs for the generation of amplicons.

Concentrations of Macronutrients and Chlorophyll-a
The results of hydrochemical and chlorophyll-a concentration analysis demonstrate that these parameters greatly vary among the parts of the lake. In both Bolshaya Imandra and Babinskaya Imandra, the NH 4 + /NO 3 − balance was shifted towards NO 3 − ions, while in Jokostrovskaya Imandra, the NO 3 − ions were depleted. This can be explained by ongoing diatom blooms. Comparative studies of the uptake of NH 4 + /NO 3 − have shown that diatoms prefer NO 3 − uptake over NH 4 + assimilation and have developed efficient mechanisms for the storage and reduction of NO 3 − ions [29,30]. The reduction of the stored NO 3 − ions to ammonium is used by diatoms to survive dark and anoxic conditions in the bottom sediment and may occur through two groups of metabolic pathways. The first group comprises the so-called dissimilatory pathways that reduce nitrate to ammonium, and the second group consists of several pathways that reduce nitrate to NO 2 , N 2 O, or N 2 . NO 3 − uptake by diatoms occurs strictly in the presence of oxygen and thus takes place in oxygenated subsurface layers. This correlates with the depleted NO 3 concentrations in the surface water layers and with the increased concentrations near the bottom. The sinking of the diatoms may thus contribute to the vertical transport of accumulated NO 3 ions to bottom sediment [31]. The increase of chlorophyll-a concentrations in the water samples normally begins in early June and correlates with the proliferation of various Diatom or Dinobryon species. The contribution of cyanobacteria to a total planktonic biomass is insignificant at the beginning of the summer, as their biomass rarely exceeds 0.1 g/m 3 . The extensive proliferation of the cyanobacteria begins at the end of July and may reach 31 g/m 3 during the August or September blooms if the conditions are favourable [2,3]. In the current study, we demonstrated for the first time that Diatoms and cyanobacteria may form a combined bloom in an arctic lake.

HABs in Arctic Waters
The algal blooming in arctic lakes is a recent phenomenon. In temperate climates, the blooming may be accompanied by the release of the toxins or may be formed by a non-toxin producing species. The blooms in Lake Imandra are formed by both pro-and eukaryotic microalgae or their simultaneous proliferation. All blooming events recorded by our group in Lake Imandra in 1990-2017 coincided with fish death of various intensities. The magnitude of the mortality was the highest during cyanobacterial blooms, while blooming of EPMs had a less noticeable impact on the fish population in general. The majority of cyanobacterial communities in pristine arctic waters exist in the form of algal mats attached to abiotic substances or moss filaments. Suspended native cyanobacterial species predominantly belong to Picocyanobacteria, whose surface-to-volume ratio allows them to form up to 30-60% of the planktonic biomass in oligotrophic arctic waters [32]. Bloom-forming species such as Anabaena, Microcystis, and Aphanizomenon are considered invasive species in arctic and are only present in eutrophic water basins [32]. The stoichiometric combination 106C:16N:1P, referred to as the Redfield ratio, is considered to be optimal for algal proliferation [18]. We have demonstrated that blooming in Lake Imandra occurs outside of the optimal Redfield ratio and takes place under conditions that are defined as N-and P-depleted for plankton. Imandra Lake undergoes the complex blooming, as the proliferation of diatom Aulacoseira sp. is accompanied by the intensive increase of biomass of the diazotrophic cyanobacterium Dolichospermum sp. This correlates with the findings of other authors also showing that cyanobacterial HABs may occur under both N-and P-limited conditions. In this case, the blooming species rely on N 2fixation and P-scavenging pathways [33]. The availability of easily accessible nutrients has a significant impact on the environmental consequences of the cyanobacterial blooms. Experiments conducted on microcystin-producing Planktothrix species have shown higher levels of toxin production in response to the addition of dissolved N and the highest levels have occurred in response to NH 4 + /urea + PO 4 − combinations [34]. Increasing levels of N and P in the water of Lake Imandra will thus further promote HAB development, and increased fish mortalities should be expected in future. Individual cyanobacterial genera and strains respond differently to various N/P combinations; thus, the selective management of these macronutrients [35,36] combined with the management of blooming communities are required to control cyanobacterial blooms in Lake Imandra and other artic lakes. Controlling of the N/P concentrations alone will likely have a minor impact on the blooming, as cyanobacteria have developed nutrient scavenging mechanisms as well as mechanisms for the accumulation of storage lipids and carbohydrates [37,38]. Peculiarities in metabolism in combination with a simple genetic apparatus lead to high mutation rates, which, in turn, constantly create new genetic forms of these microorganisms [39,40]. In addition to ordinary binary fission, some species of photosynthetic bacteria are able to form baeocytes [41,42], which are often motile and are used for dispersion into new niches [41]. Some cyanobacteria form dormant akinetes, allowing them to survive for decades while experiencing starvation in lake sediment or to travel long distances following airflows or various carriers [43]. The measures for controlling baeocytes and the germination of dormant akinetes from lake sediment should also be considered for the efficient restoration of the ecosystem of the Lake Imandra.

Putative Genetic Triats and Adaptations Beneath Simultaneous Blooming
Dolichospermum is adapted to growth at low temperatures, has a low demand for nutrients and is able to accumulate phosphates [44]. Trichomes of Dolichospermum can migrate within the euphotic zone in search of optimal photosynthesis and nutrient gradients, and thus, these species will likely develop blooms if present in the environment [45]. Dolichospermum spp. form two types of nonvegetative akinetes [46] surrounded by an additional cell envelope and contain storage glycogen [47] and cyanophycin granules [48]. The early summer akinetes are formed soon after the extensive bloom and may cause subsequent blooms within the same season if the conditions are favourable [46].
Repeated summer blooming was observed in Lake Imandra starting in 2000, indicating that the formation mechanisms of summer akinetes are active and are used by Dolichospermum sp. under arctic conditions. The formation of late-summer akinetes [43,49] is triggered by light intensity and is used for the long-term survival and formation of next year blooms in Imandra Lake. The late summer akinetes of Dolichospermum are able to survive in sediment for decades. Salmaso et al. successfully germinated akinetes of Dolichospermum from 1989 sediment cores and determined that the establishment of the species coincided with the beginning of a rapid increase in total phosphorus concentration in the studied lake [45]. The germination of dormant akinetes does not require external N and P, as these elements are initially supplied by stored cyanophycin and glycogen. Studies conducted on the akinetes of two Dolichospermum species blooming in warmer climates have shown that the optimal temperature for their germination is approximately 20-25 • C. Incubation in this temperature interval synchronized the germination of akinetes and directly resulted in blooming [50]. Temperatures of 20-25 • C are rarely observed, even in the surface water layers of Lake Imandra, and have never been recorded for the bottom layers. We have not found any studies describing the low-temperature germination of Dolichospermum sp. to date, but cyanobacteria with low germination temperature optima indeed exist. Anabaena flos-aquae flourishes at water temperatures from 10 to 29 • C, with the most active akinetes germinating between 8 and 12 • C [51]. Light intensity is the second factor that triggers germination. The abundance of akinetes in the uppermost centimetre of the sediments of profundal sites is approximately 20 times higher than in the littoral zone, thus indicating that germination occurs in the shore areas. The akinetes formed over deep areas or that drift to profundal sediments are excluded from the seed bank [52]. Under normal conditions, the deep sediment layers in Lake Imandra are never exposed to sufficient light and temperature, but the reduction of the water level due to hydropower generation returns deposited akinetes into the active pool and promotes blooms in the lake.
Diatoms are surrounded by extracellular silicon dioxide shells [53]. Plastids of autotrophic diatoms lack thylakoid subdomains and have a random distribution of photosystems I and II, providing direct energy transfer between the photosynthetic complexes. The photosynthesis efficiency is thus higher in diatoms than in the chloroplasts of plants [54]. Diatoms quickly absorb nitrogen in NO 2 , NH 4 + , NO 3 − , and organic forms though various transporters but have no N 2 -fixing apparatus and thus rely on intracellular N recycling or symbiotic relations under N-limited conditions [30]. The phenomenon of diatom diazotroph associations (DDAs) has been widely described for marine diatoms, as nitrogen starvation is more prevalent in marine environments [55]. N 2 -fixing bacteria may occur as ectosymbionts, partial symbionts or endosymbionts, and their location is regulated by the host [56]. The benefits of symbiosis with prokaryotes are not limited to the supply of N-rich substrates, as diatoms are dependent on external vitamin sources [57,58]. It has been shown that Daphnia actively grazing on Aulacoseira spp. [59] and thus simultaneous blooming with toxin producing Dolichospermum species might provide additional advantages protecting Aulacoseira spp. from potential grazers [12,60]. The lifecycle of Aulacoseira spp. does not include the formation of dormant spores, but instead, the cells enter the altered resting state with reduced chloroplasts. The resting state as well as revitalisation into a fully active state is triggered by the light conditions and may take as little as 8 h in Aulacoseira subarctica. The turbidity of the water and the thickness of the water level over the sediment containing the resting cells are thus important for the revitalisation of the Aulacoseira spp. into the active state [61].
We have studied repeated summer HABs in Lake imandra since their beginning in 2000 using manual taxonomic techniques. While both species were detected in water samples in different time points, we were not able to depict fine interplay between the diatom Aulacoseira sp. and cyanobacteria Dolichospermum sp. during the earlier HABs due to the limitations of the manual methods. Advances in the development of molecular tools not only allowed the discovery of the interplay between two species, but also allowed us to examine their genomes. The discovery of the presence of the key vitamin B 1 and B 12 biosynthetic genes in the genomes of Dolichospermum sp. reported in our study as well as the availability of experimental evidence on the symbiotic relations between diatoms and prokaryotes published by other authors allow us to suggest that the simultaneous succession of Aulacoseira sp. and Dolichospermum sp. during summer HABs in lake Imandra is not incidental and relies on symbiotic relationships between the two species.

Conclusions
In the current study, we performed a long-term analysis of the eutrophication dynamics of Lake Imandra, a large subarctic lake, and performed metabarcoding analysis of a typical summer HAB accompanied by high fish mortality. Our data demonstrate that the lake is enriched with nutrients and that the annual input of macronutrients from the API and municipal wastewater exceeds the theoretical permissible nutrient load. We also found that summer HABs in the lake are formed by the simultaneous blooming of diatom Aulacoseira sp. and a microcystin-producing Dolichospermum sp. Both blooming species have mechanisms for the accumulation of nutrients and entry into dormant stages, allowing long-term survival as well as dispersion into neighbouring bodies of water. Both organisms are able to form symbiotic relationships, likely laying beneath their coordinated succession. Simultaneous blooming of various groups of photosynthetic microorganisms capable of forming dormant stages in an initially oligotrophic arctic lake indicate that its ecosystem has dramatically transformed, and in the future, we should expect an increase in HAB magnitude in Lake Imandra.