Ecosystem regimes and responses in a coupled ancient lake system from MIS 5b to present: the diatom record of lakes Ohrid and Prespa

. We reconstruct the aquatic ecosystem interactions since the last interglacial period in the oldest, most diverse, hydrologically connected European lake system, by using palaeolimnological diatom and selected geochemistry data from Lake Ohrid “DEEP site” core and equivalent data from Lake Prespa core, Co1215. Driven by climate forcing, the lakes experienced two adaptive cycles during the last 92 ka: “interglacial and interstadial” and “glacial” cycle. The short-term ecosystems reorganizations, e.g. regime shifts within these cycles substantially differ


Introduction
The structure and functions of the continental ecosystems constantly changed throughout the Quaternary, driven by a series of successive climate variations at different temporal scales.As dynamic entities, the ecosystems developed through these fluctuations with a "panarchy", defined as a series of hierarchical adaptive cycles interacting across multiple scales of space and time (Holling, 2001;Delcourt and Delcourt, 2004;Allen et al., 2014).Four functional phases, controlled by top-down and bottom-up processes, have been identified within these cycles, growth (r), conservation (k), release ( ), and reorganization (α).Change in the processes controlling the ecosystem's structure and functions can lead to regime change and creation of new selforganized structures (Holling, 1986).Regime shifts may occur within one adaptive cycle or between adaptive cycles at different scales, and cause reorganization or even sometimes a complete ecosystem collapse (Allen et al., 2014).Such catastrophic regime changes, marked by species communi-Published by Copernicus Publications on behalf of the European Geosciences Union.
ties loss and consequently, collapse of the ecosystem's functions have been recorded in many of the world's oldest lakes, like Baikal in Russia, Hövsgöl in Mongolia and Malawi in East Africa (Karabanov et al., 2004;Cohen et al., 2007;Burnet et al., 2011;Scholz et al., 2011).Based on the intensity of forcing, and their physical, chemical and biological properties these lakes had different capacities to adapt, restore and/or reorganize their panarchy.The adaptive capacity of an ecosystem is intrinsically related to its resilience (Allen et al., 2005;Allen and Holling, 2008), and therefore, research on sediment records from ancient lakes is important for understanding the past adaptive capacities, regime changes, and resilience of these unique ecosystems.
In Europe, such Quaternary archives are neighbouring lakes Ohrid and Prespa, located on the Balkan Peninsula.Since the Pliocene origin of their tectonic basins, more than 1 Ma co-existence and co-evolution (Stanković, 1960) resulted in the oldest and most diverse permanent lake systems in Europe (Albrecht and Wilke, 2008;Levkov and Williams, 2012).Understanding the external and internal drivers for these basins is critical for determining the high level of biological endemism of especially Lake Ohrid.At present, the lake system has unique spatial interactions due to hydrological connectivity via a karstic system in Mt.Galičica (Anovski et al., 1997(Anovski et al., , 2001)).About 20 % of the Lake Ohrid water inflow originates from Lake Prespa, and, therefore, the much shallower and more nutrient rich Lake Prespa can potentially drive changes in Lake Ohrid (Matzinger et al., 2006).By applying a system analytical approach, Matzinger et al. (2006) assessed the anthropogenic influence on Lake Prespa, quantified the underground water and nutrient flow and developed a linear phosphorus (P) model.By using the 1960s' water-level high-stand of Lake Prespa (∼ 855 m a.s.l.) as a baseline for the calculation, the model indicates that a 20 m water level decrease of Lake Prespa can cause a five-fold increase of its P concentration, which may lead to a 30 % increase of the P load from Lake Prespa to Lake Ohrid (Matzinger et al., 2006).Such a scenario can potentially trigger a regime shift in the physico-chemical properties of Lake Ohrid and affect its fragile endemic community (Watzin et al., 2002;Matzinger et al., 2006Matzinger et al., , 2007)).Therefore, parallel analysis of long-term temporal variations under changing climate boundary conditions from Quaternary sedimentary records of both lakes is important for identifying the potential influence of Lake Prespa on Lake Ohrid.
As part of the project Scientific Collaboration on Past Speciation Conditions in Lake Ohrid (SCOPSCO), a successful International Continental Scientific Drilling Program (ICDP) drilling campaign was carried out at Lake Ohrid in spring 2013.A sediment sequence spanning > 1.2 Ma (Wagner et al., 2014a) was recovered at the "DEEP site" in the central part of the lake.Multi-proxy analyses of this sequence currently proceed towards understanding the influence of the past geological and environmental events on the biological evolution of the lake taxa, as one of the main project objectives.
Aiming to assess the potential interactions between the lakes, we compare Lake Prespa's dynamics, as already inferred from core Co1215 diatom, pollen, and geochemistry data with the diatom and selected bio(geochemical) data from the same age interval (ca.92.0 ka to present) of the DEEP site sediment sequence.Here, we apply the concept of panarchy and consider lakes Ohrid and Prespa as lower-scale entities of the higher hierarchical, landscape scale, as justified by their hydrological connection.Trends in the multiproxy data set are used to reconstruct and compare the adaptive cycles, regime shifts and assess the response mechanisms to external disturbances on ecosystem scale (e.g.lakes Ohrid and Prespa).At higher hierarchical level, we aim to assess the potential interconnectivity between the lakes in order to understand the cross-scale linkages and recognize the capacity of Lake Prespa as a driving force of regime shifts in Lake Ohrid.

