Introduction

The biological diversity across Earth is characterized by a heterogeneous spatial distribution, which is shaped by the complex interplay of ecological and evolutionary processes. In terms of biodiversity at species level, ancient lakes are considered as diversity “hot-spots” (Rossiter & Kawanabe, 2000). Their prolonged isolation from other ecosystems resulted in extraordinary biodiversity and high numbers of endemic lineages (Brooks, 1950; Martens, 1997; Cristescu et al., 2010). Revealing the geographical and environmental factors that influence the evolutionary and assembly processes that regulate this unique diversity is of a recent fundamental research interest. Moreover, global change during the 21st century and beyond is highlighted as a potential cause for increased extinction risks for a broad variety of organisms in many ecosystems, including ancient lakes (IPCC, 2014). Several studies already provide evidence for the negative effects of global change and anthropogenic pressure over the biodiversity in ancient lakes Tanganyika (Cocquyt, 2000), Ohrid (Matzinger et al., 2007) and Baikal (Moore et al., 2009).

Lake Ohrid, located at the Balkan Peninsula on the border between Macedonia and Albania, is the oldest lake in Europe and well recognized for its unique biodiversity. Recent age assessments place the onset of the lake formation between 1.9 and 1.3 million years before present (Lindhorst et al., 2015; Wagner et al., 2017). Biological investigations available so far describe more than 300 endemic eukaryotic species in Lake Ohrid (Föller et al., 2015). Taking its small surface area of 358 km2 into account, Lake Ohrid may therefore have the highest endemicity index of all lakes worldwide (Albrecht & Wilke, 2008). The question, which mechanisms drive this extraordinary biological diversity, endemism and community structuring, has attracted many scientists ever since the landmark publication of Stanković (1960). Earlier hypotheses, mainly based on observational evidence, emphasized the presence of bathymetric and horizontal barriers in the lake as primary driver for the geographic modes of speciation (Stanković, 1960; Radoman, 1985; Albrecht & Wilke, 2008). More recent studies applying speciation models on the gastropods also reveal the importance of ecological speciation due to the presence of physical barriers in the lake (Schreiber et al., 2011; Hauffe et al., 2016).

While much effort has been devoted to understand the drivers of present-day diversity and community patterns of the main invertebrate groups in Lake Ohrid, relevant studies addressing these questions in the contemporary algal groups are still missing. The intensive ecological studies, arising from the increasing concerns about habitat modification and human impact on Lake Ohrid, are mostly focused on the overall planktonic communities, their seasonal succession and response to increased nutrient input in the lake (Kozarov, 1960, 1974; Allen & Ocevski, 1976; Mitić, 1985, 1992). Diatoms (Bacillariophyta), single-celled siliceous protists, as the most species rich algal group in the lake have the potential for testing the role of the abiotic and biotic factors in structuring communities over spatial scales. The groups’ species richness, endemicity and relictness are well studied and available through the intensive floristic and taxonomic diatom research of the lake (Fott, 1933, 1935; Hustedt, 1945; Jurilj, 1948, 1954; Levkov et al., 2007a; Levkov, 2009; Levkov & Williams, 2011, 2012; Cvetkoska et al., 2012, 2014; Jovanovska et al., 2013; Pavlov & Levkov, 2013). However, the environmental determinants of diversity and community patterns in planktonic and benthic diatoms still remain largely unknown. Hence, we here aim to determine which environmental factors influence the present-day spatial structure of the diatom communities in Lake Ohrid. We specifically targeted the objectives of this study to (i) identify potential spatial subdivision of the diatom communities; (ii) recognize distributional trends of the endemic diversity; (iii) infer the relative quantitative importance of water depth and water chemistry in structuring spatial community patterns; and (iv) provide implications of the observed patterns and their environmental drivers for the diatom ecology and paleoecology in Lake Ohrid.

Methods

Sampling site

