DEVELOPMENT OF COCCOLITHOPHORE-BASED TRANSFER FUNCTIONS IN THE WESTERN MEDITERRANEAN SEA: A SEA SURFACE SALINITY RECONSTRUCTION FOR THE LAST 15.5 KYR

DEVELOPMENT OF COCCOLITHOPHORE-BASED TRANSFER FUNCTIONS IN THE WESTERN MEDITERRANEAN SEA: A SEA SURFACE SALINITY RECONSTRUCTION FOR THE LAST 15.5 KYR B. Ausín (b_ausin@usal.es); I. Hernández-Almeida (ivan.hernandez@giub.unibe.ch); J-A. Flores (flores@usal.es); F-J. Sierro (sierro@usal.es); M. Grosjean (martin.grosjean@oeschger.unibe.ch); G. Francés (gfrances@uvigo.es); B. Alonso (belen@icm.csic.es).


Introduction
Coccolithophores are one of the major components of marine phytoplankton.They are sensitive to changes in many environmental variables, such as nutrients, temperature and salinity, and are widely used in qualitative paleoenvironmental studies (Baumann et al., 2005;Guerreiro et al., 2013Guerreiro et al., , 2014)).These studies provide general insight into the response of coccolithophores to environmental variables, but quantitative studies (e.g.transfer functions) allow assessing these relationships in a more rigorous and clear manner.Transfer functions are based on the calibration of the modern relationship between organisms and environmental conditions, and this information is in turn used to reconstruct past environmental variables.Different statistical approaches based on coccolithophores have been proposed in order to generate quantitative paleoreconstructions of different ecological variables.Giraudeau and Rogers (1994) used factor analyses and multiple regressions to estimate chlorophyll a from coccolithophore census counts in surface sediment samples Published by Copernicus Publications on behalf of the European Geosciences Union.
in the Benguela upwelling area.Several authors (Beaufort et al., 1997(Beaufort et al., , 2001;;Incarbona et al., 2008) have calibrated the relative abundance of the coccolithophore Florisphaera profunda in surface sediment samples with respect to primary productivity and reconstructed past variations of this parameter in the Indian and Pacific oceans and in the central Mediterranean Sea.Saavedra-Pellitero et al. (2011, 2013) used linear regression methods to derive past sea surface temperature (SST) estimates in the southeast Pacific Ocean from coccolithophore census counts and accumulation rates.Bollmann et al. (2009) and Bollmann and Herrle (2007) applied multiple linear regressions to morphometric measurements of the coccolithophore Emiliania huxleyi from globally distributed core-top and plankton samples to obtain modern and past sea surface salinity (SSS) estimates.
To date, no coccolithophore-based transfer function has been applied in the western Mediterranean Sea, a semienclosed basin situated at midlatitudes (Fig. 1a).In this region evaporation exceeds precipitation plus runoff, such that water budgets tend to be balanced by the advection of relatively less saline Atlantic water (AW) through the Strait of Gibraltar (Bèthoux, 1979).The AW flows eastward while mixing with Mediterranean water to form the Modified Atlantic Water (MAW) at the surface (100-200 m) (Millot, 1999).This distinctive feature affects the spatial distribution of some environmental parameters such as SST and SSS, leading to the development of well-defined longitudinal gradients between the Atlantic Ocean and the western Mediterranean in annual terms.In this confined basin, the estimation of changes in these environmental parameters is essential for determining Atlantic-Mediterranean water mass exchange through the Strait of Gibraltar in the past (Rohling and Bigg, 1998;Schmidt, 1998).This exchange depends on variations in the hydrological cycle, ice-volume effects, and Mediterranean circulation patterns, which have a thermohaline origin (MEDOCGROUP, 1970).
The aim of this study is to explore the potential of coccolithophores for the development of quantitative reconstructions in the western Mediterranean Sea.We study the response of coccolithophore assemblages from surface sediment samples from the Atlantic Ocean and Mediterranean Sea to environmental variables.The resulting calibration model (transfer function) for salinity was used to reconstruct SSS changes at a high resolution in the Alboran Sea (Fig. 1a) for the last 25 kyr.The reliability of the reconstruction was assessed by analysis of the similarity between fossil and modern coccolithophore assemblages and fossil ordination scores.Finally, centennial and millennial SSS changes are described, discussed, and compared with regional records of SST and organic-matter preservation.

Surface sediment samples
Initially, 117 core tops located around a horizontal transect along the western Mediterranean Sea and near the Gulf of Cadiz in the Atlantic Ocean were selected.They had been retrieved at varying water depths ranging from 70 to 2620 m during several oceanographic surveys and were stored at the University of Vigo and at the Core Repository of the Institute of Marine Sciences (CSIC) in Barcelona.The first centimeter (or the second if the first was unavailable) of the 117 core tops was sampled, assuming that it essentially represents present-day conditions.

Fossil data set
The fossil coccolithophore data set used for the reconstruction comprises coccolithophore census counts from core CEUTA10PC08 (36 • 1 22 N, 4 • 52 3 W; 914 m b.s.l.), located in the Alboran Sea, previously published by Ausín et al. (2015).Fossil assemblages show a good to moderate degree of preservation.The location of this core lies under the modern path of the AW at the surface, near the Strait of Gibraltar (Fig. 1a).Sediment core chronostratigraphy was based on 15 radiocarbon ages and covered the time span from 25 to 4.5 ka calibrated BP at a ∼ 65 years temporal resolution (Ausín et al., 2015).All dates reported in this study are given in calibrated ages BP.

Micropaleontological analyses
Both modern (surface sediment) and fossil (down-core) samples were prepared for coccolithophore analyses according to the techniques proposed by Flores and Sierro (1997).A polarized-light microscope at 1000× magnification was em-  (Schlitzer, 2014).ployed to identify and count at least 500 coccoliths in each sample, belonging to 21 different taxa.Species whose relative abundance was < 1 % in the first count were considered later in 20 visual fields in order to estimate their abundance accurately.The final relative abundance of each species in each sample was then recalculated.Gephyrocapsa specimens smaller than 3 µm were lumped together and designated "small Gephyrocapsa" (Flores et al., 1997).The "medium Gephyrocapsa" group was made up of Gephyrocapsa whose size was between 3 and 5 µm.Two sizes of morphotypes of E. huxleyi (< 4 and > 4 µm) were considered, owing to their different ecological and biostratigraphic significance in the study area.Similarly, G. oceanica was split according to a size criterion of < 5 and > 5 µm, owing to their comparable potential ecological significance.Other taxa identified in this study were Calcidiscus leptoporus, F. profunda, Gephyrocapsa cf.caribbeanica, Gephyrocapsa muellerae, Helicosphaera spp., and Syracosphaera spp.(as dominant taxa).The rare taxa identified were Braarudosphaera bigelowii, Calciosolenia spp., Coccolithus pelagicus subsp.braarudii, Coccolithus pelagicus subsp.pelagicus, Oolithotus fragilis, Pontosphaera spp., Rhabdosphaera clavigera, Umbilicosphaera spp.and Umbellosphaera spp.
Twenty-nine samples were finally eliminated from the initial modern data set owing to their high content (> 10 %) in obviously reworked nannofossils.These taxa belong to older stratigraphic levels (consistently older than the Pliocene in this study), meaning that they were resuspended and transported from their original location to the sample site, and they lack any relationship with modern environmental conditions.A maximum of 10 % of reworked specimens was chosen as an acceptable threshold below which the sample could be retained in the modern training set, after ruling out these reworked specimens, without compromising the statistical representativeness of the major species (Fatela and Taborda, 2002).Later examination of the spatial distribution of reworked specimens in the retained samples revealed that those with the highest percentages were close to river mouths, relating reworked specimens to river discharges and suggesting that the rest of the assemblage could be considered autochthonous.Thus, the final training set (Supplement) comprised 88 surface samples (Fig. 1b): 78 from the western Mediterranean (58 from the Balearic Sea and 20 from the Alboran Sea) and 10 from the Atlantic Ocean.Coccolithophore relative abundances were square-roottransformed to stabilize their variances.The species Braarudosphaera sp., Calciosolenia spp., Coccolithus pelagicus subsp.braarudii, Coccolithus pelagicus subsp.pelagicus, and Pontosphaera spp.were excluded from the modern (and consequently from the fossil) coccolithophore assemblages since their maximum relative abundance was not > 1 % in at least two samples.Detrended correspondence analysis (DCA) was then performed on the modern coccolithophore assemblage to estimate the length of the environmental gradient.A length of the first DCA axis > 2 standard deviations (SD) indicates the unimodal responses of the organisms (Birks, 1995;ter Braak and Prentice, 1988), while shorter lengths indicate linear responses.
Akaike's information criterion (AIC) was used in an ordination analysis to identify the minimum number of variables (subset) that, being statistically significant, explained the maximum variation in the modern coccolithophore assemblage.Canonical correspondence analysis (CCA) was used to evaluate the influence of this environmental subset to explain coccolithophore distribution in the modern training set.
The ratio between the first constrained axis and the first unconstrained axis (λ 1 /λ 2 ) was used as a diagnosis to test the strength of a single environmental variable when the effects of those remaining were excluded from the analyses (ter Braak and Juggins, 1993).If λ 1 /λ 2 ≥ 1, the variable under examination is considered to be important for explaining the distribution of the species.The proportion of the variance in the coccolithophore training set explained uniquely by each significant environmental variable was calculated using variance partitioning.

Transfer function
Calibration models were calculated for the variable of interest (and each variable by means of exploratory analysis) using the weighted-averaging partial least squares (WA-PLS) method (ter Braak and Juggins, 1993;ter Braak et al., 1993) and the modern analog technique (MAT) (Prell, 1985), both implemented in the C2 version 1.4.3 software (Juggins, 2007).All models were calculated for the crossvalidation predictions by bootstrapping (999 permutation cycles) (Birks, 1995).In the AT, the number of analogs resulting in the maximum coefficient of determination (R 2 boot ) between the observed and predicted values and the lowest root-mean-square error of prediction (RMSEP) (Telford et al., 2004) was calculated using an optimization set together with the usual training and test sets implemented in the "analog" package for R (R Core Team, 2015).In WA-PLS, a decrease of 5 % or more in RMSEP was required to retain the next component (Birks, 1995;ter Braak et al., 1993).
Many coccolithophore species inhabit depths within a specific range of the photic zone and are subject to environmental seasonality (Winter et al., 1994).Therefore, the depth and season considered for calibration and reconstruction should be those that most influenced the coccolithophore fossil assemblage.Following the procedure described by Telford et al. (2013), we reconstructed the variable of interest based on summer-, winter-and annually averaged data at nine different depths of the photic zone from 10 to 200 m using the "pale-oSig" package v.1.1-1(Telford, 2012) for R (R Core Team, 2015).The reconstruction that explains the highest proportion statistically significant of variance in the fossil data reflects the depth and season that most influenced the coccolithophore fossil assemblage, and it hence provides the most suitable calibration.
Outliers may reduce the power of prediction of the calibration model as well as introduce undesirable effects in model coefficients (Birks, 1995).Potential outliers were determined as those whose absolute residual was higher than the mean SD of the observed values (Edwards et al., 2004).
A combination of the highest R 2 boot and the lowest RMSEP was used as a criterion for the quality prediction of the model.The graphical representations of the observed values against the values predicted by the model and the residuals against the predicted values were used as a diagnosis of the model.

Derived reconstruction and evaluation
The MAT and WA-PLS were applied to the fossil coccolithophore assemblages of core CEUTA10PC08, which were previously square-root-transformed, using the C2 version 1.4.3 software (Juggins, 2007).Sample-specific reconstruction errors under bootstrapping were derived automatically by the C2 software, considering the prediction error due to (i) errors in estimating species coefficients and (ii) errors in the calibration function (further details may be found in Birks et al., 1990).In order to assess the quality of the modern analogs for the fossil (down-core) samples, the squared chord distance between each fossil sample and each sample in the modern training set (Overpeck et al., 1985) was calculated with the MAT by the C2 software.A squared chord distance below the 10th percentile would be considered good, while values above this cutoff would represent assemblages with poor analogs (Simpson, 2007).
The first axis of the PCAs of the fossil data set (PC1 fossil ) shows the most important changes in the composition of the fossil coccolithophore assemblage.Comparison between PC1 fossil and the reconstructed variable of interest was used to assess whether the reconstruction could be considered representative of the major ecological changes in the fossil assemblage (Juggins, 2013).

Geographical distribution of coccolithophores
The small placoliths (small Gephyrocapsa and E. huxleyi < 4 µm) are the dominant taxa (Fig. 2b, c) constituting on average 83 % of coccolithophore assemblages.Small examples of Gephyrocapsa show higher abundances near the Spanish coast and southeast of the Balearic Islands.E. huxleyi < 4 µm is more abundant in the Balearic Sea and around the Ebro River delta (Fig. 2c).G. muellerae (Fig. 2d) is concentrated southeast of the Balearic Islands and a patch of 2 % appears in the northern Alboran Sea. C. leptoporus and Helicosphaera spp.(Fig. 2e, f) are almost absent in the Alboran Sea and show similar patchy distributions between the Catalan and the Balearic fronts and east of the Balearic Islands.F. profunda (Fig. 2g) is more abundant in the Atlantic Ocean (up to 16 %) and gradually decreases eastward.It appears in two patches (up to 4 %) south of the Ebro River mouth.G. oceanica (< 5 µm) (Fig. 2h) is mostly distributed near the Strait of Gibraltar.It also appears in a patch (up to 3 %) around the Andarax River mouth.

Relationship between coccolithophores and environmental variables
The PC1 explains 56.1 % of the variance within the environmental data set (Fig. 3a) and is highly correlated with CO 2− 3 , salinity, pH and T ALK .PC2 explains 22.3 % of the total variance and primarily summarizes the information on temperature and phosphate.
The ordination based on the AIC revealed that only salinity, nitrate, phosphate, silicate and oxygen are needed to explain the maximum variation in the modern coccolithophore assemblage and are significant at the 95 % level when added individually to the model via a forward selection procedure.The first axis of the DCA performed on the modern coccolithophore assemblage was 2.6 SD.Accordingly, unimodal methods were followed.The CCA (Fig. 3b) revealed sites and species distribution along this environmental subset.The others were also plotted as passive variables to avoid overfitting.The vectors show that salinity exhibits the longest gradient and is strongly correlated with the CCA1, indicating a strong relationship with coccolithophore distribution.Some sites from the Alboran and Balearic Seas and the taxa medium Gephyrocapsa and E. huxleyi (> 4 µm) were found to be distributed along the CCA2.Individual CCAs (Table 1) to calculate λ 1 /λ 2 showed that salinity was the most important variable among those found to be significant.Variance partitioning revealed that these significant variables accounted for 38.9 % of the cumulative variance in the coccolithophore training set and salinity explained a large proportion of this variance (15.5 %).

Transfer functions
Salinity explained the largest amount of variation in the coccolithophore assemblages and was therefore chosen to develop the coccolithophore-based transfer function.Additionally, comparison among the R 2 boot from preliminary calibration models for each variable confirmed the best predictive power for salinity (Table 1).
Among the WA-PLS models for salinity, the twocomponent model (WA-PLS2) was chosen as the most suitable since it afforded a reduction of 6.4 % in the RMSEP.The ideal number of analogs for the MAT was six.
The analyses of the amount of down-core variance explained by the summer, winter, and annual salinity reconstructions at nine different depths and their statistical significance revealed that the mean annual reconstruction at 10 m explained the highest variance.Hence, the reconstruction for core CEUTA10PC08 was based on the mean annual salinity data at 10 m depth and referred to as SSS reconstruction.
Five samples showed higher residuals than the SD of salinity and were preliminary identified as potential outliers (Supplement).However, only one of these samples (CO-81-2/TK-2) was identified as an outlier in both the MAT and WA-PLS regression methods.This had a bright yellowish color under the microscope, likely due to the effect of diagenetic processes.In order to retain the maximum number of observations representing modern environmental conditions, only this sample was removed from subsequent model implementations, leading to an improvement of the MAT and WA-PLS2 R 2 boot coefficient of 3.4 and 6.6 %, respectively, and reducing both Max_Bias boot and RMSEP (Table 2).
The final MAT and WA-PLS2 models showed similar quality predictions (Table 2).The salinity values in the modern training set vary from 36.2 to 38.2 psu.Intermediate values (37.1-37.6 psu) are less well represented by the observations (Fig. 4a).MAT-and WA-PLS2-predicted values are shown in Fig. 4b, c.The predicted vs. observed values from both models approach the diagonal of slope one (which indicates perfect predictions) reasonably well (Fig. 4d, e).The residuals for the MAT and WA-PLS2 models (Fig. 4f, g) are equally distributed around zero and show no apparent trends.

SSS reconstruction
SSS trends and values reconstructed for the CEUTA10PC08 core derived from both the MAT and WA-PLS2 are very similar (Fig. 5a, b).These only differ during the stadials associated with Heinrich events 2 and 1 (H2 and H1), when the WA-PLS2-estimated SSS shows more pronounced salinity decreases.
The SSS reconstructions obtained from core CEUTA10PC08 (Fig. 5a) can be divided into three intervals: i.The period from 25.5 to 15.5 ka is characterized by higher values that oscillate between 37.8 and 37 psu.
Lower values are found from 20 to 18 ka, followed by a drop of 0.8 psu at 17.3 ka.
ii.The period from 15.5 to 9 ka shows fast, largeamplitude changes.An abrupt decrease from 37.9 to 36.9 psu can be recognized at 15 ka, followed by large peaks of high values at 12.8, 11.1, and 10.2 ka.iii.The period from 9 to 4.5 ka records the lowest values, which vary between 37 and 36.5 psu, and shows a general decreasing trend.
On average, the errors associated with both SSS reconstructions are of a similar magnitude: ±0.15 psu for the MAT and ±0.17 psu for WA-PLS (Fig. 5a).Squared chord distances between fossil and modern assemblages (Fig. 5b) revealed that many samples from 25.5 to 16 ka were above the 10th percentile.A comparison between PC1 fossil and the SSS reconstruction is depicted in Fig. 5c, showing generally good agreement, especially for the last 16 kyr.

Discussion
4.1 Geographic coccolithophore distribution and SSS E. huxleyi (< 4 µm) and small Gephyrocapsa are widespread in the western Mediterranean, as previously reported for surface sediment and water column samples (Álvarez et al.,   Knappertsbusch, 1993;Oviedo et al., 2015).These taxa, especially E. huxleyi (< 4 µm), are cosmopolitan and tolerate wide ranges of temperature and salinity (Winter et al., 1994).G. muellerae abundance is higher southeast of the Balearic Islands, where the MAW encounters more saline and warmer Mediterranean waters, and close to the Alboran Front, possibly reflecting the species' preference for nutrient-rich waters, as reported for sediment trap samples in the Alboran Sea (Bárcena et al., 2004;Hernández-Almeida et al., 2011).C. leptoporus and Helicosphaera spp.(Fig. 2e, f) show similar spatial distributions and abundances.Interestingly, the CCA suggests that Helicosphaera spp.have a preference for more saline waters (Fig. 3b).By contrast, in paleoceanographic works this species has been linked to fresher and turbid waters in the Mediterranean Sea (Ausín et al., 2015;Colmenero-Hidalgo et al., 2004;Grelaud et al., 2012).Helicosphaera spp.abundance in surface sediments from the northeastern Balearic Island has also been related to upwelling events (Álvarez et al., 2010).Similarly, the abundance of C. leptoporus in the Alboran Sea has been linked to nutrient-rich waters (Bárcena et al., 2004).The similar patchy pattern shown by both species may be related to the temporary upwelling of nutrient-rich waters associated with frontal structures in the area limited by the Balearic and Catalan fronts (Font et al., 1988).In agreement with this interpretation, the co-occurrence of both species in other Mediterranean locations has already been linked to highly productive coccolithophore periods and pre-upwelling events (Hernández-Almeida et al., 2011;Ziveri et al., 2000).The CCA (Fig. 3b) suggests that F. profunda and G. oceanica (< 5 µm) would be associated to less saline waters.This notion may partly be a consequence of their higher abundance in Atlantic waters (Fig. 2d).G. oceanica has already been proposed as a tracer for AW influx in the western Mediterranean Sea (Álvarez et al., 2010;Bárcena et al., 2004;Knappertsbusch, 1993;Oviedo et al., 2015).Similarly, the F. profunda spatial distribution reflects the path of the Algerian current (Fig. 2a, d), formed by recent and fresher MAW (Fig. 2a,  d).Low percentages of F. profunda patchily distributed south of the Ebro River and in the Catalan-Balearic Sea suggest that this species may also be affected by the influence of river discharges (Álvarez et al., 2010).These results suggest that F. profunda and G. oceanica proliferate mainly in waters of Atlantic origin but not exclusively, as indicated by their presence in the eastern Mediterranean (Knappertsbusch, 1993;Malinverno et al., 2008) where Atlantic influence becomes diluted.
Salinity was highly correlated with CO 2− 3 and pH (Fig. 3a,  b).From the study of coccolithophore distributions in water column samples and in situ environmental measurements in the Mediterranean Sea, Oviedo et al. (2015) found exactly the same variables as being the most important factors when accounting for changes in heterococcolithophore assemblages.In our study, multivariate analyses revealed that salinity was significant and was the most important variable of those studied to explain the variance in coccolithophore data in this modern training set.However, the individual importance and proportion of variance explained by each of the significant variables was not assessed in the study of Oviedo et al. (2015).Despite this, the authors discarded salinity as a final explanatory variable, arguing that E. huxleyi, the most abundant and ubiquitous extant coccolithophore (Cros and Fortuño, 2002), inhabits a wide salinity range, suggesting a negligible ecological effect of salinity on coccolithophores.Contrary to this reasoning, the direct relationship between varying salinities and the morphology of E. huxleyi has been demonstrated by several authors (Bollmann and Herrle, 2007;Bollmann et al., 2009;Fielding et al., 2009;Green et al., 1998;Paasche et al., 1996;Schouten et al., 2006) in both culture experiments and marine surface sediment samples.Oviedo et al. (2015) later explained the strong and negative relationship that they found between salinity and G. oceanica, G. muellerae and E. huxleyi morphotype B/C distributions as being a consequence of their carryover by the AW through the Mediterranean.Instead of this, however, we interpret the AW influx as promoting the optimal conditions for these species to thrive in the Mediterranean Sea.Therefore, the coccolithophore relationship with salinity would reflect the different water masses which coccolithophore species prefer to inhabit.
It is worth mentioning that salinity influences the solubility of CO 2− 3 via several pathways: the solubility of free carbon dioxide in water, the solubility product constants, the concentration of hydrogen ions, and the quantity of calcium in the water (Trask, 1936).Accordingly, salinity could influence coccolithophores through coccolith calcification processes.In contrast, Bollmann and Herrle (2009) have proposed that salinity influences coccolithophores through cell turgor regulation linked to osmotic processes.Although there is no clear consensus about the mechanism through which salinity influences coccolithophores, many other studies point to a strong influence of this variable on molecular compounds only produced by coccolithophores and on specific species.In the Japan Sea, salinity has been proposed to have an ecological or physiological influence on the production of alkenone and alkenoates, which are organic compounds mainly produced by the genera Emiliania and Gephyrocapsa (Fujine et al., 2006).In the Baltic Sea, alkenone unsaturation ratios have been found to be significantly correlated with salinity (Blanz et al., 2005).In the Mediterranean Sea, Knappertsbusch (1993) found that G. oceanica distribution was linearly correlated with salinity.Based on such evidence, we propose that the assemblage composition may be conditioned by the optimum salinity range preferred by each species.Moreover, salinity has proved to be important to other marine unicellular planktonic groups such as diatoms (Jiang et al., 2014;Li et al., 2012) and dinoflagellate cysts (Jansson et al., 2014, and references therein), reinforcing the hypothesis of salinity as an important variable for planktonic communities in semi-enclosed basins.

Transfer function quality
A general good fit can be deduced for both models, although the MAT was seen to perform slightly better from a higher R 2 boot , a lower RMSEP (Table 2), and the comparison predicted values compared with observed values (Fig. 4).Intermediate salinity values (37.1-37.6 psu) are less well represented than the more extreme values (Fig. 4d, e).Unevenness can bias the RMSEP leading to overestimation of the predictive power of the model (Telford and Birks, 2011).While an even distribution would always be desirable, unevenness is a feature inherent to most training sets from oceanic environments.In this case, it is not severe, and the observations, although unevenly distributed along the salinity gradient, do not leave gaps.The distribution of the residuals (Fig. 4f, g) indicates the adequacy of the model.

Down-core SSS reconstruction
The derived MAT and WA-PLS2 SSS reconstructions (Fig. 5a) are very similar.Nevertheless, WA-PLS2 shows more marked salinity decreases than the MAT during the H2 (25.2-23.7 ka) and H1 (17.4-15.9ka).Unlike WA-PLS, the MAT does not consider the entire data set when cal-culating the species optima, only the most taxonomically similar analogs, and is more sensitive to local conditions (Telford and Birks, 2009).Fossil samples lack good analogs for the H2 and H1, coinciding with large peaks of E. huxleyi (> 4 µm) (Fig. 5b).H2 and H1 have been linked to the entry of cold and fresher water originating from the North Atlantic ice melting in the western Mediterranean Sea (Cacho et al., 1999;Melki, 2011;Sierro et al., 2005), suggesting the preference of E. huxleyi (> 4 µm) not only for cold waters (Colmenero-Hidalgo et al., 2002, 2004) but also fresher waters in the past.By contrast, Bollmann and Herrle (2007) reported a current positive correlation between the size of E. huxleyi up to 4 µm and increasing salinities from the study of globally distributed core-top samples.These authors used this relationship to estimate salinity values during the Last Glacial Maximum (LGM).Interestingly, they observed several overestimations with regard to other published values in samples characterized by high relative abundances of larger specimens of E. huxleyi (> 4 µm).These discrepancies suggest that E. huxleyi (> 4 µm) in ancient sediments lacks an analog in modern assemblages, as indicated by the high dissimilarity between fossil samples with high percentages of this species and modern samples (Fig. 5b).
Because the MAT is strongly dependent upon on the analogs selected (Telford and Birks, 2009) and since the WA-PLS2 reconstruction for H2 and H1 is more coherent with a freshwater inflow scenario, it seems that WA-PLS2 affords more reliable values than the MAT.Consequently, WA-PLS2-estimated SSS was chosen for our final interpretations.
Transfer functions assume that the ecological response of organisms to either the environmental variable of interest or to the linear combination of this important variable with others has not changed significantly over the time span represented by the fossil assemblage (Birks, 1995).The good agreement observed between PC1 fossil and the reconstructed SSS patterns from 16 ka onwards (Fig. 5c) suggests that the SSS transfer function fulfills this assumption back to 16 ka.Larger differences are observed from 25 to 16 ka, possibly promoted by the lack of analogs during this time span, discussed above.Consequently, the SSS reconstruction from 25 to 16 ka will not be discussed further.A decrease in salinity of about 0.6 ± 0.15 psu occurred from 15.4 to 14.6 ka (Fig. 6a).The global sea-level rise of ∼ 20 m during meltwater pulse 1a (mwp-1a) has been dated between 14.6 and 14 ka (Stanford et al., 2006, and references therein).Since this section covers 3000 years with no control point (Fig. 6a), it could be an artifact of poorly constrained chronology for this time interval.Nevertheless, this seems unlikely because other authors (Duplessy et al., 1992;Emeis et al., 2000;Kallel et al., 1997) have reported SSS decreases in different regions of the Mediterranean Sea and Atlantic Ocean at this time from a combination of oxygen isotope (δ 18 O) and SST data.These salinity decreases are larger than that observed for the CEUTA10PC08 core.For instance, Duplessy et al. (1992) identified a salinity drop of about 2.5 psu in an Atlantic core west of the Strait of Gibraltar.It is worth mentioning that the salinity changes estimated by this method depend strongly on the accuracy of the SST record (Schmidt, 1998) and the unknown salinity-seawater δ 18 O relationship in the past (Rohling, 1999), being sensitive to several deviations and uncertainties that are difficult to assess (Rohling, 2000;Rohling and Bigg, 1998;Schmidt, 1999).Although the uncertainty in the chronology prevents a robust correlation, the smaller SSS decrease identified in the SSS reconstruction could be related to the Laurentide ice sheet melting and retreating at ∼ 15.5 ka (Clark et al., 2001).This event has already been proposed to be the cause of the freshwater input identified at 15.3 ka south of Iceland via advection within the North Atlantic Current (NAC) and subsequently its northern branch (Thornalley et al., 2010).Similarly, the southeastern branch of the NAC could have advected freshwater to the study area.

Bølling-Allerød (B-A)
The SSS values are generally low for the B-A, the Bølling being fresher than the Allerød (Fig. 6a).Owing to the global sea-level rise during the B-A, a greater volume of AW would have entered through the Strait, decreasing the average SSS.This period of reduced salinity also coincides with the highest values of total concentration of C 37 alkenones, a proxy of organic-matter preservation, from a nearby core located off the coast of Malaga (Ausín et al., 2015) (Fig. 6b).This accumulation of high amounts of organic matter resulted in the formation of the so-called organic-rich layer (ORL-1) (Cacho et al., 2002) in the western Mediterranean, although its origin is still under debate (Rogerson et al., 2008;Rohling et al., 2015).The joint effect of a salinity reduction of 0.8 psu and a temperature increase of 3 • C (Cacho et al., 2001) (Fig. 6c) would have led to a significant reduction in sea surface density, possibly prompting stagnation of the upper water column.This, along with increased organic-matter export to the seabed (Ausín et al., 2015) and reduced deep-basin ventilation (Martínez-Ruiz et al., 2015), would have hampered organic-matter mineralization, reinforcing the formation of the ORL-1 in the Alboran Sea.According to Rohling et al. (2015), the origin of ORL-1 lies in hydraulic changes in the Strait of Gibraltar (Bernoulli aspiration depth) and/or the inhibition of deep-water formation in the Gulf of Lion, both resulting from a drastic reduction in seawater density.Those authors have shown that the mwp-1a and the monsoon flooding into the eastern Mediterranean were insufficient to trigger these mechanisms and demonstrated that the Alpine meltwater input into the NW Mediterranean at this time (Ivy-Ochs

Younger Dryas (YD) and the Holocene
The YD exhibits a shift from higher to lower SSS values, decreasing by a total of 0.6 psu along its two phases: YDa and YDb (Fig. 6a).Several large short-term SSS fluctuations occurred from the onset of the YD throughout the early Holocene (up to 8 ka).This time span coincides with a sea-level rise of ∼ 30 m (Peltier and Fairbanks, 2006) due to short-lived freshwater inputs associated with residual melting of the Northern Hemisphere ice sheets (Andrews and Dunhill, 2004;Elmore et al., 2015;Seidenkrantz et al., 2013;Tornqvist and Hijma, 2012).Six brief periods of an SSS decreasing trend were identified at 12. 77-12.06, 11.95-11.71, 11.24-11.00, 10.09-9.83, 9.30-9.12, and 8.95-7.90ka (Fig. 6a).REDFIT spectral analysis reveals a periodicity of 770 ± 40 years (Fig. 6d), very similar to the 730 ± 40 years cycle found by Cacho et al. (2001) in an SST record in the Alboran Sea, which was punctuated by the so-called Alboran cooling events (AC_events) (Fig. 6c).Although this similarity does not necessarily imply a causal relationship, the timing of SSS decreases is comparable to that of the AC_events (Table 3), suggesting a common origin.Cacho et al. (2001) have associated the AC_events with influxes of cold Atlantic waters in the Alboran Sea during icerafted debris discharges (so-called Bond events) (Bond et al., 1997) (Table 3).These latter authors noted that the oxygen isotopic record showed no evidence of any of the coolings found for each Bond event during the Holocene and argued that the cooler surface waters may also have been fresher, offsetting the expected temperature-driven δ 18 O enrichment in their records.Similarly, the highly resolved δ 18 O profile reported by Cacho et al. (2001) does not show any of the expected oxygen isotopic enrichments associated with the AC_events, supporting the presence of fresher waters at those times.We suggest that freshwater advection (FA) events (as well as AC_events) would have resulted from the influx of fresher and colder Atlantic waters in the Alboran Sea related to the southeastward drifting of meltwater from the Labrador, Greenland and Iceland seas (Bond et al., 1997).FA events only occurred during the early Holocene, while AC_events and Bond events have also been identified through the middle and late Holocene.Wanner et al. (2014) concluded that, unlike the events occurring later, early Holocene Bond events originated from changes in the meridional overturning circulation due to meltwater pulses from the Northern Hemisphere ice sheets.It is likely that FA events would only have been noticeable when this mechanism was operating (i.e. the early Holocene), since very little meltwater was present after that period (Elmore et al., 2015).
An SSS increase of 0.87 ± 0.15 psu is observed from 10.7 to 10 ka.Because the western Mediterranean is a semienclosed basin, local conditions may have played a role as additional feedbacks for this rapid high-amplitude variability.For this brief period, Frigola et al. (2008) have demonstrated the most pronounced weakening of the Mediterranean thermohaline circulation for the last 50 ka.The consequent reduction in Atlantic-Mediterranean water exchange, along with the maximum summer insolation and inland aridity (Fletcher et al., 2010), would have led to more saline surface waters.
FA1 includes the 8.2 ka event (Alley et al., 1997), which has been linked to a sub-thermocline freshening of 0.5 psu in the North Atlantic (Thornalley et al., 2009).However, no distinctive SSS changes are observed in relation to this event, suggesting that it would have had a negligible effect on surface salinity in the Alboran Sea.Minimum SSS values are recorded at 7.8 ka, possibly related to maximum high-stand conditions reached at 7.4 ka (Zazo et al., 2008), along with the influence of the African Humid Period (AHP; 11-5.5 ka) over the study area, especially up to its decline at 7.4 ka (de-Menocal et al., 2000).From 7.8 to 4.5 ka, salinity values level off around 36.6 psu, close to present SSS values.

Conclusions
Multivariate statistical analyses show that the distribution of modern coccolithophore assemblages in the Atlantic Ocean, west of the Strait of Gibraltar, and the western Mediterranean was mainly influenced by annual mean salinity at 10 m depth.MAT and WA-PLS2 calibration models show similar outcomes.These models were applied to coccolithophore assemblages from a fossil core to reconstruct SSS at a high resolution for the last 25 kyr in the Alboran Sea.Statistical analyses reveal assemblages lacking good modern analogs in relation to the species E. huxleyi > 4 µm during H2 and H1 and part of the LGM, preventing further interpretations for these periods.A low SSS was found for the B-A, possibly due to the post-glacial sea-level rise.The consequent reduction in sea surface density is suggested to have reinforced the formation of the ORL-1.During the YD and Holocene, six brief, abrupt SSS decreases at 12. 77-12.06, 11.95-11.71, 11.24-11.00, 10.09-9.83, 9.30-9.12, and 8.95-7.90ka were linked to the advection of fresher and colder AW related to the southeastward drifting of meltwater in the North Atlantic.No evidence of the 8.2 ka event is found in the reconstructed SSS, which reached its lowest values at 7.8 ka, coinciding with high-stand conditions in the Alboran Sea and the onset of the decline of the African Humid Period.SSS remained low from 7.8 to 4.5 ka, close to its present values.
A broader understanding of the ecological link between coccolithophore species and environmental parameters would be desirable in order to be able to place coccolithophore-based transfer functions within an ecological context in future works.Nevertheless, the diverse statistical tests performed in this study and the strong emphasis placed on assessing the validity and reliability of both the model and the reconstruction do reveal the potential of coccolithophores for developing transfer functions.The derived transfer function provides a potential independent proxy for quantitative reconstructions of SSS changes in other locations of the western Mediterranean Sea over the last 15.5 kyr.
The Supplement related to this article is available online at doi:10.5194/cp-11-1635-2015-supplement.

Figure 1 .
Figure 1.Maps of the study area.Panel (a): location of the CEUTA10PC08 core (red star).Black arrows trace general surface circulation.Legend: AW -Atlantic water; MAW -Modified Atlantic Water; AC -Algerian Current; NC -Northern Current.Panel (b): location of the 88 core-top samples used for final calibrations.Maps generated with Ocean Data View software(Schlitzer, 2014).

B.
Ausín et al.: Development of coccolithophore-based transfer functions

Table 1 .
Multivariate analyses results.λ 1 /λ 2 : individual CCA.Preliminary model coefficients from the MAT and WA-PLS2.R 2 boot : bootstrapped coefficient of determination between the observed and predicted values.RMSEP: root-mean-square error of prediction.* Variables determined by ordination based on AIC.