Environmental setting
Harbouring more than 300 endemic species (Föller et al., 2015), Lake Ohrid, Macedonia/Albania (41 • 01 N, 20 • 43 E, 693 m a.s.l., Fig. 1) has the highest index of endemism of all ancient lakes, when its surface area of 358 km 2 is taken into account (Albrecht and Wilke, 2008).With a mean and maximum water depth of 155 and 293 m, respectively (Lindhorst et al., 2015), it is a deep, calciumbicarbonate Ca(HCO 3 ) 2 dominated, oligotrophic lake.The hydrological balance of the lake is regulated by 37.9 m 3 s −1 inflow, of which ∼ 25 % derives from direct precipitation, ∼ 25 % from river input and ∼ 50 % from the karst aquifers, and the outflow through the river Crni Drim, (∼ 60 %) and evaporation (∼ 40 %, Matzinger et al., 2006).The estimated water residence time is ca.70 years.The average monthly air temperature of the surrounding area is 26 • C during summer and −1 • C in winter and the annual precipitation averages to 750 mm year −1 (Wagner et al., 2009).Recent investigations revealed an increase in total P concentration from ∼ 1.3 mg P m −3 in historic times (ca.150 years ago) to 4.5 mg P m −3 today (Matzinger et al., 2006(Matzinger et al., , 2007)).(Wagner et al., 2009(Wagner et al., , 2010(Wagner et al., , 2012) ) and the DEEP site core 5045-1, retrieved during the SCOP-SCO deep drilling campaign at Lake Ohrid in 2013 (Wagner et al., 2014a).Modified from Francke et al. (2016).Lake Prespa, located at the border of Macedonia, Albania, and Greece (40 • 46 -41 • 00 N, 20 • 54 -21 • 07 E, Fig. 1) is situated at 849 m a.s.l.It is a tectonic lake which belongs to a former lake complex, called Dessarets (Stanković, 1960), but its age is still uncertain.Although it hosts fewer endemic species, it has been suggested that the endemic faunal species might be older than those in Lake Ohrid (Karaman, 1971).In terms of its physical properties, Lake Prespa is a large, 254 km 2 , relatively shallow lake with ∼ 14 m mean and ∼ 48 m maximum water depth.At present, the estimated hydraulic residence time is ca.11 years, but it varies highly with changing lake levels; and for example, at ∼ 19 m mean water depth during the 1980s, this number was ∼ 17 years (Matzinger et al., 2006).The water input, 16.92 m 3 s −1 , is via river inflow and catchment runoff (56 %), direct precipitation (35 %), inflow (9 %) from the nearby Lake Micri Prespa (Fig. 1), and groundwater input (no data available).The output is via evaporation (52 %), water abstraction for irrigation (2 %) and subsurface outflow through the karstic aquifers of Galičica Mountain (46 %, Matzinger et al., 2006).The average monthly temperature in the surrounding area is 21 • C in summer and 1 • C in winter, and the annual precipitation varies between 720 in the lowlands and 1200 mm year −1 on the mountains (Hollis and Stevenson, 1997).The lake is highly sensitive to external distur-bances (e.g.climate change) and water loss; major lake-level fluctuations including a lake-level decline of almost 10 m due to restricted precipitation and increased water abstraction for irrigation have been documented between 1950 and 2009 (Popovska and Bonacci, 2007).Recent monitoring programs of the lake revealed ∼ 31 mg m −3 TP concentration (Matzinger et al., 2006), and ∼ 162 t year −1 TP input from the river tributaries (Krstić, 2012).

Core recovery and chronology
Core "5045-1" was recovered at the DEEP site in the central part of Lake Ohrid within the scope of the SCOPSCO deep drilling campaign in 2013 (Fig. 1).Six parallel boreholes (A-F) at 243 m water depth were drilled down to a maximum sediment depth of 569 m below the lake floor (m blf) and a total 1526 m of sediment cores was recovered.Comprehensive details of the drilling procedure, sediment recovery and results from the preliminary analyses of core catcher samples are provided in Wagner et al. (2014a).The age model for the upper 247.8 m of the composite sequence of the DEEP site for the sediment core and downhole logging data is based on 11 tephrostratigraphic tie points (first order tie points), and on tuning (bio-)geochemical proxy data to orbital parameters (second order tie points) (Baumgarten et al., 2015;Francke et al., 2016).In order to compare the same age intervals in lakes Ohrid and Prespa, we analysed the uppermost 37.5 m of the DEEP site sequence, which represent the last ca.92.0 ka according to the age-depth model.This period is well-constrained by the occurrence of the Y5 tephra layer at ca. 39 ka, the Y3 tephra layer at ca. 29 ka, and the Mercato tephra layer at 8.5 ka in the DEEP site sediment sequence (Leicher et al., 2016).
The age-depth model of the 17.7 m long core Co1215 from Lake Prespa is based on 11 tephrochronological tie points, radiocarbon 14 C dating, electron spin resonance (ESR) dating and comparison with the NGRIP ice core data.An age of ca.92.0 cal ka BP was extrapolated for the base of the sequence; details are provided in Damaschke et al. (2013).

Analytical work
The analysis of the 37.5 m sequence from the DEEP site comprises non-destructive X-ray fluorescence scanning (ITRAX core scanner, Cox Analytical, Sweden) for potassium intensities (K), and measurements of total carbon (TC) and total inorganic carbon (TIC) by a DIMATOC 100 carbon analyzer (Dimatec Corp., Germany).The total organic carbon (TOC) content was calculated from the difference between TC and TIC.Total nitrogen (TN) was analysed using an elemental analyzer (vario MICRO cube, Elementar Corp.) after combustion at 1150 • C and used to calculate the ratio between TOC and TN (abbreviated as C / N).Bio-  (Grimm, 1987) Diatom assemblages zonation CONISS (Grimm, 1987) Multi-proxy stratigraphic diagram C2 (Juggins, 1991(Juggins, -2007) ) Ordination PCA (SD = 2.2), DCA (SD > 2.5), Canoco 5 R "stats" v. 0.8-2 (ter Braak and Šmilauer, 2012) (R Core Team, 2012) genic silica concentrations (BSi) were inferred by means of Fourier Transform Infrared Spectroscopy at the University of Bern, Switzerland following the method outlined in Vogel et al. (2008) and the calibration after Meyer-Jacob et al. (2014).All bio(geochemical) analyses are described in more detail in Francke et al. (2016).
Diatom analysis on the DEEP site core was performed on 235 samples in total (resolution 16 cm, ca.0.3-0.5 ka), prepared from freeze-dried sediment subsamples according to the method of Cvetkoska et al. (2012).Known aliquots of cleaned samples were used to estimate the absolute diatom concentrations, "DC" ([n] valves ×10 7 g −1 dry sediment), employing the sedimentation tray method of Battarbee (1986).Permanent diatom slides were prepared using Naphrax ® as mounting medium.Around 350 diatom valves per slide were counted at 1500× magnification with an Olympus BX51 microscope.As an exception, at least 100 valves per sample were counted in diatom assemblages with low species diversity and/or concentration.Species identification followed Levkov et al. (2007), Levkov andWilliams (2011, 2012), Cvetkoska et al. (2012Cvetkoska et al. ( , 2014b)), Jovanovska et al. (2013), andPavlov et al. (2013).Diatom accumulation rates, "DAR" ([n] valves ×10 6 cm −2 ka −1 ) were calculated by multiplying the dry sediment accumulation rates, "SAR" (g cm −2 ka −1 ) and the DC, Anderson (1989).To assess the preservation quality in the sequence, the F index of Ryves et al. (2001) was calculated for Cyclotella fottii Hustedt as a proportion of the number of pristine valves and the total number of classifiable valves of the species.Values of F = 0 indicate partial dissolution of all valves and F = 1 means ideal preservation.The variation in the di-atom data was explored by detrended correspondence analysis, "DCA", and with a gradient length of 2.2, principle component analysis, "PCA" was appropriate (Jongman et al., 1995).To maximize the apparent variance, PCA was run on species morphotypes at a sub-species level.The Spearman's rho rank-correlation test ("stats" version 0.8-2; R Core Team, 2012) was used to analyse the correlation between the dominant diatom species abundance, the DAR, DC and the (bio)geochemistry data (Table 3).At both sites, the constrained stratigraphic diatom zones were defined by the incremental sum of squares distance metric of the CONISS program (Grimm, 1987).
The analytical approaches of the geochemistry, isotope, lithology, pollen, and diatom analyses of Lake Prespa core Co1215 have been extensively described in Aufgebauer et al. (2012), Wagner et al. (2012), Damaschke et al. (2013), Leng et al. (2013), Cvetkoska et al. (2014aCvetkoska et al. ( , 2015)), and Panagiotopoulos et al. (2014).Diatom analysis on core Co1215 was carried out on 222 samples in total (resolution 8 cm, ca.0.1-0.4ka) and the variation in the diatom data was analysed with DCA based on the > 2.5 gradient length (Cvetkoska et al., 2015).Comparative details of the diatom analytical methods applied to both site cores, DEEP and Co1215, are presented in Table 1.

Principal diatom results and environmental controls
In total, 141 diatom species were identified for the last ca.92.0 ka of the DEEP site core.The diatom diagram (Fig. 2) displays 31 species and morphotypes, present at > 2 % in the assemblages.The diatom assemblages are dominated by planktonic taxa, mainly from the genus Cyclotella (Kützing) Brébisson.The facultative planktonic "FP" and benthic life habitat groups are represented with higher species numbers, but low overall relative proportions due to the large depth of the coring location, which is beyond the optimal depth for their distribution in the lake.Based on previous studies (Stanković, 1960;Allen and Ocevski, 1976), the vertical distribution of the phytoplankton in Lake Ohrid is primarily driven by temperature, in combination with light, nutrient availability, and water mixing.The last 92.0ka of the DEEP site diatom record are generally dominated by three planktonic species with different autecology, C. fottii, C. ocellata and C. minuscula (Jurilj) Cvetkoska.Cyclotella fottii is an oligotrophic, stenothermic species, adapted to low temperatures; viable populations of the species have been found at more than 100 m water depth, where the light availability is far below the optimal intensity for photosynthesis and temperature is permanently reduced (Stanković, 1960).This indicates that low light availability is not limiting C. fottii, while temperature is likely more important factor for its vertical distribution and abundance in the phytoplankton.Cyclotella ocellata is an ecologically successful species of wide distribution, though, in Lake Ohrid, it dominates in the epilimnion during the spring and summer season (Stanković, 1960;Zhang et al., 2016).This is the zone of warmer temperatures, higher light availability, but lower nutrients than the hypolimnion.However, in Lake Ohrid C. ocellata has been interpreted as mesotrophic and an indicator of increased productivity in comparison to C. fottii (Wagner et al., 2009;Lorenschat et al., 2014;Zhang et al., 2016).As a result of the lack of autecological data for C. minuscula we relied on a previous study by Cvetkoska et al. (2014a), the findings of Winder et al. (2009), and the comparison with the geochemical proxies.At the DEEP site, maxima in the relative abundance of C. minuscula correspond to the Y5, Y3, and the Mercato tephra layers and peaks in K intensities.This is also shown by the high-resolution study of the diatom response to the Y5 tephra impact on lakes Ohrid and Prespa by Jovanovska et al. (2016).The tephra influx/deposition increases the silica availability, which is enhancing the diatom productivity, consequently leading to strong competition for nutrients and light resources.
At present, most of the FP and benthic species in Lake Ohrid are found between the shallow littoral zone and ∼ 50 m water depth (Levkov et al., 2007).The overall low relative abundances, 5-15 %, of these species in the DEEP site record result from the coring location, and thus add little supporting evidence to the interpretation.However, the periods of their slightly increased abundances during the last ca.92.0 ka can be interpreted as an indication for intensified transport from the shallower parts of the lake due to wind-induced wave activity and lake-internal currents (Vogel et al., 2010c).
The DC and DAR (Figs. 2,4) show high positive correlation values with the other (bio)geochemistry proxies (BSi, TOC, TIC, C / N, see Table 3).In Lake Ohrid sediments, BSi and TOC are of mainly authigenic origin, indicating the primary productivity in the lake (Wagner at al., 2009;Vogel et al., 2010b;Francke et al., 2016).In addition, by adjusting the DC to the sedimentation rates, DAR can be used to qualitatively assess the overall, glacial-interglacial scale trends of diatom productivity in the record.
Regarding the interpretation of the C / N ratios, as an indication for the source of organic matter "OM", we point to the study of Lacey et al. (2014) on the Lini sequence where it was demonstrated by means of Rock Eval pyrolysis that in Lake Ohrid, the "OM" is predominantly of pure algal source.The contribution of vascular vegetation to the OM at the DEEP site is also highly unlikely due to the coring location.Nonetheless, low C / N ratios (∼ 4) during the glacials can also result from enhanced supply of clay-bound ammonium following mineral soil erosion rather than algal OM, along with ongoing organic carbon degradation (Holtvoeth et al., 2016).
In nature, the biotic response to multiple environmental factors depends on the response to the single dominant (limiting) driver, and the chance of a driver of large effect being present increases with the number of drivers (Brennan and Collins, 2015).On a glacial-interglacial scale, multiple environmental drivers change synchronously with rising and/or declining temperatures, like the light and nutrient availability.In case of Lake Ohrid, the autecology of the species (e.g.alternations between C. fottii and C. ocellata dominance) and the significant correlation of their relative abundances with the other productivity related proxies (TOC, BSi, DC, DAR) allow us to infer the orbital timescale trends in the diatom data as primarily driven by the large-scale temperature changes.Consistent with the interpretation of the diatom response at the Lini site (core Co1262; Zhang et al., 2016), we do not exclude the influence of other drivers, like nutrients, wind-induced water mixing, thermal stratification, light etc., but consider these factors as drivers that may play a more important role during millennial timescale events.

Ecosystems dynamics
In the DEEP site sequence, four major hierarchical diatom zones (DZ 1-4, Fig. 2), each divided into 4-6 subzones, were identified during the last ca.92.0 ka.The diatoms are compared to selected geochemistry data, i.e. the biogenic silica content (BSi) and potassium (K) intensities inferred from XRF core scanning in the record (Fig. 4).Diatom zones are discussed according to the Marine Isotope Stages (MIS) (Lisiecki and Raymo, 2005).For a detailed overview of the diatom data and palaeoenvironmental interpretation of Lake Prespa core Co1215, the reader is referred to Cvetkoska et al. (2014aCvetkoska et al. ( , 2015)).
The biological, chemical and physical components of the ecosystems constantly fluctuate, although sometimes at slow rates (Scheffer and Carpenter, 2003).Therefore, we apply the term regime, referring to the periods of diatom and bio(geochemistry) inferred ecosystems dynamics in terms of trophic state and/or lake level.Changes in the diatom communities' structure related to changes in the inferred lake's dynamics are considered as regime shifts.In case of Lake Prespa, multiple shifts between planktonic and benthic species-dominated assemblages mark the stratigraphic diatom data.The dominance of benthic diatom assemblages reflects the lower lake-level regimes, and the related physical changes in habitat distribution, but does not indicate higher lake productivity.Indeed, shifts in Lake Prespa trophic regimes are marked by changes within the planktonic diatom assemblages, (e.g. from Cyclotella to Aulacoseira dominated assemblages).
Regarding the comparison between both sediment sequences, the age models of Ohrid DEEP site core and Prespa core, Co1215 are independent and only few direct correlation points via the occurrences of well-known tephra layers exist (e.g., the Y5 and Mercato tephra).The dating errors of Ohrid DEEP site core are on the 2σ confidence level for the tephra layers, while each cross-correlation point to orbital parameters include an error of ±2000 years (Francke et al., 2016;Leicher et al., 2016).The dating errors of Lake Prespa core, Co1215 can be estimated to be < 1000 years for the last 40 ka, but increase significantly to several thousands of years in the lower part, where ESR dating, poorly known tephrostratigraphy, and tuning of organic matter content to marine and ice core records provide less well constrained chronological tie points (Damaschke et al., 2013).Therefore, the comparisons of the diatom data (and inferred temperature/trophic changes) need to be done carefully.This is especially the case of the single spikes, which may easily have an offset of several centuries to millennia.The MIS 5b period in Lake Prespa diatom record was interpreted as an unstable, transitional phase of moderate, but variable moisture availability and temperatures.Lowlake-level phases characterized by oligo-mesotrophic conditions in Lake Prespa occurred at 90.0, 88.4, and 85.6 ka (Cvetkoska et al., 2015).
The basal diatom zone, (DZ 4d) in Lake Ohrid DEEP site core is marked by dominance of the C. ocellata complex until 88.7 ka, when C. minuscula peaks with ca.65 % relative abundance, taking over a major part of the epilimnetic diatom productivity to the end of MIS 5b (Figs. 2, 4).The autecology of C. minuscula indicates that its dominance was favoured by the increased clastic input in this part of the record, rather than warmer temperatures, as it is also reflected by the corresponding increase of K intensities and SAR (Fig. 4) (Francke et al., 2016).The low DAR match well with low BSi, TOC and C / N ratio, suggesting less productivity and/or preservation of the organic matter (Francke et al., 2016).The overall comparison between both records indicates that both lakes responded to the similar, transitional phase of moderate temperatures, but lake-level fluctuations in Lake Prespa do not find an equivalent in Lake Ohrid DEEP site.In Lake Prespa, a climatic drought between 85.0 and 83.0 ka causing a lake-level decline and a shift to a regime of highest productivity in the record was interpreted from the diatom and pollen data (Panagiotopoulos et al., 2014;Cvetkoska et al., 2015).A second lake-level low-stand at Lake Prespa between 77.6 and 76.6 ka, related to decreased effective precipitation, occurred before the end of the last interglacial, which is marked at 70.2 ka in the diatom record (Cvetkoska et al., 2015).Sedimentological, hydroacoustic data and bivalves indicate another lake-level low-stand around 74.0 ka (Wagner et al., 2014b), but this is not reflected in the Co1215 diatom data.
In Lake Ohrid DEEP site diatom record, the period between ca.85.0 and 78.0 ka is characterized by the dominance of C. ocellata morphotypes with 4 and 5 ocelli (Fig. 2) and high overall DC and DAR (Figs. 2, 4).Similar diversification and presence of different C. ocellata morphotypes during MIS 5, MIS 3 and the Holocene was observed by Reed et al. ( 2010) and Cvetkoska et al. (2012).To this point, the hypothesis of climatically triggered species evolution seems quite plausible, even though revealing the nature of this diversification requires more data from the earlier Quaternary.The higher productivity inferred from the diatom data in this part of the record is also seen in the increased BSi content, high TOC values and C / N ratio (Fig. 4).In a phase of constant SAR and catchment vegetation dominated by www.biogeosciences.net/13/3147/2016/Biogeosciences, 13, 3147-3162, 2016  mesophilous/montane trees (Sadori et al., 2016), this implies a temperature-related increase in the diatom productivity.The diatom signal of warmer climate is also supported by the high TIC content, indicative of increased endogenic calcite precipitation and/or preservation.A small increase in the abundance of C. fottii at ca. 76.5 ka, correlates with a return to dominance of the typical C. ocellata morphotypes.
The diatom and geochemistry data show that the overall dynamics of the Lake Ohrid DEEP site during MIS 5 are primarily related to the regional trends in climate change, from cold/dry (MIS 5b) to warmer/wetter conditions during MIS 5a (Bar-Matthews et al., 2003;Martrat et al., 2004).The absence of similar abrupt peaks in both records, DEEP and Co1215, and the lack of mesotrophic species and bio(geochemical) evidence for a lake-level decline in Lake Ohrid support the notion that the periods of low lake levels and highest productivity in Prespa between 85.0 and 71.0 ka did not affect Lake Ohrid, at least in its deepest and central parts.

5.2
The last glacial, MIS 4-2 (71.0-14.0ka; DZ 4-1f) 5.2.DZ 4a,3e,part of DZ 3d) In Lake Prespa, the glacial conditions during MIS 4 triggered a progressive opening of the landscape (Panagiotopoulos et al., 2014) and a regime shift to low lake level, evident by low plankton abundances and diatom growth restricted to the shallow littoral habitats in a low productive lake (Cvetkoska et al., 2015).
Contrary, the glacial climate signal in Lake Ohrid DEEP site diatom record appears muted between 70.5 and 65.3 ka due to the high relative proportions, 30-50 %, of the epilimnetic, thermophilous, spring/summer species C. ocellata (Fig. 2, DZ 4a).At the same time, the low DAR, BSi and TOC point to low diatom productivity and/or enhanced OM decomposition.High C. minuscula abundance and K intensities in this period (Fig. 4), indicative of enhanced clastic supply due to catchment erosion, are consistent with climate cooling.Nonetheless, the high relative abundances of C. ocellata between ca.70.5 and 65.3 ka indicate that probably the spring/summer temperatures of the epilimnetic zone were still favourable for its growth.This implies a gradual transition of Lake Ohrid from interglacial to glacial conditions, unlike Lake Prespa, where the MIS 5/4 transition is marked with a sharp boundary at ca. 70.2 ka.
The low proportions of epilimnetic taxa and dominance of C. fottii between 65.3 and 60.0 ka (Fig. 2) in a zone of low preservation quality is indicative of low water temperatures and weak thermal stratification of the lake.In addition, the high proportions and/or dominance of C. fottii morphotypes > 30 µm can be related to ultra(-oligotrophic) conditions, and consequently low competition and decreased grazing pressure, reflecting the cold glacial climate conditions.Similarly, large-celled diatom populations have been noted in other ancient lakes (e.g.Baikal, Hövsgöl) and it is supposed that they reflect oligotrophy, in contrast to smaller-celled diatoms which are more competitive in conditions of limiting nutrients and light (Stoermer et al., 1989;Edlund et al., 2003).