Lake Ohrid (41°01′N, 20°43′E, 693 m a.s.l., Figure 1) is calcium bicarbonate dominated and characterized by a mean water depth of 155 m up to maximum 293 m (Lindhorst et al., 2015). Its hydrological balance is regulated via 37.3 m3 s−1 inflow of which ~ 25% origins from direct precipitation, ~ 25% from river input and ~ 50% from the karst aquifers. Outflow is regulated through the river Crni Drim (~ 60%) and evaporation (~ 40%; Matzinger et al., 2006). The average surrounding air temperature is 26°C during summer, and − 1°C during winter, while the average annual precipitation is ~ 750 mm year−1. The lake is oligotrophic and phosphorus limited (Allen & Ocevski, 1977); however, recent studies showed that the mean total phosphorus concentration has risen to 4.5 mg P m−3 as result of the global warming and local anthropogenic impact (Matzinger et al., 2006, 2007).

Fig. 1
figure 1

Location of Lake Ohrid at the Balkan Peninsula (right corner). Topographic and bathymetric map of Lake Ohrid showing the sampling locations at each transect: Kalishta, Struga, Sateska, Ohrid, Peshtani and St. Naum

Full size image

Sampling and analytical work

Sampling was conducted in November 2009 in six shoreline transects (Fig. 1) at 5–10 m distance intervals, starting from 0.5 to 50 m water depth. The transect names are referring to the nearest landmarks or settlements (Table 1). The depth of the lake at every sampling point was measured with a mini echo-sounder device SM-5 Depthmate Portable Sounder®. Two types of samples were collected at each point: a diatom sample, taken from the lake sediment surface, and a water sample from ~ 50 cm above the sediment surface, used for physical and chemical analyses.

Table 1 List of sampling locations at Lake Ohrid and explanatory variables measurements used for the ordination analysis
Full size table

The diatom surface sediments were taken with a Van Veen Grab Sampler. Each sample was acid-cleaned with potassium permanganate (KMnO4) and hydrochloric acid (HCl) to oxidize organics and remove the carbonates (Cvetkoska et al., 2014). Permanent microscope slides were prepared using Naprax® and analysed under oil immersion at 1500x magnification with Nikon Eclipse E-80 microscope, equipped with Nikon Coolpix P600 digital camera. To analyse the diatom communities, at least 400 valves were counted per slide, using the standard transect-based method (Battarbee, 1986). The only exception was Sa 1 (Sateska, 0.5 m water depth) where only 100 valves were counted due to the low diatom concentrations. Diatom taxa were identified using standard literature (Krammer & Lange-Bertalot, 1986, 1997, 2000, 2004) and the dedicated Ohrid and Prespa taxonomic work of Levkov et al. (2007a, b), Levkov (2009), Jovanovska et al. (2013) and Cvetkoska et al. (2014). All diatom samples and microscope slides are stored at the Macedonian National Diatom Collection (MKNDC) at the Faculty of Natural Sciences in Skopje, Macedonia.

Water temperature (T, °C), electric conductivity (Cond.) and pH were measured in situ with Hanna® HI98129 Combo tester. Additional physical and chemical characteristics (Table 1) of the bottom water layer were measured in laboratory on samples collected with a 1L Ruttner sampler. The underwater Photosynthetic Photon Flux Density (PPFD) as a measure of the Photosynthetically Active Radiation (PAR, 400–700 nm) was measured with a LI-192 cosine corrected Underwater Quantum Sensor connected to a LI-250A Light Meter. The diffuse attenuation coefficients (Kd) of the downward PAR were calculated following Zhang et al. (2006). The samples for the analyses of the water chemistry were also collected with 1L Ruttner sampler. The concentrations of nitrate nitrogen (NO3-N), total phosphorus (TP), sulphate (SO42−), silica (SiO2), cobalt (Co), aluminium (Al), barium (Ba), iron (Fe), magnesium (Mg), calcium (Ca), potassium (K) and sodium (Na) were measured according to APHA (1998). The concentration of ammonia nitrogen (NH3-N) was determined by the Spectrophotometric method in accordance with NEN 6472 (1983).

