Late Holocene Hydro-Climate Variability in the Eastern Mediterranean: A Spatial Multi-Proxy Approach

A total of thirteen (13) paleoclimatic coastal and hinterland archives of the broader eastern Mediterranean region were collected and examined statistically in search of underlying trends for the period 2800 to 200 BP. For each archive, a proxy record representative of hydro-climatic changes was selected, normalized using z-factors to facilitate intercomparison, and analyzed statistically. Multivariate statistical analysis was performed using a clustering analysis (HCA) and dimension reduction (PCA), which led to groupings of similar records temporally, and allowed the identification of spatially underlying modes of variability. Two main modes of variability were identified, further supporting complex trajectories of paleoclimatic evolution in the region. The first mode was identified for sites presenting a trend from a wetter to an overall drier phase, with respective changes at major phase shifts at 1400 BP and 1100 BP. All sites were from the southern and northern Balkan region, as well as southwestern Turkey. A contrasting dry to wet trend was identified for a site in the Peloponnese (Greece) and the Levant, with a major phase shift at around 750 BP. The inclusion of different proxies from very different environmental settings and the 200-year window has complicated the connection of established short-term climatic events to the study’s findings.


Introduction
Global climate change in the last decades has been a hot topic for geoscientific research and public discussion alike, mainly around human-induced warming. With this warming being unambiguously attributed to anthropogenic greenhouse emissions, quantifying it, making projections, as well as preparing for future change have brought up major challenges. The ultimate goal has been the establishment of future scenarios that can successfully assess the natural as well as social, political, and economic impacts around the globe [1].
The dissimilar effect [2] and heterogeneity of climate change in different regions have emphasized the need for regional assessments. In this regional context, the Mediterranean has been characterized as a climate change ''hot-spot" [3,4], with an urgent need for adaption to its effects [5,6]. Specifically, for the eastern Mediterranean region, which is the area of interest in this study, current warming is projected to result in a severe decrease in precipitation [7].
Instrumental observations of the climate system hardly cover the last two centuries; therefore, the understanding of climate variability has to complementarily rely on reconstructions of past climates [8,9]. Past efforts have led to a dense network of such reconstructions in the circum-Mediterranean and the eastern Mediterranean more specifically [10,11]. Synthesizing the existing data has divulged spatial and temporal variability with dominant south-east versus north-west dipoles throughout the basin's Holocene paleoclimatic evolution [12][13][14]. This evolution has been roughly divided into three phases: a humid phase from the early Holocene until 7000 BP, an intermediate transitional period from 7700 to 5500 BP, and aridification from 5500 BP onwards [15,16]; however, a concise and holistic understanding for the eastern Mediterranean and the whole basin has been very difficult, even more so in the late Holocene. With an increasing amount of information, it is even more imperative to identify trends and fit variability into a temporal context [10,11,17,18].
The main aim of this study is to investigate and detect spatial and temporal patterns underlying the very complex climatic variability in the eastern Mediterranean for roughly the last 2800 years. Potentially, results from this study can provide invaluable aid in future climatic studies and solid ground for the disentanglement of climate dynamics in this climatically forced region.

