Key Environmental Factors Controlling Planktonic Foraminiferal and Pteropod Community’s Response to Late Quaternary Hydroclimate Changes in the South Aegean Sea (Eastern Mediterranean)

Abstract

A multidisciplinary study was conducted in order to investigate the environmental factors affecting the planktonic foraminiferal and pteropod communities of the south Aegean Sea. Aspects of the Late Quaternary paleoceanographic evolution were revealed by means of quantitative analyses of planktonic foraminiferal and pteropod assemblages (including multivariate statistical approach; principal component analysis (PCA)), the oxygen (δ¹⁸O) and carbon (δ¹³C) isotopic composition of planktonic foraminifera and related paleoceanographic (planktonic paleoclimatic curve (PPC), productivity (E-index), stratification (S-index), seasonality) indices, extracted by the gravity core KIM-2A derived from the submarine area between Kimolos and Sifnos islands. Focusing on the last ~21 calibrated thousands of years before present (ka BP), cold and eutrophicated conditions were identified during the Late Glacial period (21.1–15.7 ka BP) and were followed by warmer and wetter conditions during the deglaciation phase. The beginning of the Holocene was marked by a climatic amelioration and increased seasonality. The more pronounced environmental changes were identified during the deposition of the sapropel sublayers S1a (9.4–7.7 ka BP) and S1b (6.9–6.4 ka BP), with extremely warm and stratified conditions. Pteropod fauna during the sapropel deposition were recorded for the first time in the south Aegean Sea, suggesting arid conditions towards the end of S1a. Besides sea surface temperature (SST), which shows the highest explanatory power for the distribution of the analyzed fauna, water column stratification, primary productivity, and seasonality also control their communities during the Late Quaternary.

1. Introduction

The Aegean Sea is an ideal archive to investigate climatic evolution at both global and local scale because of its intermediate position between the higher- and lower-latitude climate systems [1,2,3], high sedimentation rate marine records compared to the open Mediterranean Sea [4,5,6,7], and its paleo-latitudinal and land-locked configuration [8,9]. Such marginal seas are often more responsive to paleoceanographic and paleoclimatic changes than global oceans, with climatic signals to be recorded in an amplified fashion in Mediterranean properties such as temperature, salinity and specific elemental concentrations, because of their smaller size and partial isolation [10], and therefore can be considered miniature oceans. In addition to interactions with the Black Sea, northern Aegean, and Levantine basins with remote and local atmospheric forcing [1,11], the south Aegean Sea is characterized by intense biogeochemical contrasts in its hydrology in response to a climatic gradient from mid-latitude to subtropical regimes that appears to be very sensitive to climate changes. However, most of the current paleoclimatic and paleoceanographic studies are still limited to deep marine records [6,12,13,14,15], and consequently little is known about the continental shelf and/or coastal areas within this marginal Sea [16,17,18,19,20,21].

Environmental changes related to the different water column and/or sediment characteristics can be recorded virtually instantaneously in paleoceanographic proxy data, such as stable isotope and other geochemical ratios [22,23,24,25,26,27,28,29,30,31], and micro-fossil abundances, such as planktonic foraminifera and pteropods [15,17,28,32,33,34,35,36]. This makes them extremely valuable for both stratigraphic correlations and paleoenvironmental/paleoclimate reconstructions [15,23,28,29,32,35,37,38,39,40,41,42,43,44,45]. Their significance in the study of modern and past marine ecosystems in the eastern Mediterranean Sea is well underlined [20,31,46,47,48,49,50,51]. Particularly, they are used as indicators of temperature, salinity, density, and nutrient content of the water column, making it possible to identify past circulation through the sedimentary record [7,15,29,47,52,53] and detect long- and short-term paleoclimatic and paleoceanographic changes in the study area [6,10,12,15,16,53,54,55] during the last glacial cycle.

Pteropods are widespread and abundant in the global ocean and entirely adapted to a pelagic life cycle [56,57]. Owing to the aragonite nature of their shells which increases their weight as settling particles and hence their sinking speed [58], their deposition is expected to be close to their habitat [59]. Particularly in the Mediterranean Sea, preservation of pteropods shells is excellent as a result of the relatively shallow water, high bottom water temperatures, and probably the limited number of mud feeders [60]. A considerable number of studies (e.g., [33,34,61,62,63,64,65,66]) have shown that Late Quaternary pteropod assemblages and their distribution pattern in the world oceans have changed with temperature and the overall climatic conditions that also affect the aragonite compensation depth (ACD). Recent studies [67] have shown that modern eastern Mediterranean pteropod communities are found to be more abundant than in those at western Mediterranean Sea. Their abundances are positively correlated with the aragonite saturation state (Ωar), O₂ concentration, pH, salinity and temperature, and negatively correlated with nutrient concentrations [67,68]. However, pteropod assemblages and their distribution in the Aegean Sea during the Late Quaternary are poorly documented.

The present study focuses on identifying and describing key environmental factors that control Late Quaternary planktonic foraminiferal and pteropod distribution in the south Aegean Sea, based on marine sediments retrieved by a 2-m long gravity (KIM-2A) core. In addition, paleoclimatic data were revealed from their distribution patterns coupled with variations of oxygen and carbon isotopic signals. The combination of the above data enables speculation on the factors’ response to the climatic changes.

2. Regional and Climatic Setting