Statistical analyses

The Shannon Index (H) was calculated for each sampling point, as a measure for the α diversity taking the species richness and species’ relative abundances into account. Spatial variations in community structure among the different transects were explored by calculating the Whittaker Index of β diversity (Whittaker, 1960, 1972; Anderson et al., 2011). To test for possible subdivision of the diatom communities along the sampling locations, we applied the unweighted pair group method with arithmetic mean (UPGMA) clustering based on Bray–Curtis similarity measure. In the next step, widespread taxa were removed from the analyses and Bray–Curtis similarity was calculated with endemic taxa percentages only, to assess potential subdivisions among the endemic diversity over the sampled area. Criteria and species classification into widespread and endemic followed Levkov & Williams (2012).

To determine which individual species contributed most to the differences between the samples, we used the species contribution to similarity method (SIMPER), which measures the percentage contribution of each species to average dissimilarity between groups (Clarke & Ainsworth, 1993).

The combined effect of the different environmental parameters on the variation in species composition among the sites was modelled with redundancy analyses, RDA (Rao, 1964; Legendre & Legendre, 1998; Legendre & Anderson, 1999). The linear response in the species dataset, as suggested by the low gradient lengths (Standard deviation = 2.1), justifies the selection of RDA as the most appropriate technique. Diatom abundances were square-root transformed and centred, as this transformation tries rather to stabilize variances, and down-weight rare taxa to avoid overemphasizing their effects over the results (Legendre & Legendre, 1998). The validity of the model was tested with linear transects permutation test with 999 permutations. The explanatory power of each environmental predictor was determined with an interactive forward selection and optimal fit was achieved after 14 predictors were taken into account (Table 3).

The conditional effects of the explanatory variables were further tested with variation partitioning (Borcard et al., 1992; Šmilauer & Lepš, 2014) between two groups of explanatory variables, water chemistry (TP, NH3-N and Cond.) and water depth. Groups were predetermined based on the forward selection results, only variables with P < 0.01 were chosen for the analysis.

To analyse the response of the diatom α diversity (Shannon’s diversity index) to abiotic factors, we implemented a Generalized Linear Model (GLM) ANCOVA model using transect and sediment type as categorical predictors, and environmental variables (Table 1) as continuous predictors. To avoid collinearity, variables with VIF values above 10 were excluded. Environmental variables were selected based on a ‘best subsets’ routine. To compare the performance of the different models generated, we used Akaike’s Information Criterion (AIC). Computations were performed with Statistica V. 10.0 (StatSoft, 2011).

Redundancy analyses and Var-Part were performed in Canoco v.5.0 (ter Braak & Šmilauer, 2012), while the Bray–Curtis similarity analyses, the UPGMA clustering, the diversity indices and SIMPER analysis were implemented in the statistical package PAST (Hammer et al., 2001).

Results

Diatom diversity and community patterns

A total of 290 diatom taxa including morphological varieties were identified in the samples. Dominant in most of the samples are Amphora indistincta Levkov, Pantocsekiella ocellata complex (Pantocsek) Kiss & Ács, Pantocsekiella minuscula (Jurilj) Kiss & Ács, Encyonopsis microcephala (Grunow) Krammer, Pseudostaurosira brevistriata (Grunow) Williams & Round and Staurosirella pinnata (Ehrenberg) Williams & Round. In terms of endemic species richness, ~ 33% taxa are endemic, as compared to the total number of identified taxa.

Lowest AIC values were obtained for a GLM model using an inverse Gaussian error distribution and a power (squared) link function. The model incorporates only four variables with no significant interactions: depth (P < 0.001), sulphates (P  ≤ 0.05), transect (P  ≤ 0.01) and sediment type (P  ≤ .01) (Fig. 2a–d). There is a clear, monotonic tendency of H to decrease with increasing depth, irrespective of other factors. This structural parameter varied also largely depending on transect sampled and sediment typology but did not depend significantly on any other limnological or physiographical variable. The variability of this metric, as measured in terms of standard errors (Fig. 2c–d), shows also a great dependence on transect location and sediment type, the largest variances being observed in Sateska and on mud/sand samples, respectively.