Materials and Methods
In this study, the archives used were lagoon/lake sediments and cave speleothems. Basic criteria for accepting or rejecting several candidate archives are in some accordance with those in the review by [11] and are: (a) a well-established chronology with at least one measurement for every 200 years and one accepted dating for every 2000 years; (b) a temporal range covering the period from 200 to 2800 BP. Availability of the data is a major parameter in the number of final sites.
For the intercomparison of archives, a single proxy was used for each study. As hydroclimatic variability was the mean of comparison, only proxies reflecting hydro-climate were chosen. The selection of one representative proxy for each site, indicating wetter/drier conditions, was made mostly based on the interpretations of the original studies.
The comparison of unevenly sampled disparate archives was tackled by resampling the sequences into bins of equal time intervals [11,17,18]. A 200-year interval was chosen as a representative time period covering important hydro-climatic change, which has been used in similar [11] and pollen-based studies in the Mediterranean [19,20]. For each 200-year bin, all proxy measurements corresponding to that interval were averaged, yielding a single representative value. In the period of interest (2800 to 200 years), this led to 14 bins representing each archive. Temporal coverage sites with only one empty bin, containing no measurement in 200 years, were accepted, and with more than one rejected. If the missing bin covered the last 200 years (200 BP to present, with present being 1950), the bin was left empty; otherwise, the value was calculated averaging the previous and next bin.
To compare proxies of different environmental sensitivity, the use of z-scores was implemented. In each sequence, the proxy series mean for 2800 BP to present was subtracted from the binned value and then divided by the standard deviation of the sequence. Through this, all proxies were normalized with 0 as the mean and the standard deviation as the measurement unit of distance, and therefore variability, from the mean. Positive z-scores were associated with a relatively drier signal from the proxy, and negative z-score values with wetter signals, respectively. To ensure consistency among the datasets, in some cases z-scores were multiplied with (−1) if a wetter period was represented by positive values and vice versa. The interpolation at intermediate intervals, in time, ultimately filters the high-frequency variation and therefore blurs any short-term events with a duration smaller than the bin's. Short, extreme events can disproportionally affect a bin, resulting in a false interpretation for the 200-year interval.
To detect similar temporal fluctuations, a hierarchical cluster analysis was performed on the z-score normalized data (Table S1). For the dissimilarity measure, the Spearman distance was chosen. A correlation-based distance was used as the grouping was based on overall fluctuations regardless of magnitude; in this case, positive z-scores (a dry signal) and vice versa. For the intercluster dissimilarity and the dendrogram, the average linkage was used.
Following the methodology in [17], we strived to find a composite record for each cluster in the hope of being representative of a temporal trend and useful as a means of comparison to the subsequent PCA. Within each cluster, an average z-score was calcu-lated for every 200-year interval period, resulting in a single representative record for the whole cluster.
Lastly, principal component analysis was used to infer trends in our data. The principal components representing the majority of variance in the data were used together with the cluster composites in determining trends underlying the proxy variability. A total number of 13 archives were considered in the present study and are presented in Table 1.

Clusters
Using hierarchical clustering analysis, all sequences were attributed to three distinct clusters based on their long-term variation ( Figure 1). Each cluster consists of sites with temporally similar patterns of relatively wetter (−) or drier (+) periods, according to the z-score-transformed data. For each cluster, a composite sequence is presented, representing the mean trend throughout the cluster in time ( Figure 1).
The first cluster includes four (n = 4) sites: (1) Gialova lagoon [21], (8) Lake Cubuk [28], (11) Jeita cave [32], and (12) Lake Kinneret [33]. The averaged time series (composite series) for the first cluster indicates drier periods from 2800 to 2600 BP, a less dry and closer to the mean period from 2600 to 2200 BP, and an increase in relatively drier conditions reaching the cluster's maximum from around 1400 to 1200 BP. From 1200 to 800 BP, a stable, less dry climate is followed by transition to a wetter period from 800 to 200 BP. The composite is characterized by an overall drier to wetter trend ( Figure 2) The second cluster includes two (n = 2) sites: (2) Lake Lerna [22] and (13) Soreq cave [34][35][36]. The averaged time series for this cluster indicates wetter conditions relative to the mean from 2800 to 2600 BP, more dry conditions from 2600 to 2400 BP, and considerably drier conditions from 2400 to 1800 BP. From 1800 to 1000 BP, wetter conditions are apparent, with the wettest period of the composite being from 1600 to 1400 BP. From 1000 to 800 BP, dry conditions are evident and are followed by a considerably wetter period from 800 to 400 BP. An abrupt change to wetter conditions in the bin spanning from 600 to 400 BP is also apparent in the first composite. In the period 400 to 200, BP a shift to a drier climate is recognized. The first cluster includes four (n = 4) sites: (1) Gialova lagoon [21], (8) Lake Cubuk [28], (11) Jeita cave [32], and (12) Lake Kinneret [33]. The averaged time series (composite series) for the first cluster indicates drier periods from 2800 to 2600 BP, a less dry and closer to the mean period from 2600 to 2200 BP, and an increase in relatively drier conditions reaching the cluster's maximum from around 1400 to 1200 BP. From 1200 to 800 BP, a stable, less dry climate is followed by transition to a wetter period from 800 to 200 BP. The composite is characterized by an overall drier to wetter trend ( Figure 2) The second cluster includes two (n = 2) sites: (2) Lake Lerna [22] and (13) Soreq cave [34][35][36]. The averaged time series for this cluster indicates wetter conditions relative to the mean from 2800 to 2600 BP, more dry conditions from 2600 to 2400 BP, and considerably drier conditions from 2400 to 1800 BP. From 1800 to 1000 BP, wetter conditions are apparent, with the wettest period of the composite being from 1600 to 1400 BP. From 1000 to The last cluster consists of seven sites (n = 7): (3) Aliki lagoon [23], (6) Lake Ohrid [26], (4) Lake Trichonida [24], (5) Lake Butrint [25], (7) Lake Dojran [27], (9) Sofular cave [29], and (10) Lake Golhisar [30,31]. From 2800 to 1200 BP, generally wetter composite mean conditions are inferred, with a considerably wetter bin from 1600 to 1400 BP. This is followed by a drier phase until 200 BP. An overall trend from relatively wetter to drier conditions is apparent for this cluster (Figure 2). climate is recognized.