The Aegean Sea is in the northern sector of the eastern Mediterranean (Figure 1a), between the Turkish coastline to the east, the Greek mainland to the north and west, and bounded on the south by the island of Crete and Cretan Arc. It is connected to the Black Sea through the Straits of Bosporus and Dardanelles, and to the Levantine Sea through several larger and deeper straits between Peloponnesus, the islands of Crete and Rhodes, and south-western Turkey (Figure 1b). It is characterized, in general, by a cyclonic water circulation, although the most active dynamic features are the mesoscale cyclonic and anticyclonic eddies, either permanent and/or recurrent [8]. It is separated into two major sub-basins with different climatic conditions: the “north” and the “south” Aegean Sea. The north is more humid than the semiarid south [8].

The south Aegean (extend in between 35°N and 37°N) is one of the most oligotrophic areas in the Mediterranean Sea [8], with its surface water circulation mostly affected by arid climatic conditions, while it is also modulated by the effect of the Cretan gyre [69]. It mainly consists of the Cretan basin and the shallow shelf of the Cyclades Plateau, along with the Myrtoan Sea at the NW part of the region (Figure 1b). Milos and Kimolos islands lie in the westernmost sector of the Cyclades plateau. Both islands are part of the south Aegean Volcanic Arc; the most important geological structure of the Aegean Sea. The submarine area between Kimolos and Sifnos Islands is characterized by a rather complex relief [70] as it is related to the volcanic arc. The sedimentary and Quaternary tectonic evolution of the aforementioned submarine area have been studied previously [71], and it is separated into three parts (the southern, the eastern, and the northern; [70]). The core studied derives from a submarine depression located at the northern part.

3. Materials and Methods

3.1. Location and Sampling Strategy of Core KIM-2A

A gravity core (KIM-2A; 200 cm length) was recovered with R/V Aegaeo in May 2015 from the south Aegean Sea and covers the Late Quaternary (the last ∼21 ka BP). Core KIM-2A was collected from the Cyclades plateau, from a submarine depression (640 m water depth) located north of Kimolos Island (36°95.464′N, 24°06.354′E), as shown in Figure 1b. After splitting, the half core was stored as archival material, whereas the working half has been sampled for multiple analyses (micropaleontology, sedimentology, and geochemistry).

Fifty-eight samples were taken and used for paleontological (planktonic foraminifera and pteropods) and geochemical (Total Organic Carbon; TOC, stable oxygen and carbon isotopes; δ¹⁸O, δ¹³C) analyses between 196 cm and 4 cm. Samples were taken throughout the core but with different sampling resolutions. For the first gray interval corresponding to S1a, the sampling resolution was 1 cm; beyond these intervals the sampling was 2 cm, whereas the top 40 cm were sampled at a mean 8 cm resolution.

3.2. Micropaleontological Analyses

All samples were prepared following standard micropaleontological procedures. For the faunal analyses, approximately 3 cm³ of dried sediment was washed and sieved through a 63 μm screen, and residues were dried in an oven at 50 °C. Qualitative and quantitative analyses have been performed on both planktonic foraminiferal and pteropod assemblages for the >125 μm size fraction, split into aliquots, each one containing at least 300 specimens. It should be noted that the 125 μm fraction was selected since is the most common studied fraction in relevant investigations within the Aegean Sea, which analyze the Late Quaternary foraminiferal record [10,18,21,32,36], and implement a paleoclimatic analysis [6,15,20,72]. All shells were handpicked, identified following [73], counted in each sample, and then converted into percentages, based on the extrapolation of a counted split. Planktonic species with phylogenetic affinities and similar ecological characteristics [35,74] were counted together and grouped to better interpret distribution patterns. On this regard, all Globigerinoides ruber morphotypes (“Normal,” “Platys,” “Elongate,” and “Twin” types; [32]) were plotted together, distinguishing only the “alba” and “rosea” varieties due to their different ecological characteristics [75]. Furthermore, Globigerina bulloides group includes the species G. bulloides and Globigerina falconensis, the Globigerinoides sacculifer group includes Globigerinoides trilobus and G. sacculifer, and the Globigerinella siphonifera group includes the species Globigerinella aequilateralis, Globigerina calida, and Globigerina digitata. The species Globigerinita glutinata includes the morphotypes with and without bulla. Within the group of Neogloboquadriniids, two types were discerned: Neogloboquadrina pachyderma and Neogloboquadrina dutertrei.

3.3. Total Organic Carbon and Stable Isotopes

Total organic carbon (TOC) was determined based on the [76] methodology with the [77] adaptation at the Laboratory of the Hellenic Survey of Geology and Mineral Exploration (H.S.G.M.E.). For stable oxygen and carbon (δ¹⁸O, δ¹³C) isotope measurements, 30 specimens of the planktonic species G. ruber f. alba were picked from the 250–300 μm size fraction. In particular, we exclusively used the morphotype “Normal” of [32] (equivalent to G. ruber sensu stricto [75]) in order to minimize potential morphotype-specific differential responses in stable isotope compositions [24,78,79]. This narrow size fraction was used to minimize ontogenetic and growth rate effects on shell geochemistry [80]. The analyses were carried out at the Laboratory of Geology and Geophysics at Edinburgh University. Foraminiferal δ¹⁸O and δ¹³C data were calibrated to National Bureau of Standards 19 (NBS19), and the isotope values are reported in ‰ relative to Vienna Pee Dee belemnite scale. The external standard errors of the stable carbon and oxygen isotope analyses are <0.06‰ and 0.08‰, respectively.

3.4. Chronology