Fig. 2
figure 2

A–D Shannon’s index plotted against the environmental predictors selected by the GLM ANCOVA model

Full size image

The α diversity shows overall minimum and maximum values in the transect at city of Ohrid (code: Oh), Hmin = 1.82, at Oh 10 and Hmax = 3.44 at Oh 3. Generally, the Struga (code: St) transect has the highest average diversity H = 3.13, while the Sateska transect (code: Sa) displays the lowest average of H = 2.64. However, an overall trend of bipartite variation in diversity is notable at each transect (Figs. 2, 3a) for its littoral [0–10(15) m] and sub-littoral zone [10(15)–50 m]. In the littoral zone, the diversity increases towards ca. 5 m water depth, and subsequently decreases to its minima around 10(15) m water depth. Similarly, in the sub-littoral zone, diversity values increase between 10(15) and 25(30) m and successively decline after 30 m water depth. The β diversity values (Fig. 3b) show lowest values (1.67) at the St transect, whereas Kalishta (code: Ka) displays the highest value (2.3) of all transects.

Fig. 3
figure 3

Bathymetric map of Lake Ohrid displaying the α (A) and β (B) diversity and their values

Full size image

The Bray–Curtis similarity analysis of the widespread and endemic taxa shows diatom assemblages cluster by depth as well as transect location, thus supporting delineation of two major groups, littoral and sub-littoral, at a similarity threshold value, S = 0.45. Each group consists of two sub-groups, at S = 0.48 and 0.55, respectively (Fig. 4a), and hence we here define four diatom assemblage zones in Lake Ohrid, excluding the sample from Sa (0.5 m) as most distinct:

Fig. 4
figure 4

a Bray–Curtis similarity analysis of the diatom communities in Lake Ohrid. b Bray–Curtis similarity analysis of the endemic diatom communities

Full size image
  1. (i)

    north-eastern (NE) and eastern (E) littoral, including samples from 0.5 to 15 m water depth collected from the St, Sa, Oh, Peshtani (code: Pe) and St. Naum (code: sN) transects;

  2. (ii)

    north-western to eastern (NW-E) littoral (samples between 0.5–10 m at Ka; at 5 and 15 m from St; and 1–5 m from Sa);

  3. (iii)

    NW sub-littoral (10–50 m, Ka, St and Sa); and.

  4. (iv)

    NE sub-littoral (10–50 m, Sa, Oh, Pe and sN).

The Bray–Curtis similarity applied to the endemic communities (Fig. 4b) delineates five zones at S = 0.2: (i) samples mainly from the lake littoral (1–15 m); (ii) samples from the sub-littoral (15–50 m); (iii) samples from Pe 0.5 m, sN 0.5 m, sN 1 m; (iv) Sa 0.5 m, Oh 0.5 m; and (v) samples from the NE littoral of Ka, St and Oh (Ka 10 m, St 5 m, Oh 10 m, St 10 m, Ka 5 m) which has lowest similarity index as compared to the other locations. The number of endemic taxa per water depth for each transect is shown in Fig. 5.

Fig. 5
figure 5

Endemic richness at each sampling location. The sample collected at 0.5 m water depth at the Sa transect has been removed from the diagram due to possible bias resulting from the low diatom counts

Full size image

The SIMPER analysis shows that ~ 70% of the cumulative dissimilarity between the Bray–Curtis groups is owed to 13 taxa (Table 2). Among them, Amphora pediculus (Kützing) Grunow, Staurosirella sp. and P. brevistriata dominate the littoral zone, while the P. ocellata complex and P. minuscula dominate the sub-littoral, all together already contributing with ~ 50% to the cumulative group dissimilarity.