5.2.2
The interstadial, MIS 3 (57.0-29.0ka; part of DZ 3d, DZ 3c-a, DZ 2d, c) In Lake Prespa, the increased moisture availability at the interstadial onset prompted an increasing lake level, but the productivity remained low until 38.2 ka.Aridification and/or cooling were inferred from the diatom and pollen data at ca. 29.0 ka (Panagiotopoulos et al., 2014;Cvetkoska et al., 2015).The start of the interstadial warming in the DEEP site diatom record is documented between ca.57.5 and 52.0 ka by a successive increase in the relative proportions of C. ocellata (25-35 %, Fig. 2) and BSi concentrations (Fig. 4), as well as decreasing K intensities.However, low BSi, TOC, DAR and DC in assemblages dominated by C. fottii and C. minuscula between 57.5 and 35.0 ka suggest low productivity and/or enhanced decomposition, and similarly, as in Prespa indicate that the interstadial warming was insufficient to enhance the dominance of the epilimnetic, thermophilous species.The decline of C. ocellata to < 5 % and increased abundances of C. fottii (Fig. 2) indicate the return to colder conditions after ca.33.0 ka, culminating in the return to glacial climatic conditions in Lake Ohrid at ca. 28.3 ka.
During MIS 3, a change in the structure and productivity of the diatom assemblages in Lake Prespa core Co1215 was related to the joined effect of the CI/Y5 tephra and the inferred climate aridity/cooling at ca. 39.3 ka (Wagner et al., 2012;Leng et al., 2013;Panagiotopoulos et al., 2014;Cvetkoska et al., 2015).
Similarly, the most notable features of MIS 3 in Lake Ohrid are two maxima in DAR and DC between 39.2 and 38.8 ka and at ca. 29.0 ka, prompted by ca.85 and ca.55 % relative abundance of C. minuscula, respectively (Fig. 2).Both events coincide with high K intensities (Fig. 4), and correspond to the Y5 and Y3 tephra layers in the DEEP site sequence (Leicher et al., 2016).The diatoms indicate that the Y5 and Y3 tephra deposition disturbed the diatom communities by creating conditions of low transparency and high silica availability that enhanced the dominance of small celled Cyclotella species, due to the strong competition for nutrients and light resources.In a case study with leaching experiments on fresh ash samples from Etna eruptions, D'Addabbo et al. (2015) showed that there is only a small impact of tephra on Lake Ohrid water, except for Si and F saturation.Here, for example, the over-dominance of C. minuscula after the Y5 tephra impact in both lakes (this study, Cvetkoska et al., 2015;Jovanovska et al., 2016), is compensated by its small valve size (3-7 µm) as superior adaptation enhancing its competitive strength during such conditions.