The chronology of the core KIM-2A is based on five accelerator mass spectroscopy radiocarbon (AMS ¹⁴C) dates (on G. ruber tests at Beta Analytics laboratories), supplemented by tie-points which correspond to well-dated bio-litho-stratigraphic horizons (e.g., sapropel and planktonic foraminiferal biozone boundaries) from nearby Aegean cores (

Table 1. Calibrated radiocarbon dates (AMS ¹⁴C) and chronostratigraphic control points’ dates.

AMS & Chronostratigraphic Control Points	Depth (cm)	Conventional Radiocarbon Age (BP)	Two Sigma Calibrated Age Range (BP)	Mean Calibrated Age (ka BP)	References
Beta—425634	14.5	4890+/−30	4845–5325	5.08
Ia/Ib boundary	20			5.2	[84]
Beta—425635	28	5320+/−30	5444–5855	5.65
S1b top	40			6.4	[7,16]
Beta—425636	50.25	6790+/−30	7036–7292	7.16
S1b base	52.5			7.3	[7,16]
S1a top	65.5			7.9	[7,16]
Beta—425637	79.5	8320+/−30	8532–8883	8.71
S1a base	89			10	[7,16]
Ic/II boundary	109			11.3	[84]
II/III boundary	153			15.5	[84]
δ18OG. ruber depletion	159			15.9	[52]
Beta—425638	195	18890+/−70	21962–22508	22.24

). Conventional ¹⁴C ages were calibrated by means of the Calib version 7.0.2 software [81] and the Marine13 data set with a regional reservoir age correction (ΔR) of 139 ± 40 years for the S1 interval [82] and of 58 ± 85 years outside the sapropel S1 interval [83]. The chronology adopted in this study for core KIM-2A was derived from a polynomial fit through the calibrated dates and the chronostratigraphic control points’ dating as shown in Table 1. Hereafter, ages in this study are reported in calibrated thousands of years before present, notated ka BP.

3.5. Multivariate Statistical Analyses

Principal component analyses (PCA) is used to reduce the dimensionality of a multivariate data set to a few principal factors that determine the distributions of species. For this analysis all raw data for the totality of the samples and specimens were used. Raw data were processed using PAST (2.17) multivariate statistical software package of [85]. The resulting factor scores show the contribution of each factor in every sample, and therefore the down-core contribution of each factor. The total number of factors was defined by minimizing the remaining “random” variability, and by the possibility to relate the factors to modern hydrographic conditions and planktonic foraminiferal and pteropod ecology.

3.6. Paleoceanographic Indices

The planktonic foraminiferal relative distributions were used as a first-order estimate of sea-surface temperature (SST) variations. An index of the SST variations was constructed based on the down-core variation of planktonic foraminiferal abundances, referred to as planktonic paleoclimatic curve (PPC). The PPC was obtained by the formula 100 × (w − c)/(w + c), where w represents the warm-water indicators (G. ruber f. alba, G. ruber f. rosea, Orbulina universa, G. sacculifer gr., Globoturborotalita rubescens, and G. siphonifera gr.), and c the cold-water indicators (Globorotalia inflata, Globorotalia truncatulinoides, Turborotalita quinqueloba, G. glutinata, and Globorotalia scitula). The eutrophication index (E-index; [7]) was estimated using the sum of the eutrophic species (N. pachyderma, N. dutertrei, G. bulloides, T. quinqueloba, and G. inflata) versus the sum of the eutrophic plus oligotrophic species (G. ruber alba, G. ruber rosea, G. rubescens, G. sacculifer, O. universa, and G. siphonifera). The down-core ratio between G. bulloides and G. ruber was also estimated showing the degree of the stratification of the upper water column [86]. Following [7], this ratio is referred to as S-index, and its values reflect periods of strong summer stratification of the water column where oligotrophic taxa dominate (low values) and/or periods of strong winter mixing of the water column where eutrophic taxa dominate (high values).

4. Results

4.1. Lithological Description, Time Stratigraphic Framework, and Sedimentation Rates

Lithologically, the study core contains a distinct organic-rich dark interval, divided into two separate sub-units (S1a and S1b respectively), representing the regional expression of the most recent sapropel S1 [87,88] (Figure 2a). From the bottom up to 87.5 cm, light gray clay can be observed (Munsell soil color 5Y 7/1). The following 24 cm are of gray color (5Y 5/1) and correspond to the lower sub-unit (S1a) of the sapropel. Between 63.5 cm and 52 cm light gray clay (5Y 7/1) can be observed, indicative of the sapropel interruption (S1i). From 52 cm to 40 cm mud of gray color (5Y 5/2) characterizes the upper sapropel sub-unit (S1b). The clay continues up to 13.5 cm with a light gray clay (5Y 7/1) color. From this point to the top of the core watery clay of light-yellow color (5Y 7/4) can be observed (Figure 2a).