Table 2 Results from the SIMPER analysis showing the average dissimilarity between the groups as delineated by the Bray–Curtis similarity analysis, percentage contribution (Contrib. %) of the species to the average community dissimilarity, cumulative percentage (Cum. %) and the mean abundance of taxa in the groups (NE, E littoral, NW-E littoral, NW sub-littoral and NE sublittoral). Only taxa with > 1% contribution are shown
Full size table

Environmental variables

The environmental parameters used in the ordination, their measurements and marked minimum and maximum values are presented in Table 1. The four RDA Axes together explain 52% of the total variation in the diatom data, of which 26.1% are explained by Axis 1 and 15.3% by Axis 2. The explanatory variables account for 71.6% of the explained variation in the diatom data. The inflation factors of all variables were lower than the critical heuristic value of 10, indicating that they are not correlated. The only exceptions were downward PAR intensity at the bottom (Edbottom), which showed high positive correlation of r = 0.9 to PAR, and dissolved oxygen (DO) with a r = 0.7 to T (°C), and these two parameters were removed from the analysis. The interactive forward selection by using unrestricted permutation test showed that only four variables have significant (P < 0.01) explanatory power, together accounting for 43% of the total explained variation in the diatom data (Table 3). The RDA results indicate that large proportion of variation in the diatom data is driven by the water depth.

Table 3 Results from the interactive forward selection of the explanatory variables used in the RDA ordination
Full size table

The NE littoral diatom zone appears related to a combination of several water chemistry parameters, NH3-N, PAR, T (°C), pH and Cond., while the variation in the NW littoral zone can be mostly explained by Cond., NO3-N and TP. Samples from the NW sub-littoral align along the TP and Magnesium (Mg) gradient, while water depth and sediment type appear as major drivers of the variation observed in the NE sub-littoral diatom zone (Fig. 6a).

Fig. 6
figure 6

a Axis 1 versus Axis 2 sample scores from the RDA analysis. The symbols conform to the Bray–Curtis similarity groups: NE and E littoral (red triangles), NW-E littoral (yellow squares), NW sub-littoral (pink triangles) and NE sub-littoral (green circles); and b Axis 1 versus Axis 2 species scores from the RDA analysis

Full size image

Variation in the group of centric species, such as Pantocsekiella ocellata (5–10 µm, 3 ocelli), P. minuscula and Cyclotella fottii Hustedt appear to be driven by the depth and sediment type at each sampling location. However, the P. ocellata morphodemes with more than three ocelli and above 10 microns seem to rather be related to the TP and Mg concentrations. The influence of the explanatory variables on the diatom taxa is visualized in Fig. 6b and conforms the SIMPER analyses (Table 2). Species abundances in the NE littoral zone, like Staurosirella spp., Amphora pediculus, Karayevia clevei (Grunow) Round and Placoneis pseudanglica (Lange-Bertalot) Cox are mostly driven by NH3-N, PAR, T and pH, while the species abundances in NW littoral and sub-littoral zones, like Staurosirella pinnata, A. indistincta, Navicula cryptotenella Lange-Bertalot, P. ocellata (> 10 µm, 5 ocelli), P. ocellata (> 10 µm, 3 ocelli), P. ocellata (5–10 µm, > 3 ocelli) are explained by conductivity, NO3-N and TP. The species from the NE sub-littoral zone are correlated to water depth and sediment type (Fig. 6b). The total variance partitioning shows that the environmental and the spatial predictors explain similar proportion of the variance, while the shared fraction is very small (Fig. 7).

Fig. 7
figure 7

Results from the variation partitioning between the explanatory variables with P < 0.01

Full size image

Discussion

The results show distinct spatial patterns of the diatom diversity and community structure within the Macedonian side of Lake Ohrid, driven by a combination of water depth and water chemistry.

Diatom diversity and community patterns

Diatom species richness in Lake Ohrid was estimated to 789 taxa, 117 of which were considered endemic for the lake (Levkov & Williams, 2012). These numbers have certainly increased due to the revisions of several diatom genera and new species descriptions, resulting from the intensive taxonomic and palaeoecological research on the lake over the past few years. In this study, we identified 210 widespread and 80 endemic diatom taxa.