Principal Components
The first principal component explains 44.2% of the variation in our data. It mainly emphasizes centennial-scale trends from progressive drying followed by progressive wetting ( Figure 3A). Temporally, it shows drier conditions from 2800 to 1800 BP ( Figure 3A). From 1800 to 1400 BP, even drier conditions are evident, followed by a change to overall wetter conditions from 1000 BP onward. PC1 resembles composite 1 by contrasting an overall drier to wetter phase. This is further indicated by the coherence of sites in the first cluster presenting positive loadings with the first component ( Figure 3B). A Spearman correlation test between composite 1 and the scores of the first component (Spearman ρ = 0.67, p-value = 0.016) further indicates substantial correlation.

Principal Components
The first principal component explains 44.2% of the variation in our data. It mainly emphasizes centennial-scale trends from progressive drying followed by progressive wetting ( Figure 3A). Temporally, it shows drier conditions from 2800 to 1800 BP ( Figure 3A). From 1800 to 1400 BP, even drier conditions are evident, followed by a change to overall wetter conditions from 1000 BP onward. PC1 resembles composite 1 by contrasting an overall drier to wetter phase. This is further indicated by the coherence of sites in the first cluster presenting positive loadings with the first component ( Figure 3B). A Spearman correlation test between composite 1 and the scores of the first component (Spearman ρ = 0.67, p-value = 0.016) further indicates substantial correlation.
The second principal component explains 20.1% of the variance in our dataset. PC2 emphasizes a progressive change to wetter conditions followed by a progressive drying ( Figure 3A). More specifically, wetter conditions cover the period from 2600 to 2200 BP, with a drier period following for the rest of the studied record until 400 BP ( Figure 3A). Some similarity occurs with the general trend of cluster 1, from wetter to drier. However, the loadings associated with PC2, and the discrepancy in the time of change from overall wet to overall dry at around 2200 to 2000 BP, point to a subcluster of cluster 3 (Aliki lagoon, Lake Ohrid, and Lake Trichonida), together with Lake Kinneret, Lake Cubuk, and Gialova (cluster 1), as well as Lake Golhisar in cluster 1 to a lesser extent ( Figure 3B). Intercluster dissimilarity is therefore evident. It should also be noted that Lake Lerna, in proximity to Aliki, Trichonida, and Gialova, shows negative loadings, with opposite signals for the period 2600 to 2200 BP ( Figure 3B). wet to overall dry at around 2200 to 2000 BP, point to a subcluster of cluster 3 (Aliki la-goon, Lake Ohrid, and Lake Trichonida), together with Lake Kinneret, Lake Cubuk, and Gialova (cluster 1), as well as Lake Golhisar in cluster 1 to a lesser extent ( Figure 3B). Intercluster dissimilarity is therefore evident. It should also be noted that Lake Lerna, in proximity to Aliki, Trichonida, and Gialova, shows negative loadings, with opposite signals for the period 2600 to 2200 BP ( Figure 3B).
The third principal component explains 13.6% of our data. It emphasizes a progressively drier climate from 2000 to 1400 BP and a wetter phase from 1400 to 600 BP, followed by an increase in arid conditions thereafter ( Figure 3A). This component was not further linked to our data.