Chronology adopted in this study for core KIM-2A derived from a polynomial fit through the five AMS ¹⁴C datings mentioned above, and the time markers correlative to the start and end of sapropel deposition, as well as the planktonic foraminiferal biozone Ia/Ib, Ib/Ic, Ic/II and II/III boundaries of [84] in the Mediterranean Sea (Table 1). The resulting age model is illustrated in Figure 2b and is consistent with the multi-proxy chronological framework of [37] for the Aegean sediment cores, thus supporting the robustness of our chronology. For the ages of the onset and termination of sapropel deposition, we have used the relevant ages from the nearby core NS-14 [7,16,17]. Additional control points related to the planktonic foraminiferal biozone boundaries were used for the age model construction, since they are useful chronologic standards for dating late Quaternary sequences in the central and eastern Mediterranean (including the Aegean Sea; [6,7,52,89] due to their Mediterranean-wide applicability and synchroneity [90,91,92].

Based on the relative abundances of the planktonic foraminifera species (Figure 3) the interval between 196 cm and 153 cm is characterized by the dominance of N. pachyderma, T. quinqueloba, G. scitula and G. glutinata with additional components the species G. ruber f. alba and G. bulloides. This glacial fauna corresponds to assemblage III and it has been recognized throughout the Mediterranean Sea (e.g., [6,52]). The interval between 153 cm and 109 cm corresponds to assemblage II and is characterized by high relative abundances of Neogloboquadrinids and G. glutinata and the presence of and G. inflata. In the lower part of this interval G. ruber f. rosea appears for the first time. Subzone Ic (109–40 cm) was identified by the warm subtropical species (G. ruber f. rosea, G. siphonifera gr. and O. universa). In addition it includes abundant G. bulloides and G. rubescens specimens. The sharp increase in the abundance of G. inflata at 41.5 cm is inferred to mark the onset of the Ib subzone. The Ia/Ib boundary (at 20 cm) is marked by the decrease of the latter species along with the decrease in the N. pachyderma abundance. In this core the Bioevent “Start of δ¹⁸O_ruber depletion T1a” of [52], was also detected at 159 cm.

According to our proposed age model (Figure 2b), the sedimentary horizons sampled in this study span the interval from the late glacial period, and the subsequent transition (Termination 1; T1) to the middle Holocene (Northgrippian stage) (i.e., ~5–22 ka BP). The average sedimentation rate is 11.86 cm/ka, and is in good agreement to those reported in the marginal Aegean basin [4,7,13,52]. These sedimentation rates were derived from the age model, assuming that the sediment accumulation has been moderately consistent throughout each interval. In particular, the average sedimentation rates are 8.10 cm/ka for the late glacial, 10.30 cm/ka for the Termination T1, and 16.46 cm/ka for the Holocene.

4.2. Planktonic Foraminifera Distribution Pattern

The qualitative analysis revealed 19 planktonic foraminiferal species lumped into 14 groups: G. ruber f. alba, G. ruber f. rosea, G. sacculifer group, G. bulloides group, G. siphonifera group, G. scitula, G. truncatulinoides, G. inflata, O. universa, T. quinqueloba, G. glutinata, G. rubescens, N. pachyderma, N. dutertrei. The down-core stratigraphic distributions of their relative abundance are shown in Figure 3.

From the bottom of the core up to 153 cm, the fauna is characterized by high relative abundances of the species G. ruber f. alba (22–34%), N. pachyderma (10–40%), T. quinqueloba (3–30%), and G. bulloides gr. (8–33%). Additional components of the fauna are N. dutertrei, G. glutinata, and G. scitula. After the 153 cm, an abrupt decline in Neogloboquadrinids and T. quinqueloba can be observed, while G. bulloides follows an opposite trend with relative abundance of 33%. From that point up to 127 cm N. pachyderma and N. dutertrei reach their maximum abundances (45% and 15% respectively), whereas T. quinqueloba is still present but with low percentages (up to 10%). The relative abundance of G. bulloides does not exceed 16% and G. ruber f. alba is also present but with equally low frequency (lower than 20%). Between 127 cm and 100 cm the fauna is characterized by the dominance of G. bulloides, N. pachyderma, N. dutertrei, G. ruber f. alba, G. glutinata, and G. rubescens. Additional components, with lower percentages, are the species T. quinqueloba and G. inflata. The interval between 100 cm and 41.5 cm is characterized by a shift in fauna. Most of the species dominating the previous interval are decreasing or becoming absent (T. quinqueloba, Neogloboquadrinids, G. glutinata and G. inflata). Prevailing species of this interval are O. universa (44%), G. ruber f. rosea (38%), G. bulloides gr. (27%), G. siphonifera gr. (20%), and G. sacculifer gr. (14%). In the final segment of the core, from 41.5 cm to the top, the sampling presents poor resolution, but certain significant changes can be observed. The G. ruber f. alba, G. inflata and G. truncatulinoides present a peak (57%, 31%, and 7% respectively) and the relative abundances of G. ruber f. rosea, O. universa, and G. siphonifera gr. are decreasing. N. pachyderma, G. inflata, and G. glutinata are re-appearing, but with low percentages.

4.3. Pteropod Distribution Pattern

All samples examined for foraminifera include significant amounts of aragonitic pteropods indicating excellent preservation (without any signals of dissolution) along with the strong carbonate preservation potential of the eastern Mediterranean basin [72,93]. A total of 12 species of Euthecosomata (Heliconoides inflatus, Limacina trochiformis, Limacina bulimoides, Limacina retroversa, Creseis acicula, Creseis sp., Boasia chierchiae, Hyalocyclis striata, Styliola subula, Clio pyramidata s.l., Diacria trispinosa, Cavolinia spp.) were identified. Adult specimens, when present, were fragmented (Cavolinia spp., C. pyramidata). The protoconchs (Clio, Diacria and Cavolina) made the identification of certain species and genera possible, as they were the only residue left. The down-core variation of their abundance is presented in Figure 4.

Within the basal part of the core sequence, the fauna is composed almost exclusively of the pteropod L. retroversa. An additional component is the species C. pyramidata, but with very low percentages (<3.5%). Between 153 cm to 127 cm the pteropodal fauna becomes more diverse with the species H. inflatus, D. trispinosa, C. acicula, and B. chierchiae appearing in the fauna. The relative abundance of C. pyramidata gradually increases (up to 45%), in contrast to the decline of L. retroversa (drops to ~40%). The latter disappears completely from the fauna at 105 cm. Between 109 cm and 85 cm the pteropod fauna consists mainly of H. inflatus, C. acicula, C. pyramidata, and D. trispinosa, with L. bullimoides appearing for the first time at 109 cm. The species C. pyramidata and D. trispinosa are exponentially decreasing untill, and including, the top of the core sample. In the last 85 cm of the core, Cavolinia spp. appears in the fauna, reaching its maximum abundance (58%) at 41.25 cm and 85.5 cm, with the species H. inflatus and B. chierchiae (15–58% and 3–40% respectively) as additional components. At 65.5 cm, L. trochiformis presents a short occurrence and between 65 cm and 60 cm Creseis sp. presents its maximum relative abundance (36%). Towards the top of the core L. trochiformis and S. subula present their highest percentages (~6% and 15% respectively).

4.4. Total Organic Carbon and Stable Isotopes

Total organic carbon (TOC) in KIM-2A generally exhibits values around 1%. This pattern is interrupted in two intervals related to the sapropel sub-units deposition. The first interval (S1a; 89–64 cm) is characterized by high TOC concentration, ranging from 1.5% to 3.4%, and in the second interval (S1b; 52.5–40 cm), TOC concentration ranges from 1.8% to 2.4% (Figure 2a). Thus, the deposition of S1 sapropel layer started at 89 cm and terminated at 40 cm. The interruption of sapropel S1 layer (S1i) is detected between 64 cm and 52.5 cm, as suggested by the TOC content (~1.3%).

As shown in Figure 2a, between 196 cm and 159 cm, δ¹⁸O_{G. ruber} values range from +2.9 to +3.2‰. From 157 cm to 115 cm, a depletion in δ¹⁸O_{G. ruber} values is observed (up to +0.9‰) and persists until the 34 cm with even lower values. More precisely, the depleted values were recorded from two intervals (88.5–65.5 cm and 50.25–34 cm) with an average value of −0.08‰ and +0.25‰ respectively. In the interval corresponding to 66.5–61 cm, slightly heavier δ¹⁸O_{G. ruber} values were observed (+0.6 to −0.1‰). In the final unit of the core, from 28 cm until the top, a core enrichment in δ¹⁸O_{G. ruber} is recorded. The δ¹³C_{G. ruber} values of KIM-2A core exhibit more scatter than the δ¹⁸O records. In the basal part of the core (200–88 cm) values range between +1.3 to +0.2‰. In the interval between 88 cm and 61 cm the δ¹³C_{G. ruber} ranges between 0.3 and 0.9‰ with the exception of three high-positive peaks at 80.5 cm, 73.5 cm, and 71.5 cm with values of 1.4‰, 1.2‰, and 1.3‰ respectively. From 50.25 cm to the top of the core, δ¹³C_{G. ruber} exhibits heavier values with an average of +1.29‰ (Figure 2a).

4.5. Principal Component Analysis

A standardized principal component analysis (PCA) was carried out on the total data set using the varimax method, in order to determine the impact of various environmental parameters on the planktonic distribution. The application of this statistical analysis yielded a three-factor model for both planktonic foraminiferal and pteropod communities (Supplementary Materials; including PCA scores and biplots). The interpretation of the three components in each case was based on the screen plots of eigen values, and the factor loadings of the planktonic foraminiferal and pteropod species respectively. The 3 distinguished factors were considered to account for 81.57% and 82.81% of the total variance in each category respectively (

Table 2. Principle component analysis (PCA) factors based on planktonic foraminifera and their percentages of the total variability for core KIM-2A.

PCA Factors	Eigenvalue	% Variance	Cumulative % of the Total Variance
1	432.137	50.26	50.26
2	193.038	22.45	72.70
3	76.2131	8.86	81.57
4	53.3399	6.20	87.77
5	38.4841	4.47	92.25
6	24.4952	2.85	95.09
7	15.5611	1.81	96.90
8	7.97016	0.93	97.83
9	7.63218	0.89	98.72
10	5.03693	0.58	99.30
11	2.88765	0.33	99.64
12	2.25794	0.26	99.90
13	0.666282	0.08	99.98
14	0.159802	0.02	100.00

and

Table 3. PCA factors based on pteropods and their percentages of the total variability for core KIM-2A.

PCA	Eigenvalue	% Variance	Cumulative % of the Total Variance
1	1209.14	58.66	58.66
2	273.91	13.29	71.95
3	223.817	10.86	82.81
4	104.846	5.09	87.89
5	100.727	4.89	92.78
6	63.5077	3.08	95.87
7	53.2518	2.58	98.45
8	23.7844	1.15	99.60
9	4.32791	0.21	99.81
10	2.86233	0.14	99.95
11	0.755953	0.04	99.99
12	0.236728	0.01	100.00

), with their factor loadings showing the contribution of each factor in every sample and therefore the downcore contribution of each factor (

Table 4. Ranking of the planktonic foraminiferal species and their factor loadings along the PCA factors in core KIM-2A. Bold data indicate the most important factor loadings in each factor.

Variables	Factor 1	Factor 2	Factor 3
O. universa	−0.511	0.363	0.194
G. ruber f. alba	−0.158	−0.776	0.443
G. ruber f. rosea	−0.375	0.302	0.101
G. sacculifer gr.	−0.028	−0.091	−0.122
G. siphonifera gr.	−0.251	0.099	−0.152
G. inflata	0.079	−0.020	−0.136
G. bulloides gr.	0.094	−0.199	−0.691
G. rubescens	0.008	−0.012	−0.209
N. pachyderma	0.643	0.320	0.305
N. dutertrei	0.155	0.100	0.051
T. quinqueloba	0.197	0.001	0.258
G. truncatulinoides	0.007	−0.060	0.004
G. glutinata	0.129	−0.014	−0.110
G. scitula	0.032	−0.009	0.064

and

Table 5. Ranking of the pteropod species and their factor loadings along the PCA factors in core KIM-2A. Bold data indicate the most important factor loadings in each factor.

Species	Factor 1	Factor 2	Factor 3
H. inflatus	−0.387	0.693	0.453
L. bulimoides	−0.015	−0.012	0.010
L. retroversa	0.895	0.175	0.304
L. trochiformis	−0.005	0.000	0.008
B. chierchiae	−0.154	−0.334	0.286
C. acicula	−0.064	0.227	−0.328
Creseis sp.	−0.027	−0.031	−0.038
H. striata	−0.004	−0.002	0.008
S. subula	−0.006	0.015	0.004
C. pyramidata	0.009	−0.024	−0.437
D. trispinosa	0.051	−0.020	−0.420
Cavolinia sp.	−0.131	−0.570	0.380

; Figure 5). In the case of a bipolar factor, which has extremes of positive and negative loadings, high positive factor scores are related to the positive pole and high negative scores to the negative pole, respectively.

4.5.1. Planktonic Foraminifera

The first varimax factor (PCA1; Figure 5a) accounts for 50.26% of the total variance (Table 2) and is interpreted as temperature indicator. Species with positive loadings (N. pachyderma, T. quinqueloba, N. dutertrei, G. glutinata, G. inflata, G. bulloides, and G. scitula; Table 4) thrive in cold–water masses, while the species with negative loadings (O. universa, G. ruber f. rosea, G. siphonifera gr., and G. ruber f. alba; Table 4) thrive in warm–water conditions. Factor PCA2 (Figure 5b) accounts for 22.5% of the total variance (Table 1), with the positive pole being expressed by species living in a highly stratified water column (O. universa, G. ruber f. rosea and Neogloboquadrinids) and negative pole by species typical of the weak development of these conditions (G. ruber f. alba and G. bulloides gr.). The third varimax factor (PCA3) describes 8.9% of the total variance (Table 2), and also display a bipolar character, with its positive pole to be represented mainly by G. ruber f. alba and N. pachyderma and the negative pole by G. bulloides. The above species that characterize the PCA3 factor are the main exponents of the seasonal contrasts governing planktonic foraminiferal assemblages in the Mediterranean Sea during the last glacial cycle [94,95,96]. Thus, the third varimax factor (PCA3; Figure 5c) is referred to as a seasonality factor.

4.5.2. Pteropods

The first varimax factor (PCA1; Figure 5e) accounts for 58.66% of the total variance (Table 3) and was interpreted as a temperature indicator. Negative loadings consist mainly of the warm–water species H. inflatus, whereas positive loadings consist mainly of the subarctic species L. retroversa (Table 5). The second factor (PCA2; Figure 5e) describes 13.29% of the total variance (Table 3) and was interpreted as a productivity factor, as its positive pole is represented mainly by the mesopelagic oligotrophic H. inflatus. The third varimax factor (PCA3; Figure 5f) explains 10.27% of the total variance (Table 3). The positively loading taxa expressed by the epipelagic Cavolinia sp., L. retroversa and B. chierchiae, as well as the mesopelagic and tolerant to low oxygen concentration H. inflatus, whereas the negatively loading taxa (mesopelagic C. pyramidata and D. trispinosa) prefer a well-ventilated water column (Table 5). Thus, the third factor (PCA3) can be regarded as a stratification factor.

5. Discussion

5.1. Factors Controlling Planktonic Fauna Distribution in the Aegean Sea

Of the oceanographic factors typically considered, SST shows the highest explanatory power for the distribution of the planktonic fauna during the Late Quaternary. This is in accordance with previously published studies showing temperature as the dominant factor controlling the biogeography of planktonic foraminifera and pteropods at both global and local scales [61,67,97]. However, the 2 remaining factors (PCA-2, PCA-3) exhibit a bipolar character and could be considered as indicators of the annual stability of the water column. For the planktonic foraminifera fauna, they show that the faunal composition in the south Aegean Sea was not only controlled by SST, but also seems to be affected by the degree of development and location of a permanent or seasonal thermocline/pycnocline. Its vertical placement in the water column is a direct consequence of changes in sea surface salinity (SSS) and productivity (SSP), which ultimately reflected the seasonal fluctuations of the periods of vertical mixing of water in the periods of intense stratification. The interpretation of the second axis focused on the appearance depths of pycnocline and deep chlorophyll maximum (DCM) and the thickness of mixed layer. The interpretation of the third axis focused on upwelling currents and/or river inputs (e.g., G. bulloides), parameters which primarily control the food availability and reproductive cycles of foraminifera [73] and are directly correlated to the seasonal fluctuations they present [94,98,99,100]. Therefore, a useful additional dimension of planktonic foraminifera ecology that is underlined by the PCA conducted in this study is the degree of vertical stratification of the water column and the way it is recorded (seasonal presence/absence pycnocline and DCM and upwellings and runoff), which are inextricably linked with the factors of primary productivity and seasonality. Similarly, the second and third factor (PCA2, PCA3) of pteropods are focused on the hydrological conditions and the overall oxygenation of the water column. Particularly, the down-core scores of the second factor coincide with the δ¹³C and E-index values, indicating that variations in primary productivity have an impact on pteropod abundances. Even though nutrient concentrations are not a limited factor for their distribution [67], our data suggest that fluctuations in nutrients and salinity due to the increased freshwater inputs during the sapropel deposition favor the flourishment of some species (Cavolinia spp., B. chierchiae; Figure 4). Additionally, the third factor suggests that oxygen concentration, and thus the intensity of the oxygen minimum zone (OMZ), are parameters that affect pteropod distribution and particularly the mesopelagic species [66].

5.2. Paleoceanographic Reconstruction

The results of the multivariate statistical analyses, in combination with paleoceanographic indices and isotopic data (Figure 6), reveal a succession of Late glacial to Holocene paleoclimatic and paleoceanographic changes. The evidence of these changes is interpreted and discussed in terms of the events that mainly accompanied the transition out of the late glacial period and the deposition of sapropel S1 during the Holocene Climatic Optimum (HCO).

5.2.1. Late Glacial

During the late glacial period (21.1–15.7 ka BP), the heaviest δ¹⁸O values (2.49–3.26‰), accompanied by relatively low PPC values (−32% to +4%), suggest a cold upper water column (Figure 6a,b). Particularly, this interval was characterized by high percentages of the cold water foraminifera species T. quinqueloba (~30%), accompanied by G. glutinata (9%), and G. scitula (8%), and significant percentages of the warm-water G. ruber f. alba (~34%), that are suggestive of milder climate after the Last Glacial Maximum (LGM), which is in accordance with other records in the Mediterranean [101,102]. Pteropod fauna is composed mainly by the cold water species L. retroversa (98%) and the cold-tolerant mesopelagic C. pyramidata in very low percentages (3%) which is consistent with relevant late glacial Mediterranean records [33,61,103]. Eutrophic species are also abundant in this interval (N. pahyderma, N. dutertrei, T. quinqueloba and G. bulloides) and are associated with the high values of E-index (Figure 6e). Notably, the high abundance of N. pachyderma (26%) indicates shallowing of the pycnocline and the formation of a DCM layer. In addition, δ¹³C values, around +1‰, and the trend to higher S-index values (Figure 6f,j), also support the development of eutrophicated waters. The moderate abundance of G. bulloides (5−17%; Figure 3) suggests little to no upwelling and/or runoff contribution in primary productivity. Therefore, the injection of nutrients into the euphotic zone can be attributed to the intensification and southward shift of westerly winds, as indicated by atmospheric circulation models for this time interval [104].

5.2.2. Deglaciation

At 15.7 ka BP an abrupt shift of PPC to positive values, accompanied by lighter δ¹⁸O values (+2.2‰) (Figure 6a,b), are indicative of the climatic amelioration that occurred during the last deglaciation. These warmer conditions are also supported by the occurrence of the warm water species H. inflatus and D. trispinosa (up to 60% and 50% respectively), and the decrease in L. retroversa percentages (between 33% and 85%). This warming trend is in agreement with relevant paleoclimatic records from the eastern Mediterranean, and is attributed to the Bølling–Allerød (B-A) interstadial [7,13,14,105,106]. The increased SST and humidity are also recorded by the higher abundance of the terrestrial biomarkers in the south Aegean [17], and by a change in the benthic faunas from oxic to dysoxic indicator species [13]. In the beginning of this interval, Neogloboquadrinids were temporarily replaced by G. bulloides (26%), suggesting local upwelling. Though, later on the eutrophic N. pachyderma and N. dutertrei present their highest abundance (42% and 14% respectively). Additional components of this interval are G. ruber (both variants), G. bulloides, T. quinqueloba, and G. inflata, suggesting temperate and meso- to eutrophic waters, with strong seasonal mixing and local upwelling.

This state persisted until the onset of the Younger Dryas (YD) at about 12.9 ka BP, which is depicted in the abrupt decrease in PPC (from +28% to −10%) and in heavier values of δ¹⁸O (1.0–2.5‰) (Figure 6a,b). The planktonic foraminiferal fauna shows an increase in cold water species (N. pachyderma, T. quinqueloba, G. inflata and G. glutinata). Pteropod fauna is composed mainly of the cold-water species epipelagic L. retroversa and the mesopelagic temperate to warm-water species C. pyramidata and D. trispinosa, while the warm-water H. inflatus presents a decreasing trend (Figure 3). This climate response of south Aegean depression to the YD event (12.9–11.7 ka BP) seems to be in accordance with relevant signals from other Aegean sub-basins [6,7]. Towards the end of YD (12.6–12.2 ka BP) the species G. rubescens, C. acicula, and B. chierchiae are added to the fauna suggesting that mild climatic conditions prevailed for a brief time interval of about 400 years within the YD event. This brief climatic amelioration in the mid YD, which has been attributed to the displacement of the polar front by a few degrees north [26], has been also observed in north-central Aegean marine records [7,106] and coincides with the pattern of GRIP and NGRIP ice-core records [107,108], pollen-based reconstructions from the Jura [109] and the Balkans [110], as well as chironomid-based reconstructions from North Italy [111]. The increases in the S-index and E-index at around 12.2 ka BP (Figure 6e,j) coincide with this amelioration, reflected by the improved ventilation and eutrophication of the water column.

5.2.3. Holocene

With the ending of the YD, a general climatic amelioration is seen in the records (increase in PPC, lighter values in δ¹⁸O, decline of temperature factors PCA1; Figure 6a–d) marking the beginning of the Holocene. Planktonic foraminifera fauna consists mainly by G. ruber f. alba that increases in abundance towards the onset of sapropel deposition, along with G. bulloides and N. pachyderma. The two latter species present an opposite trend, decreasing towards the onset of sapropel deposition. This trend suggests the gradual development of stratified and oligotrophic surface waters. Pteropod fauna follows a similar pattern with the decreasing abundances of the species D. trispinosa and C. pyramidata, which are indicative of a well ventilated water column [112], and the increase in the epipelagic B. chierchiae and C. acicula (Figure 4). The reduction of E-index along with the values of δ¹³C and the S-index (Figure 6e–j) are further evidence for the development of these conditions. The replacement of G. inflata by G. glutinata (Figure 3) may be tentatively explained in terms of increased seasonality, and of increased freshwater input, which reduced surface buoyancy loss and hence suppressed mixing during the beginning of the Holocene [7]. This was also enhanced by the presence of the epipelagic pteropods (Creseis spp.) that proliferate in low-salinity waters [113], and the seasonality factor (PCA3 of planktonic foraminifera; Figure 6k).

At around 9.4 ka BP, the deposition of sapropel S1 as witnessed by their high organic carbon content (C_org: 1.8–3.3%, Figure 2a), coincided with the start of the overall δ¹⁸O depletion (+0.6‰ to −0.9‰; Figure 6b). The changes of the planktonic foraminifera fauna are characterized by an increase in G. ruber f. rosea, G. ruber f. alba, G. siphonifera and O. universa that are suggestive of extremely warm and stratified conditions (PPC ~100%; PCA1 low values; Figure 6b–d). The lower values of G. bulloides/G. ruber ratio in this interval are also indicative of strongly stratified water column and are in agreement with the stratification factor for both faunas (PCA2 of planktonic foraminifera and PCA3 of pteropod; Figure 6h–j). An increase in temperature and humidity around this time has been documented in all marine and terrestrial pollen records in the eastern Mediterranean region [6,7,14,15,16,114]. This paleoclimate change coincides with the Holocene summer precession-related insolation maximum in the Northern Hemisphere [115], and the monsoon intensification that resulted in a widespread increase in humidity over the Mediterranean region and concomitant increase of freshwater input to the Mediterranean Sea [116,117]. Pteropod fauna is characterized by the dominance of the warm oligotrophic H. inflatus, and the warm epipelagic B. chierchiae, C. acicula, and Cavolinia spp. (Figure 4). Mesopelagic species (C. pyramidata and D. trispinosa) are decreasing dramatically due to the enhanced stratification of the entire water column. Mesopelagic pteropods are affected by the OMZ alterations, which are climatically controlled [63,112]. In the humid and warmer conditions that persisted during the formation of S1, the subsequent stratification of the water column favored a strong and well developed OMZ that probably led to the reduction of mesopelagic species. The presence of the mesopelagic H. inflatus into the sapropel sublayers can be explained by its habitat. More explicitly, this species adopts a variable depth habitat during its growth stages, and it is more susceptible to the low oxygen concentration in the OMZ [63,118,119]. The presence of L. bulimoides in the upper part of S1a, with peaks at 8.4 ka BP, 8.1 ka BP, and 7.9 ka BP (~6%) and its absence in the S1b, suggest that during the end of S1a (8.6–7.7 ka BP) the conditions were more arid than during the onset of S1a (9.4–8.6 ka BP) and the interval of S1b (6.9–6.4 ka BP). In these two phases of S1a, δ¹³C present a decreasing trend with an average value of +0.8‰ in the first phase and +0.6‰ in the second (Figure 6f).

The warm and stratified conditions favorable for the sapropel deposition were interrupted between 7.7 ka BP and 6.8 ka BP. This interval (S1i) is marked by the decrease in PPC (from 90% to 60%) and the heavier values of δ¹⁸O (+0.5‰) and δ¹³C (+0.8‰) as shown in Figure 6a,b,f. The subsequent cooling is also reflected in significant faunal changes, such as the increase in abundance of G. inflata, T. quinqueloba, G. bulloides, and N. pachyderma (Figure 3). These species are associated with relatively cold temperatures and increased food supply, suggesting high primary production and stronger mixing of the water column [46,52]. In addition, the pteropod L. trochiformis, which is related to the mixed layer of the water column and thrives in upwelling conditions [120,121,122,123], presents a peak at the beginning of S1i (Figure 4).

From 6.4 ka BP to the top of the core (~5.0 ka BP), a trend to heavier δ¹⁸O values (from 0‰ to +1.3‰) and the drop of PPC (~60%) are recorded (Figure 6a,b). In planktonic foraminifera fauna an increase in G. inflata, G. truncatulinoides and N. pachyderma along with the reduction of G. ruber f. rosea, G. siphonifera gr. O. universa indicate lower SST, and stronger seafloor oxygenation due to vertical mixing. This latter is also suggested by the increase in the mesopelagic pteropod D. trispinosa. The peak of L. trochiformis at 5.7 ka BP (Figure 4) is suggestive of upwelling conditions. At this point, and up to 5.0 ka BP, increased percentages of G. bulloides and G. sacculifer are indicative of increased productivity. This is also indicated by the high values of E-index, the heavy δ¹³C (up to +1.7‰) and the PCA2 factor of pteropods (Figure 6e–g).

6. Conclusions

The results of the analysis of planktonic foraminifera and pteropods, combined with the results of principal component analysis and oxygen and carbon isotopic signals, describe the key factors controlling the formation of the micropaleontological assemblages, and therefore the paleoenvironmental and paleoclimatic changes during the last ~21 ka BP. Our data suggests that the sea surface temperature, stratification of the water column, seasonality, and productivity are the main controlling factors of the faunal distribution. During the Late Glacial, all the records indicate the occurrence of cold and eutrophicated waters. The time interval of 15.9–11.7 ka BP, corresponding to the deglaciation phase, was characterized by gradual climatic amelioration. During this interval, the warm interstadial Bølling–Allerød was identified, while low SST records at 12.9 ka BP revealed the onset of the YD. During this event, mild climate conditions were observed for around 400 years, and are attributed to the displacement of the polar front by a few degrees north. With the onset of the Holocene, a general climatic amelioration can be observed, with the gradual development of stratified oligotrophic surface waters. These conditions were intensified during the sapropel S1 deposition, which appears in two layers (S1a and S1b). The faunal and isotopic data suggest that conditions were more arid towards the end of S1a than at the onset of the S1a and the S1b. The interruption of sapropel deposition (S1i) and the post-sapropel interval are characterized by lower SST and stronger seafloor oxygenation due to vertical mixing.