The results yielded distinct vertical and horizontal patterns in diatom diversity and community structure reflected in (i) changes in relative abundances of the widespread taxa along the depth gradients; (ii) presence of different widespread taxa in the NW and NE-E part; and (iii) presence of different rare and endemic taxa at different depths and sampling locations.

The presence of vertical and horizontal patterns of distribution and speciation due to differences in substrate, sedimentation, vegetation, physical and chemical characteristics, was already recognized within the invertebrate groups in Lake Ohrid (Stanković, 1960; Radoman,1985; Albrecht & Wilke, 2008: see Figs. 9, 10). Our results indicate that vertical (bathymetrical) and horizontal zonation can be applied to the diatom communities in Lake Ohrid. Bathymetrically, there are four sub-zones (a) upper littoral (0.5–5 m); local variations in diversity, but generally increasing up to highest values around 5 m water depth; (b) lower littoral [5–10(15) m, or the “Chara belt”]; markedly sharp decline in diatom diversity; (c) upper sub-littoral [10(15)–30 m, the “Dreissena belt”]; increase in the diatom diversity; and d) lower sub-littoral (30–50 m); successive decline in diversity to lowest values at the transition towards the profundal zone. The horizontal zone can be divided in three sub-zones (a) NW, Kalishta zone; the most distinct community, supported by the highest β diversity; (b) NE, Struga and Sateska zone; β diversity reflects high community turnover between each site; and (c) E, Ohrid, Peshtani, and St. Naum zone, less pronounced differences between the community composition at sampling site.

The vertical and horizontal patterns in diatom diversity and community structure closely resemble the spatial variability in sediment structure in the lake (e.g. Vogel et al., 2010). The sediment structure varies along the whole basin, due to the complex geodynamic setting of the lake. The sediments from the western and eastern parts differ in the physical and chemical properties due to the presence of Mesozoic (Apulian) and Paleozoic rocks (Pelagonian), respectively (Robertson, 2004). Such variability may have important influence on water chemistry and consequently, microhabitat structure in the lake. The observed horizontal patterns in diatom diversity and community structure imply a presence of different habitats along the NW-E gradient indeed, probably associated to the complex sedimentary characteristics in the basin. Similarly, to Lake Ohrid, the benthic diatom community from the neighbouring Lake Prespa showed that different species dominated in the littoral zone of the eastern and western parts of the lake, related to differences in geology and sediment structure (Levkov et al., 2007b).

Concerning the endemic diatom taxa, earlier studies (Levkov et al., 2007a; Levkov & Williams, 2012) noticed several spatial levels of their distribution in the whole watershed area, including Lake Prespa where taxa are (i) restricted to only one of the lakes, (ii) common to both lakes and iii) restricted to the spring area and/or Lake Ohrid. Our results show that some endemic taxa, like Amphora ohridana Levkov, Caloneis acuta Levkov, Sellaphora krsticii Levkov, Nakov & Metzeltin appeared to have limited distribution in the littoral zone, while others, such as Caloneis biconstrictoides Levkov, Cymbopleura juriljii Levkov & Metzeltin and Sellaphora macedonica Levkov & Metzeltin were restricted to the sublittoral zone of the lake. The comparison of the number of endemic taxa per water depth (Fig. 5) yielded that in general the endemicity is higher in the sublittoral zones, similarly as observed for the gastropod communities (Hauffe et al., 2011). An exception is the shallow littoral zone (0.5–1 m depth) of Peshtani and St. Naum, where the endemic taxa comprised up to 25% of their diatom community. Overall, the results from this study covering Lake Ohrid only, imply that the presence of horizontal and vertical barriers might act as factors limiting the distribution of the endemic species in the lake.

Relative quantitative importance of environmental variables

The results from the RDA analysis indicated that water depth and water chemistry almost equally influence the spatial diatom community dynamics in Lake Ohrid.