Discussion
The clustering process resulted in groupings partly based on similar temporal fluctuation to the 2800-year mean in each study (z-scores with a 200-year step are presented in Figure S1). The unifying trend in the first cluster seems to be a dry to wet millennial trend with a major shift to wetter conditions at around 800 to 600 BP (Figure 4). The studies in this cluster are spatially very distant, grouping parts of SW Greece, northern Turkey, and the northern Levant. At Gialova lagoon, overall stable warmer (2600 to 2300 BP) to increasingly drier conditions from 2000 to around 1300 BP are reported [21]. Similar aridity trends are generally inferred for the Levant until 1700 BP and for the eastern Mediterranean in general from 2400 to 1600 BP [10]. More specifically, at Jeita cave, the drier conditions last until 1200 BP while the 2800-1400 BP dry period is described as ''a severe millennial event'' [32] of the Holocene. At Lake Kinneret, the drier period extends from 2300 to 750 BP [33]. At Lake Cubuk, drier conditions are, however, indicated only until 2300 BP (Near East aridification phase) and 1600 to 1100 BP [28], and this location therefore stands out in the cluster. Proximate Lake Iznik has extended aridification but only The third principal component explains 13.6% of our data. It emphasizes a progressively drier climate from 2000 to 1400 BP and a wetter phase from 1400 to 600 BP, followed by an increase in arid conditions thereafter ( Figure 3A). This component was not further linked to our data.