5.2.3
The last glacial maximum, MIS 2 (29.0-14.0ka; DZ 2b, a, 1f) During the LGM, Lake Prespa experienced low productivity and low lake levels until 15.7 ka, and the diatom communities were dominated by facultative planktonic and benthic species (Cvetkoska et al., 2014a(Cvetkoska et al., , 2015)).In Lake Ohrid, minima in DAR, DC, and low BSi content between ca.29.0 and 13.6 ka (Fig. 2) indicate low productivity and low temperatures.The low TOC and TIC (Fig. 4) support this interpretation and suggest enhanced decomposition of the organic matter due to cold/dry glacial regime (Wagner et al., 2009;Vogel et al., 2010b;Francke et al., 2016).Similar conditions with winter-ice cover at least in the littoral zone, frequent complete overturn of the lake, substantial opening of the landscape and reduction of vegetation belts were interpreted from multi-proxy analyses of other Lake Ohrid cores Lz1120 and Co1202 (Wagner et al., 2009;Vogel et al., 2010b), and from the DEEP site pollen record as well (Sadori et al., 2016).
www.biogeosciences.net/13/3147/2016/Biogeosciences, 13, 3147-3162, 2016 In terms of diatom response, the similar pattern of glacial assemblages dominated by large, robust valves of C. fottii, was also observed in the previous diatom records from Lake Ohrid (Wagner et al., 2009;Reed et al., 2010;Cvetkoska et al., 2012).Such similarity indicates that this is a "glacialclimate" driven regime of low productivity and diversity, and not just an artifact of taphonomic bias due to weak preservation.In general, both lakes seem to have buffered climate change during the LGM, since no major changes of the diatom communities, and/or species extinctions have been observed in the diatom records.
5.3 Termination I and the Holocene, MIS 1 (14.0 ka to present; DZ 1e-a) Unlike the strong response of Lake Prespa to the wetdry-wet phases between 15.7 and 12.3 ka due to the Bølling/Allerød (B/A) and Younger Dryas (YD) (Panagiotopoulos et al., 2013;Cvetkoska et al., 2014a), a first signal of climate warming subsequent to the LGM in Lake Ohrid derives from an unusual 15 % abundance of "Janus", e.g., heterovalvate cells of C. ocellata at ca. 14.0 ka (Fig. 2).The formation of heterovalvate cells is a phenomenon observed in many diatoms and often related to variation in nutrients, pH or temperatures (McBride and Edgar, 1998;Stoermer, 1967;Teubner, 1995).Their increased abundances at Lake Ohrid DEEP site core can be tentatively related to improved temperature and light conditions at the end of the LGM, as supported by somewhat higher DAR, DC, TOC and BSi between 13.6 and 11.2 ka (Figs. 2, 4).At the same time, low TIC implies frequent mixing of the water column and dissolution of calcite during cold/dry winters (Vogel et al., 2010b;Francke et al., 2016).Indeed, Lacey et al. (2014) and Zhang et al. (2016) interpreted the changes between 12.3 and 11.8 ka in the bio(geochemical) and diatom data set of the "Lini" core Co1262 as reflection of the YD climate reversal.However, this record dates back to 12.4 ka only and does not provide complete insights into the MIS 2/1 transition.Due to the sampling resolution and presence of hiatuses, cores Lz1120 and Co1202 from Lake Ohrid (Fig. 1) only provided incomplete or inconclusive results for this period (Wagner et al., 2009(Wagner et al., , 2012;;Reed et al., 2010;Vogel et al., 2010a;Cvetkoska et al., 2012).We thus tentatively relate the weakness of response observed at the DEEP site to a combination of the depth at the coring location, low temperatures, and low nutrient availability prior to 11.2 ka.
While a rapid water-level increase in Lake Prespa was inferred between 10.0 and 8.0 ka (Panagiotopoulos et al., 2013;Cvetkoska et al., 2014a), in Lake Ohrid DEEP site record, the presence of Stephanodiscus species, indicative of increased nutrient availability, and slightly increased DC and DAR during low SAR imply enhanced productivity between 11.2 and 9.0 ka.A comparable trend of increased productivity was observed at the Lini site (Lacey et al., 2014), while Zhang et al. (2016) found similar trends of water temperatureinduced lake productivity with an interruption between 10.2 and 9.5 ka.
Based on core Co1215 diatom data, sustained moisture availability, high lake levels, and oligo-mesotrophic regime prevailed in Lake Prespa between 7.9 and 1.9 ka (Cvetkoska et al., 2014a(Cvetkoska et al., , 2015)).At Lake Ohrid DEEP site, maximum BSi content, high TIC and TOC concentrations, and diatom assemblages co-dominated by C. ocellata and S. transylvanicus between ca.8.9 and 4.8 ka (Figs. 2, 4) indicate high lake productivity.These conditions are rather indicative of warm mid-Holocene climate than high nutrient supply from Lake Prespa.Rising water-levels at Prespa could have led to higher water supply to Lake Ohrid, but rising lake level would have also led to a decrease in the nutrient concentrations in the water (Matzinger et al., 2006).
The mid-Holocene at Lake Ohrid is punctuated by two synchronous peaks of C. ocellata (ca.5-7 µm) and C. minuscula, with sum abundances of ca.50 and 75 % at ca. 8.5 and 7.9 ka, respectively, co-occurring with peaks in K (Figs. 2, 4).These pronounced signals are likely related to the Mercato tephra layer in the DEEP site sequence (Leicher et al., 2016;Francke et al., 2016) and the "8.2 ka" cold event.But, delineating the exact reasons for the delay of the peaks in diatoms by ca.300 or 400 years requires a high-resolution study.
Around 1.0 ka, the combination of human impact, climate aridification and decreased summer precipitation during the "Medieval Climate Anomaly", MCA was considered as the tipping point that led to significant lake-level decline at Lake Prespa (Cvetkoska et al., 2014a and references therein).In fact, Prespa's decline of ∼ 20 m at the time of the MCA, as also evident in the position of the historic settlements around the lake (Sibinoviç, 1987), is a close analogue to the modelled lake-level decline of 20 m that can potentially increase the P load in Lake Ohrid by almost 30 % (see Fig. 11 in Matzinger et al., 2006).However, the lack of meso-eutrophic species and low TOC, in combination with the high water depth at the DEEP site coring location and the sample resolution do not support enhanced nutrient levels and, by extension, significant influence from Lake Prespa.Moreover, the co-dominance of C. ocellata and C. minuscula between 1.5 ka and present can be related to a combination of warm climate and/or anthropogenic influence at the DEEP site, as also shown from the diatom study of core Co1262 from the Lini site (Zhang et al., 2016).
6 Ecosystems internal dynamics and interactions 6.1 Adaptive cycles: thresholds, regimes, resistance and resilience From the above interpretation and comparison between the lakes, shifts in the diatom communities, and by inference the lakes' response, are triggered by critical thresholds at differ-ent scales of intensity and by different external factors.As pointed out by Cvetkoska et al. (2015), while the diatom response in Lake Prespa is primarily driven by moisture availability, the multi-proxy approach in this study shows that regime shifts in Lake Ohrid are tipped by temperature thresholds, inducing changes in wind activity, light and nutrient availability.
Apart from the encountered differences, two adaptive cycles can be identified at orbital scale: "interglacial and interstadial" and "glacial" cycle.The glacial cycle of Lake Prespa is characterized by much lower productivity than the interglacial/interstadial cycle and dominated by facultative planktonic and benthic species (Cvetkoska et al., 2015).The lake experienced multiple regime shifts on sub-orbital timescale, resulting in low or high lake levels, and changes between (oligo-) meso-and eutrophic regimes.Lake Ohrid's glacial cycle is mainly dominated by the hypolimnetic planktonic C. fottii and is characterized by a lower productivity than the interglacial/interstadial cycles (Fig. 2).Within these cycles, short-term regime shifts in the diatom communities occur at sub-orbital timescale, but the lake remained ultra-to oligotrophic.Based on the presented proxy data, there is no indication for lake-level changes throughout the entire studied period.
More importantly, all diatom data show that during the last 92.0ka both lakes did not undergo catastrophic ecosystem collapses imprinted by loss of the diatom communities, and/or ecosystems functionality, such as happened at other ancient lakes during the LGM, like Baikal and Hövsgöl (Karabanov et al., 2004;Khursevich et al., 2006).
Assessing the capacity of both ecosystems to absorb disturbances and reorganize in order to retain their functional and structural characteristics implies that both lakes have different response mechanisms.The concept of the "ecosystems resilience" (Holling, 1973) is widely accepted, but also controversial because of its definition of mechanisms and components (Carpenter et al., 2001;Scheffer et al., 2001;Walker et al., 2004;Lake, 2013).Here, we consider the "resistance" and "resilience" as different components of the ecosystems response to disturbance; by interpreting the resistance in this case as the capacity of the ecosystem to absorb and withstand the disturbance and the resilience as the capacity to reorganize and/or restore or "recover" from the disturbance (Westman, 1978;Webster et al., 1983).As evident from the diatom and bio(geochemical) data, and the comparison of the lakes (Fig. 4), the amplitude of Lake Ohrid's response, appears smaller and in some cases shifting more gradually into the new regime, like for example at the MIS 5/4 transition (Figs. 2, 4).
This implies a high capacity of Lake Ohrid to absorb the disturbances and retain, and/or reorganize its structure, in order to maintain its functionality.Similarly, once pushed into a new regime, a tipping point of much higher intensity is required so the system can reorganize and/or restore a similar dynamic regime.In the DEEP site record, this pattern is ev-ident from the gradual, prolonged "recovery" periods after climate disturbances, as the gradual transition and muted response to the two-step deglaciation pattern of Termination I (cf.Lowe et al., 2008).
Relatively low intensity disturbances in Lake Prespa, however, are occasionally already sufficient to trigger highamplitude changes in lake levels and productivity (DCA Axes 1 and 2 in Fig. 4), emphasizing a lower lake capacity to absorb disturbance.An example is the abrupt shift from oligo-mesotrophic to eutrophic conditions during MIS 5b (Fig. 4).Interestingly, the relatively short periods of "recovery" after these climate-induced disturbances, like for example the rapid lake-level increase after the LGM, demonstrates that Lake Prespa is a resilient ecosystem.These lakes' specific differences can be attributed to their different physical properties and hydrological balances.

Lake ecosystem interactions
Understanding the interactions between lakes Prespa and Ohrid requires assessment of several parameters which influence the dynamics of their hydrological connection through the karst system of Mt.Galičica: (i) Lake Prespa water balance, (ii) the P load and TP concentration in the lake, (iii) the underground outflow from Lake Prespa, and (iv) the P load from Lake Prespa to Lake Ohrid.Matzinger et al. (2006) assessed these parameters, demonstrating that currently about 20 % of Lake Ohrid water inflow and 7 % of the P load originate from Lake Prespa.However, a large part, ∼ 65 %, of the P load from Lake Prespa is retained in the karst aquifers.Lake Prespa is highly sensitive to water loss and a 5 m lakelevel decrease corresponds to a ∼ 25 % loss of the total lake volume, with significantly increasing concentration of dissolved nutrients (Matzinger et al., 2006).
In order to understand the role of Lake Prespa as a possible driver of lake-level and/or nutrient shifts in Lake Ohrid over the past 92.0ka, we explored the panarchy within the "sister lake system" by identifying possible cross-scale linkages between the adaptive cycles and regime changes of both lakes as lower-scale entities of the landscape dynamics.
From the comparison of the multi-proxy inferred lake levels and productivity of both lakes (see Ecosystem dynamics, Fig. 4), it is evident that the low lake-level and eutrophic conditions at Prespa during the interglacials, such as during MIS 5a and the MCA, do not counterpart with similar regimes at Lake Ohrid DEEP site.On the contrary, Lake Ohrid recharge was probably reduced during the glacial, when lake level at Prespa was low and discontinuous permafrost prevailed in the catchment (Belmecheri et al., 2009).Indeed, Lacey et al. (2014) already discussed the potential scenario for reduced water input from Prespa during the glacials resulting in lower lake water δ 18 O.The palynological analyses of the last ca.0.5 Ma from the DEEP site sequence (Sadori et al., 2016) and pollen based reconstructed LGM temperatures (Peyron et al., 1998) support the existence of discontinuous permafrost and suggest that the karst aquifer system was not completely shut off during the glacials.Distinct water-level reductions are not evident in the DEEP site diatom record from the last ca.92.0 ka, nor have been observed in the geochemistry and diatom data from the previous studies covering the last interglacial-glacial cycle (Wagner et al., 2009;Reed et al., 2010;Cvetkoska et al., 2015;Zhang et al., 2016).
As pointed by Matzinger et al. (2006) the P load to Lake Ohrid changes with decreasing Lake Prespa water level, but it highly depends on the following factors: (i) a smaller subaquatic outflow is reducing the P input, (ii) declining volumes of Lake Prespa lead to increase in P concentration in the lake and (iii) shorter hydraulic residence time of Lake Prespa decreases its P concentration.On glacial-interglacial scale, this would imply that potentially higher P supply to Ohrid due to lower lake level in Prespa is probably compensated by smaller subaquatic outflow.
On a Quaternary scale, it can be argued that climate disturbances leading to changes in lake level and trophic regimes in Lake Prespa were probably insufficient to affect Lake Ohrid's productivity.This decoupling results from the high resistance of Lake Ohrid, or at least the DEEP site, against the disturbances which could have been introduced from Lake Prespa.Moreover, the data show that from MIS 5b to present, Lake Ohrid did not undergo regime changes introduced by Lake Prespa and did not experience a phase of collapse, or release ( ) as defined by Holling (1986).

Conclusions
From the presented data, it is evident that apart from the climatically triggered regime shifts, the lack of similar abrupt peaks in the lake levels and/or productivity in the sediment records of lakes Ohrid and Prespa, there is no clear indication for a direct link between the lakes' regime changes.While Lake Prespa shifted from shallow eutrophic to deep (oligo-)mesotrophic regimes and vice versa, Lake Ohrid changed between ultra oligo-and oligotrophic, deep lake system over the last glacial-interglacial period.
In general, both lakes display different response mechanisms to external disturbances.While Lake Ohrid is more resistant to climate change, Lake Prespa is very sensitive, but a resilient ecosystem, and recovers in a relatively short time.Nonetheless, both lakes did not experience loss of the diatom communities and/or catastrophic ecosystem collapses.
The overall comparison provides sufficient evidence to disregard the theory of Prespa-dependent regime shifts in Ohrid, at least in its central and deepest parts.This may give future research on Lake Ohrid a more confident interpretation of the earlier parts of the DEEP site sediment sequence in terms of climate change.
The Supplement related to this article is available online at doi:10.5194/bg-13-3147-2016-supplement.

Figure 1 .
Figure 1.Map of the northern Mediterranean region showing the location of lakes Ohrid and Prespa.Marked in purple are the locations of the cores: JO2004 (Belmecheri et al., 2009), Lz 1120, Co1202, Co1262 and from field campaigns at Lake Ohrid between 2005 and 2012, cores Co1204 and Co1215 recovered from Lake Prespa during field campaigns in 2009 and 2011;(Wagner et al., 2009(Wagner et al., , 2010(Wagner et al., , 2012) ) and the DEEP site core 5045-1, retrieved during the SCOP-SCO deep drilling campaign at Lake Ohrid in 2013(Wagner et  al., 2014a).Modified fromFrancke et al. (2016).
The most abundant FP and benthic genera are Amphora Kützing, Cocconeis Ehrenberg, Gomphonema Ehrenberg, Diploneis Ehrenberg ex Cleve, Navicula Bory, Staurosira Ehrenberg and Staurosirella D. M. Williams et F. E. Round.Following previous studies of Reed et al. (2010), Cvetkoska et al. (2012), and Wagner et al. (2014a), the dominant planktonic taxa were split in morphotypes, useful to test the link between the longterm ecosystem dynamics and species evolution in the DEEP site record.Cyclotella fottii was initially split in four morphotypes, based on valve size, and size and shape of the central area, while the C. ocellata Pantocsek complex was split in eight morphotypes based on valve size and number of ocelli in the central area.

Figure 4 .
Figure4.Comparison diagram between Lake Ohrid core 5045-1 and Lake Prespa core Co1215 showing the relative abundance data of selected diatom species, the PCA 1 and 2 sample scores; DC and DAR, the sediment accumulation rates (SAR) and the biogenic silica content (BSi, % weight) of Lake Ohrid core DEEP.Prespa DCA 1 and 2 sample scores represent lake level and productivity.Displayed are the diatom zones (DZ 1-4) of Lake Ohrid core DEEP; the dotted lines mark the six major diatom zones of Lake Prespa core Co1215(Cvetkoska et al., 2015).Selected geochemical proxies from Lake Ohrid core DEEP: total organic carbon (TOC), total inorganic carbon (TIC), C / N ratio and potassium (K); data fromFrancke et al. (2016); Lake Prespa core Co1215: TOC, TIC, C / N ratio, K; data fromWagner et al. (2012) and S. pinnata (%), data fromCvetkoska et al. (2015).The MIS boundaries are fromLisiecki and Raymo (2005).

Table 1 .
Locations, sediment recovery and summary of methods used for the diatom analyses of cores 5045-1 from Lake Ohrid and Co1215 from Lake Prespa.

Table 3 .
Spearman's rho (q) coefficient values for the correlations between selected diatom data, and geochemical proxies from Lake Ohrid DEEP site core.
All presented values are significant at p < 0.05." / " low values and/or no significant correlation.