Diatom diversity and community structure in the lake changed along the water depth gradients as different species dominate the littoral and the sublittoral zones. The diatom communities in the sandy littoral area, 0.5–1 m depth, are dominated by Achnanthidium sp., Amphora indistincta, A. pediculus, Cocconeis placentula Ehrenberg, E. microcephala, Gomphonema sp., Navicula antonii Lange-Bertalot, Placoneis balcanica (Hustedt) Lange-Bertalot, Metzeltin & Levkov, Pseudostaurosira sp. and S. pinnata. Since the nearshore area is characterized as the most dynamic part of lakes, it is expected to reflect the presence of distinctive habitats, associated with specific biogeochemical and/or metabolic processes, hence species communities (Lewis, 2009). The shallow water diatom communities in Lake Ohrid represent a mixture of “pioneering” taxa, which prefer sandy or macrophytic substrates (Levkov et al., 2007a). These communities appeared associated with a combination of multiple environmental factors, among others the concentration of ions, light availability, temperature and pH. The influence of the water depth over diatom community structure in the upper littoral and the sublittoral zone was evident from the higher abundance of planktonic diatom taxa such as the Pantocsekiella ocellata complex. At greater water depths, the recourses for benthic species become limited (e.g. light availability, temperature), to which species often respond with decreased abundances. In contrast, these resources remain relatively stable in the epilimnetic zone, possibly promoting high abundances of planktonic species (Passy, 2008; Cantonati & Lowe, 2014). The strong influence of water depth in regulating the spatial diatom community patterns has also been shown in other lakes studies (Yang & Duthie, 1995; Stone & Fritz, 2006; Laird et al., 2010; Wang et al., 2012).

The decline of diatom diversity in the lower littoral part, however, can be associated with the substrate structure rather than any specific change of the environmental factors. The lower littoral part of Lake Ohrid is characterized by dense growth of the macrophytic algae from the genus Chara Linnaeus, forming submerged meadows extending between 5 and 10(15) m water depth (Trajanovska et al., 2014). The low diversity and communities dominated by epiphytic diatoms, like C. placentula and E. microcephala, imply that the “Chara belt” might act as dispersal barrier for the diatom communities in the lake. The importance of the Chara belt as habitat for some species and as barrier for others has been also stressed by Radoman (1985).

The group of environmental variables indicated that diatom distribution in Lake Ohrid is also driven by water chemistry (Fig. 6), particularly TP, NH3-N and conductivity. The influence of TP and NH3-N over the diatom communities in the area of Kalishta, Struga and Sateska reflects ongoing eutrophication and wastewater pollution in the northern part of the lake. Lake Ohrid is a typical oligotrophic lake, however, increasing pollution and nutrient loads from the surrounding settlements and the river tributaries, Sateska and Koselska, were already recognized as the major source leading to eutrophication of the northern lake area (Matzinger et al., 2007; Veljanoska-Sarafiloska et al., 2011; Trajanovska et al., 2014). Start of the human impact and eutrophication at  ~ 1955 AD was inffered from the analysis of several shallow gravity cores recovered from the Kalishta spring area where changes in coastal sedimentation rate and increase in organic matter content was determined (Matter et al., 2010).

The TP concentrations that we measured in the northern lake area were up to ten times higher than those measured in the eastern area (Table 1), and the diatom communities were accordingly dominated by nutrient tolerant (e.g. Kelly et al., 2008) and/or meso-eutrophic species (e.g. Van Dam et al., 1994), such as E. microcephala, N. cryptotenella, Nitzschia linearis Smith and P. ocellata (> 10 µm).

Implication for diatom ecology and palaeoecology in Lake Ohrid

Several studies showed that Lake Ohrid is under increased human pressure which poses a threat to its unique biota (Matzinger et al., 2007; Albrecht & Wilke, 2008), and our results further corroborate these findings. The significant contribution of total phosphorus and nitrogen ammonia concentrations in explaining the diatom distribution in the northern lake area should raise serious concerns about the effects of the ongoing eutrophication and pollution over the lake diversity.