Discussion
The clustering process resulted in groupings partly based on similar temporal fluctuation to the 2800-year mean in each study (z-scores with a 200-year step are presented in Figure S1). The unifying trend in the first cluster seems to be a dry to wet millennial trend with a major shift to wetter conditions at around 800 to 600 BP ( Figure 4). The studies in this cluster are spatially very distant, grouping parts of SW Greece, northern Turkey, and the northern Levant. At Gialova lagoon, overall stable warmer (2600 to 2300 BP) to increasingly drier conditions from 2000 to around 1300 BP are reported [21]. Similar aridity trends are generally inferred for the Levant until 1700 BP and for the eastern Mediterranean in general from 2400 to 1600 BP [10]. More specifically, at Jeita cave, the drier conditions last until 1200 BP while the 2800-1400 BP dry period is described as ''a severe millennial event" [32] of the Holocene. At Lake Kinneret, the drier period extends from 2300 to 750 BP [33]. At Lake Cubuk, drier conditions are, however, indicated only until 2300 BP (Near East aridification phase) and 1600 to 1100 BP [28], and this location therefore stands out in the cluster. Proximate Lake Iznik has extended aridification but only until 2000 BP [37]. The period following 1300 to 700 BP is often termed as the medieval climate anomaly (MCA). This period in the Gialova record is divided into a wetter phase (1300 to 1000 BP) and a drier phase (until 700 BP). For the Levant, the MCA is characterized by wetter conditions. The period of extended drought in the cluster comes to an end with an overall shift to wetter conditions after 750 BP, coinciding with the Little Ice Age (LIA) period. Although, this is not the case for Lake Cubuk, which seems to be a misfit in this cluster. until 2000 BP [37]. The period following 1300 to 700 BP is often termed as the medieval climate anomaly (MCA). This period in the Gialova record is divided into a wetter phase (1300 to 1000 BP) and a drier phase (until 700 BP). For the Levant, the MCA is characterized by wetter conditions. The period of extended drought in the cluster comes to an end with an overall shift to wetter conditions after 750 BP, coinciding with the Little Ice Age (LIA) period. Although, this is not the case for Lake Cubuk, which seems to be a misfit in this cluster. The second cluster consists of Lake Lerna [22] and Soreq cave [34][35][36] (Figure 5). A similarity in hydro-climatic variability among these sites has been outlined by the authors in [22] and has been successfully pointed out by the clustering process. Rather more interesting is the division of the Peloponnesian peninsula, between the first and the second cluster, and likewise of the northern and southern Levant [32]. The division of the peninsula has been attributed to changes in the North Sea-Caspian atmospheric circulation pattern (NCP) [38,39], with (+) phases leading to more precipitation in the north-eastern Peloponnese and drier conditions in the south-west and vice versa. The second cluster consists of Lake Lerna [22] and Soreq cave [34][35][36] (Figure 5). A similarity in hydro-climatic variability among these sites has been outlined by the authors in [22] and has been successfully pointed out by the clustering process. Rather more interesting is the division of the Peloponnesian peninsula, between the first and the second cluster, and likewise of the northern and southern Levant [32]. The division of the peninsula has been attributed to changes in the North Sea-Caspian atmospheric circulation pattern (NCP) [38,39], with (+) phases leading to more precipitation in the north-eastern Peloponnese and drier conditions in the south-west and vice versa.
In contrasting states to the first cluster and an intermediate second, the sites in the third cluster present a trend from wetter to overall drier conditions ( Figure 5). The sites in this cluster geographically cover the southern Balkans (except the Peloponnese) and northern and south-western Turkey, contrasting the first cluster and resembling a northsouth/east-west seesaw in the region. PC2, Lake Trichonida [24], Lake Ohrid [26], and Aliki Lagoon [23] are similar and show intercluster variability. In the same manner, similarities are shown between Lake Butrint [25] and Lake Dojran [27] and to a lesser extent Sofular cave [29] and Lake Golhisar [30]. For the first three sites, an overall wetter phase shifts to broader drier conditions around 1400 BP, whereas for the other sites in the cluster, this transition happens at around 1100 BP. In contrasting states to the first cluster and an intermediate second, the sites in the third cluster present a trend from wetter to overall drier conditions ( Figure 5). The sites in this cluster geographically cover the southern Balkans (except the Peloponnese) and northern and south-western Turkey, contrasting the first cluster and resembling a northsouth/east-west seesaw in the region. PC2, Lake Trichonida [24], Lake Ohrid [26], and Aliki Lagoon [23] are similar and show intercluster variability. In the same manner, similarities are shown between Lake Butrint [25] and Lake Dojran [27] and to a lesser extent Sofular cave [29] and Lake Golhisar [30]. For the first three sites, an overall wetter phase

Conclusions
In this study, a total of 13 paleoclimatic reconstruction sites were used from the southern Balkans, Turkey, and the Levant to disentangle the paleoclimatic variability of the eastern Mediterranean for the last 2800 years. According to the authors, the proxy best representing palaeohydrological change (wetter/drier conditions) was used from each study. All proxy sequences were filtered to 200-year bins and normalized. Clustering analysis and dimension reduction were implemented in order to identify long-term unifying trends. The study's findings can be summarized as:

•
Two main long-term modes of variability were identified in the studied sites, indicating that no single trajectory of climate variability is evident for the whole eastern Mediterranean. • A wet to dry millennial trend with a major phase shift at 1400 BP and 1100 BP was identified mainly for the northern and northwestern part of the study area, including the Balkans (except the Peloponnese) and Turkey. A contrasting dry to wet millennialscale with a major phase shift coinciding with the LIA at 750 BP can be extracted for southern Greece and the Levant. This might fit in the overall north-south/east-west seesaw representation of Mediterranean spatial variability. • Complex interregional variability is evident and widespread at the centennial level throughout the Mediterranean area. Intercluster variability is apparent in all clusters, indicating microclimates or site-specific responses.

•
The use of hierarchical clustering analysis and principal component analysis has helped identify long-term unifying trends in dissipating climatic archives. • The number of sites, incorporating different proxies and different environments, has led to difficulties in short-term correlations, even at a 200-year interval level.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.