In the context of the Water Framework Directive ‘WFD’ guidelines, diatoms have been suggested as the most appropriate bioindicators to be used in regular monitoring programs of the freshwater bodies in Macedonia (Krstić et al., 2007). Although the information diatoms provide is rapid, reliable and cost efficient, the biological monitoring of Lake Ohrid still largely relies on macrophytes, invertebrates and fish, while diatoms were included only in short survey studies (e.g. Schneider et al., 2014). Despite the extensive floristic, taxonomic and palaeoecological diatom research, a systematic diatom-based ecological study on Lake Ohrid has not yet been performed. Our study provides initial data for interactions between environment and spatial distribution of diatom communities and thereby a baseline for establishing a more accurate biomonitoring strategy for the lake. Systematic diatom sampling and interpretation need to take the non-homogenous diatom distribution within the lake into account, hence a sampling strategy capturing the communities across the vertical and horizontal zones.

Contrary to the diatom neo-ecology, Lake Ohrid has been in the spotlight of palaeoecological research after its inclusion as a target area of the International Continental Scientific Drilling Program, and the start of the interdisciplinary research project Scientific Collaboration on Past Speciation Conditions in Lake Ohrid, SCOPSCO in 2009. Several sediment sequences, covering the last glacial–interglacial period, were recovered within the project, prior an ICDP deep drilling campaign took place in 2013 (Wagner et al., 2014). A ~ 600-m-long master sediment sequence was retrieved from the deepest part of the lake. Being one of the most abundant microfossils in the sedimentary archive of the lake, diatoms were intensively studied as one of the most important biological proxies which can provide information about the past in-lake paleo-processes and the past environmental and climate conditions in the area (Reed et al., 2010; Cvetkoska et al., 2012, 2016; Jovanovska et al., 2016; Zhang et al., 2016).

However, in the lack of systematic diatom-based ecological data all these studies were limited to qualitative palaeoenvironmental interpretations. The high degree of endemism and the lack of modern analogue species, which represent a common challenge for Quaternary palaeoecological interpretations of fossil diatom data from other ancient lakes (Mackay et al., 2010, Wilke et al., 2016), are no exception in Lake Ohrid. Furthermore, while all these studies were performed at high-resolution taxonomic level, an extraordinary phenotypic diversity was observed in some of the species complexes e.g. Pantocsekiella ocellata (Cvetkoska et al., 2012) but, the drivers behind it remain uncertain. The data presented here demonstrate that populations with different morphological traits are driven by different environmental factors (e.g. P. ocellata morphs), corroborating the possibility for environmentally driven phenotypic diversity. While the question how persistent and genetically different these phenotypes are remains open for future studies, their presence in the fossil record renders important palaeoecological information. Considering the here-provided information on the environmental drivers of diatom community and diversity patterns in Lake Ohrid, future studies should focus on constructing a reliable ‘within-lake’ diatom calibration dataset for quantitative palaeoecological reconstructions. Overall, a systematic assessment integrating ecological and palaeoecological aspects is necessary to determine boundary conditions and to improve the action plans towards protection and conservation of Lake Ohrid with its exceptional biodiversity.

Conclusions

In this study, we demonstrated distinct spatial patterns of community and diversity in diatoms from Lake Ohrid. The widespread and endemic species have limited distribution across the northern and eastern part of the lake, which appear to be related to the presence of vertical and horizontal barriers. The multivariate analyses reveal that the distribution of diatoms is related to variation in water depth and water chemistry gradients. Specifically, both variables appear to equally influence the diversity and the extant community structure. The significant influence of total phosphorus concentration on the diatoms in the northern parts of the lake reflects the ongoing eutrophication as a result of the increased anthropogenic nutrient loads.

While this study provides initial insight into the diatom distribution from the Macedonian side of Lake Ohrid, future studies should expand the spatial scale by including the whole lake and its spring area. Moreover, integrating biotic and abiotic variables in spatial and temporal dimension should be considered in order to unravel the assembly processes driving extant diversity patterns in Lake Ohrid.