http://orcid.org/0000-0001-6596-5303
Figures
Abstract
We present paleo-water depth reconstructions for the Pefka E section deposited on the island of Rhodes (Greece) during the early Pleistocene. For these reconstructions, a transfer function (TF) using modern benthic foraminifera surface samples from the Adriatic and Western Mediterranean Seas has been developed. The TF model gives an overall predictive accuracy of ~50 m over a water depth range of ~1200 m. Two separate TF models for shallower and deeper water depth ranges indicate a good predictive accuracy of 9 m for shallower water depths (0–200 m) but far less accuracy of 130 m for deeper water depths (200–1200 m) due to uneven sampling along the water depth gradient. To test the robustness of the TF, we randomly selected modern samples to develop random TFs, showing that the model is robust for water depths between 20 and 850 m while greater water depths are underestimated. We applied the TF to the Pefka E fossil data set. The goodness-of-fit statistics showed that most fossil samples have a poor to extremely poor fit to water depth. We interpret this as a consequence of a lack of modern analogues for the fossil samples and removed all samples with extremely poor fit. To test the robustness and significance of the reconstructions, we compared them to reconstructions from an alternative TF model based on the modern analogue technique and applied the randomization TF test. We found our estimates to be robust and significant at the 95% confidence level, but we also observed that our estimates are strongly overprinted by orbital, precession-driven changes in paleo-productivity and corrected our estimates by filtering out the precession-related component. We compared our corrected record to reconstructions based on a modified plankton/benthos (P/B) ratio, excluding infaunal species, and to stable oxygen isotope data from the same section, as well as to paleo-water depth estimates for the Lindos Bay Formation of other sediment sections of Rhodes. These comparisons indicate that our orbital-corrected reconstructions are reasonable and reflect major tectonic movements of Rhodes during the early Pleistocene.
Citation: Milker Y, Weinkauf MFG, Titschack J, Freiwald A, Krüger S, Jorissen FJ, et al. (2017) Testing the applicability of a benthic foraminiferal-based transfer function for the reconstruction of paleowater depth changes in Rhodes (Greece) during the early Pleistocene. PLoS ONE 12(11): e0188447. https://doi.org/10.1371/journal.pone.0188447
Editor: Fabrizio Frontalini, Universita degli Studi di Urbino Carlo Bo, ITALY
Received: July 7, 2017; Accepted: November 7, 2017; Published: November 22, 2017
Copyright: © 2017 Milker et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All data needed to replicate our analyses are stored and freely available from DRYAD under doi: 10.5061/dryad.349h0.
Funding: The sampling campaigns and sample preparation benefitted from funds of the DFG project Fr 1134/7 to Andre Freiwald. The German Research Foundation (DFG) supported for financial support through the Cluster of Excellence CliSAP (Integrated Climate System Analysis and Prediction) at the University of Hamburg.
Competing interests: The authors have declared that no competing interests exist.
1 Introduction
Since first quantifications of past sea surface temperatures (SST) by Imbrie and Kipp [1], the number of studies that applied transfer functions to microfossil assemblages have increased substantially. Transfer functions (TF) are now widely used to reconstruct paleo-environmental conditions, as for instance past relative sea levels, paleo-SSTs, paleo-precipitation and -salinity [2–10]. The majority of studies which use foraminiferal-based TFs to reconstruct relative paleo-sea levels are restricted to intertidal and shelf areas, and are biased towards younger time intervals such as the Holocene or latest Pleistocene periods. One exception in marine settings is the study of Avnaim-Katav et al. [11], who used a TF to reconstruct relative sea-level changes in the shallow Eastern Mediterranean Sea over a longer time interval, i.e., over the last 1 Myr. However, this study is limited to reconstructions for a few interglacial periods due to the shallow location of the studied sediment cores, lacking suitable glacial sediments. A number of studies estimated changes in paleo-water depths on longer time-scales, e.g., during the Miocene [12, 13] or in the Neogene [14, 15], Neogene to Holocene [16, 17] and early Pleistocene [18]. These studies used different approaches. Some applied semi-quantitative estimates based on water-depth ranges of modern faunas [16, 18] or weighted arithmetic and geometric means of the depth ranges of selected benthic foraminifera to be applied to fossil sequences [12, 13, 17]. Others used the modern analogue technique (MAT; [19]) [14] or MAT in combination with paleo-water depth estimates based on modified plankton/benthos ratios [15].
A prerequisite for accurate paleo-environmental reconstructions is that microfossils have a close relation to the environmental variable of interest [20]. Various studies have shown that benthic foraminifera, particularly in intertidal environments, can be used as accurate relative sea-level indicators [21–28] because they show a strong vertical zonation with respect to elevation in natural intertidal environments [29–31]. This is best explained by the sensitivity of intertidal foraminiferal assemblages to the duration and frequency of tidal flooding and related abiotic factors such as salinity [24, 29, 32]. Other studies indicate that shallow-marine foraminifera can also be used to reconstruct former relative sea levels with different quantitative methods [11, 33–36]. However, in these shallow marine settings (i.e., in the Mediterranean Sea), the foraminiferal distribution is not directly influenced by water depth but by abiotic and biotic factors changing with water depth, in particular by food flux to the sea floor and by substrate type, related to bottom water currents [37–41]. Benthic foraminifera TFs have in contrast never before been used to reconstruct water depths of deeper, i.e., outer shelf to upper slope habitats. This is unfortunate, because here they bear the potential to aid considerably in paleo-water depth reconstructions as water depth (and related gradients in food supply) is likely the major influential factor on the assemblage composition in those deeper environments.
For accurate paleo-environmental reconstructions it is, however, of further importance that the applied TF approach can accurately reproduce the underlying species-environment relationship in a data set. For instance, one often applied method for relative sea level reconstructions is the weighted averaging-partial least squares (WA-PLS) method [42] that is based on unimodal species distribution along the targeted environmental variable. Indirect and direct constrained ordination techniques such as detrended correspondence analysis (DCA, [43]) or detrended canonical correspondence analysis (DCCA, [44]) are widely used to determine whether species have a linear or unimodal response in a data set or along a specific target variable. However, these techniques assume a symmetrical unimodal response shape which is often not the case [45–47]. Alternative approaches can analyze the data directly, by comparing the shape of species distributions with respect to the targeted environmental factor, and determine their response by fitting a reaction curve [48, 49]. Species response curves (coenoclines) towards an environmental factor can be estimated by a variety of methods, including generalized additive models (GAM) and LOESS smoothing [50–53]. These methods have the advantage that the decision is based on real observations of the species response, instead of a derived value as in ordination analyses. Other issues that need to be addressed when using TFs to reconstruct past environmental conditions are unevenly sampled environmental gradients [54], spatial autocorrelation in the data set used to develop the TF [55, 56], the lack of modern analogues [57, 58], and the significance of the reconstructions [59]. Only a thoroughly tested TF can in the end be used to reliably reconstruct the target parameter, by making sure that no other parameters are leaking through or even dominating the derived reconstruction.
The major goal of our study is thus to test the applicability of a regional TF based on modern shelf-to-slope foraminiferal assemblages from the Western Mediterranean and the Adriatic Seas to quantify changes in paleo-water depths in an early Pleistocene sediment record from the island of Rhodes (Eastern Mediterranean) (Fig 1). We investigate the Recent species-environment relations and test the accuracy and robustness of our TF that is then applied to the fossil record. We further determine the reliability of the paleo-water depth estimates, and evaluate its usefulness for future quantitative reconstructions.
A. The Mediterranean Sea showing the location of the modern samples (Adriatic Sea [37, 60], Western Mediterranean Sea [61]) used for the development of the transfer function. B. The island of Rhodes with location of the Pefka E section and main Plio- to Pleistocene siliciclastic and carbonate dominated outcrops after [62, 63]. Maps created with Ocean Data View [64].
2 Study area
The basement of Rhodes consists of Mesozoic to Cenozoic (Jurassic to Paleocene) rocks, mainly carbonates, which are overlain by terrestrial Mio–Pliocene clastic sediments [62]. The Pliocene clastic fluvial and lacustrine sediments and marine Pleistocene sediments were deposited in isolated basins, especially along the eastern coast of Rhodes [65] (Fig 1B). According to the most recent lithostratigraphic scheme [66], the Plio–Pleistocene deposits are subdivided into three synthems (Trianda, Rhodes and Lindos–Acropolis Synthems), reflecting tectonically-driven transgression–regression cycles over multiple timescales, that are bound by major unconformities. The Lindos Bay Formation (LBF) mainly studied here is part of the Rhodes Synthem that can be subdivided into the transgressive Kolymbia Lindos Bay and St. Paul’s Bay Formations [67], and the regressive Cape Arkhangelos and Tsampika–Ladiko Formations. According to Hanken et al. [63], the up to 20 m thick Kolymbia Formation consists of fossil-rich grain- and packstones with increasing content of terrigenous silt and clay towards the top. The up to 30 m thick Lindos Bay Formation [62, 63] consists of predominantly homogenous (bioturbated) blue–gray marls with a rich microfossil content [18, 63, 68]. Moissette and Spjeldnæs [68] reported an initial increase in water depth followed by a shallowing upwards towards the top of the LBF based on bryozoan zonations. The Cape Arkhangelos Formation consists of a calcarenite that contains bryozoans, mollusks and brachiopods in the lower part, and bivalves and red algae (including large rhodoliths) in the middle and upper part, representing a shallow-marine carbonate environment with a shallowing upward trend [63, 69].
The 16.4 m long Pefka E section is situated on the SE coast of Rhodes (36°3’50”N; 28°3’58”E; Fig 1B), 1.5 km south-east and 3.5 km south of the villages Pefka and Lindos, respectively. The section covers the uppermost part of the Kolymbia Formation, the Lindos Bay Formation, and the lowermost part of the Cape Arkhangelos Formation deposited on an erosional surface on top of the LBF (Fig 2). The Kolymbia Formation, at the base of the section (from that the upper 40 cm were sampled), consists of sandy to clayey brown bioclastic limestones. It gradually fines upward into the LBF within 27 cm. The LBF in the Pefka E section has a total thickness of 15.4 m and mainly consists of gray-brown and olive-green homogenous and bioturbated marls with regular intercalations of 5 to 20 cm thick brownish to slightly reddish-gray laminated marl intervals. Within the section, these laminated marl intervals are mainly restricted to the lower and middle parts of the LBF. The upper part of the LBF contains a thin horizon with iron-impregnated cavities and the uppermost 50 cm shows intense bioturbation. The burrows of Thalassinoides penetrate the LBF from the erosional unconformity on top of the LBF and are filled with calcarenites of the overlying Cape Arkhangelos Formation.
3 Material and methods
3.1 Chronostratigraphic framework and stable isotope analysis of the Pefka E section
The chronostratigraphic framework of the Pefka E section is based on a combination of biostratigraphic information and cyclostratigraphy. Biostratigraphic information came from the distribution of calcareous nannofossils. For example, Frydas [70] suggested a deposition within nannoplankton zone NN 19 (see also Titschack et al. [66]), and recently Quillévéré et al. [71] provided a detailed chronostratigraphy of the section, suggesting a deposition from NN 19b to NN 19d based on the highest occurrence (HO) of Calcidiscus macintyrei, lowest occurrence (LO) and HO of large Gephyocapsa, and HO of Helicosphaera sellii. Nannoplankton samples analyzed at 106–109 cm and 1106–1108 cm (from the top) suggest a deposition during nannoplankton zone NN19b based on the occurrence of Gephyrocapsa aperta, Rhabdosphaera claviner, and Pseudoemiliana spp., and the absence of Discoaster brouweri and C. macintyrei (Erlend Martini, pers. comm., 2005). For the cyclostratigraphic framework, we have used the first common occurrence (FCO) of Hyalinea balthica (in the >125 μm fraction), which was dated at 1.492 Ma in the Vrica section in Calabria, Italy [72]. The cyclostratigraphic approach is further based on correlation of the filtered time series of the relative abundance of epifaunal species in the >125 μm fraction and the Northern Hemisphere (65°N) summer insolation (June 21) [73]. The relative abundance of epifaunal species was used because they respond sensitively to climate-driven (i.e., precession-driven) changes in oxygen concentrations and food supply at the sea-floor due to sapropel (black shale) formation [74]. Blackman–Tukey spectral analyses of the relative abundance of epifaunal species revealed maxima at 167 and 104 cm, likely representing obliquity and precession cycles and reflecting climate-related changes in trophic conditions. Subsequently, a Gaussian band pass filter, centered at the period of 104 cm, was calculated for comparison with the filtered precessional component (23 kyr cycle) of Northern Hemisphere insolation. Spectral analyses, band pass filtering, and correlations were performed with the software AnalySeries v. 2.0 [75].
Stable oxygen isotope analyses were measured on the epibenthic foraminifer Cibicidoides pseudoungerianus with the exception of one sample, where Cibicidoides pachyderma s.l. was used due to the lack of C. pseudoungerianus. Only pristine and clean specimens were picked from the 250–355 μm fraction. We disregarded samples where the calcareous species showed distinct signs of dissolution (i.e., two lowermost samples of the section). The isotopic composition of two to three specimens per sample was measured at the University of Leipzig. The samples were reacted with 105% phosphoric acid at 70°C in a Kiel IV online carbonate preparation device connected to a Finnigan MAT 253 mass spectrometer. Results are given in δ-notation, normalized to the Vienna Pee Dee Belemnite scale (VPDB). External precision (1σ), as determined by replicate analysis of the NBS19 carbonate standard (1 standard per 5 samples), was on average ±0.055‰ for δ18O.
3.2 Fossil foraminiferal investigations in the Pefka E section
In the years 2000 (lower part of the section) and 2001 (middle and upper parts of the section), 162 samples of 2 cm thickness each were obtained mainly every 10 cm from the 16.4 m long Pefka E section (Fig 1B). Prior to sampling, the weathered surface was thoroughly removed in order to get access to fresh material. All samples were treated with 10% hydrogen peroxide for 24 hours to disaggregate the sediment and were subsequently wet-sieved over a 63 μm screen, dried at 40°C, and then dry-sieved over a 150 μm screen. Foraminiferal investigations were carried out on representative splits containing at least 200–300 benthic individuals of the >150 μm sediment fraction. The >150 μm fraction was used because foraminifera in the modern data used in this study were counted from the same size fraction. Four samples from the Pefka E section (at 617, 1057, 1317, and 1427 cm depths) contained less than 100 specimens per sample, and were thus excluded from further analyses. We further excluded the two lowermost samples (at 1630 and 1640 cm depths) because we noted strong calcite dissolution in these samples, and the six uppermost samples (between 0 and 43 cm depth) because this part of the section is intensely bioturbated and the majority of foraminiferal tests are badly preserved. For the taxonomic identification, preferentially on species level, we used standard illustrations and descriptions [76–80]. When identification on species level was not possible, individual species were grouped into their genus. Potentially reworked specimens (i.e., broken tests, test fragments) were not counted.
3.3 Compilation and testing of a combined modern data set for paleo-water depth estimates
For the paleo-water depth estimates, we used two sets of modern surface samples from the Adriatic Sea (269 samples ranging from 8 to 1216 m water depth [37, 60]) and the Western Mediterranean Sea shelf (45 samples ranging from 20 to 235 m water depth, [61]) (Fig 1A). The samples from the Western Mediterranean Sea were re-counted in the >150 μm fraction to ensure compatibility with the Adriatic Sea data set. We compiled these data to cover a wide range of environmental variability in the modern Mediterranean Sea and to ensure that the main fossil species are present in the modern data set. The integration of samples from distant locations further may help to limit the issue of spatial autocorrelation [55, 56]. In the Adriatic Sea dataset, we removed 57 samples with >500 specimens/g sediment, because these samples correspond to lag deposits, as is shown by the large amounts of abraded, broken, and re-crystalized tests. Such samples, in which foraminifera have been strongly concentrated, often contain a mixture of Recent and glacial faunas, so that they are no longer representative of modern conditions and water depth. Next, we excluded all modern samples (n = 55) taken from water depths <20 m from the Adriatic Sea data set, because from the fossil faunal composition in the Pefka E section we did not expect to reconstruct paleo-water depths shallower than this. This further helps to reduce the effect of uneven sampling along the water depth gradient, because it eliminates a relatively oversampled depth range [54]. For consistent species assignment of the modern and fossil species, we investigated the individual taxonomic concepts by comparing SEM illustrations and descriptions [37, 60, 80]. The combined modern data set was finally reduced to species with a relative abundance of ≥ 3% in any one sample to exclude very rare species, at the same time ensuring that both data sets contain the majority of abundant species. The modern final data set contains 199 samples and 70 individual taxa (S1 File). In the fossil data set, we excluded all species that are not present in the combined modern data set. This reduced fossil data set contains on average 79% of individuals in the full data set. All datasets (fossil and modern) have been reclosed (i.e., standardized to 100% per sample) using the species that were finally retained in our analyses. In this way, the subcompositional coherence of the datasets is maintained, as is recommended practice for invariable analyses in subcompositional datasets [81].
We applied DCA to the combined modern data to test whether the species have a linear or unimodal response using R (version 3.3.0; [82]) and the R-package “vegan”, version 2.3–3. [83]. We further applied the reaction norm-test written in R (S2 File) to determine the individual species’ response to the water-depth gradient by using four different methods. In all cases, the relationship between the abundance of each species with the environmental factor (here: water depth) is analyzed for its shape (i.e., linear vs. non-linear behavior).
The first method (goodness-of-fit) fits a linear and a quadratic function to the data and tests whether the quadratic function describes the data significantly better, using the F-distribution [84], which would imply a non-linear relationship between species abundance and environmental factor. The second method uses piecewise robust linear regression along a moving window to determine the species responses to water depth. Here, the data set is divided into bins of equal size along the range of the environmental parameter, and a robust linear regression based on the M-estimator [85] is applied to every segment. When the slope of the regression line does not significantly differ between all segments, the relationship between abundance and environmental parameter is assumed to be linear along the entire gradient. The third method (binning) calculates mean values of abundance in several bins of equal size along the environmental gradient. When the abundance of the species across all bins either remains constant or else decreases/increases by a constant value from bin to bin (within the range of the 95% confidence interval), a linear relationship can be assumed for this species. The fourth method uses an estimate of the relationship between species abundance and environmental parameter based on a generalized additive model (GAM) on the Gaussian distribution with identity as link function [50, 52]. The Akaike information criterion [86] is used to decide, whether the response curve is linear or non-linear. The four reaction norm test methods have been evaluated using simulated coenoclines which represent truncated or full unimodal as well as linear responses. The tests have shown that the goodness-of-fit and GAM methods are most reliable, with 96.5%, 88.5% and 100%, and 95.5%, 87.0% and 100% correct classifications, respectively, for linear, truncated unimodal and full unimodal behavior. Binning, independent of the number of bins selected, is conservative against a linear response (with only 9.5% correct classifications) and piecewise robust regression is conservative against a truncated non-linear response (with only 21.5% correct classifications) (Table A in S3 File).
As a result of our comparisons between the DCA and reaction norm test results (see section 4.2.1), we applied canonical correspondence analysis (CCA;[44]), with water depth as the only environmental variable, to estimate the amount of variance in the combined modern data set that is explained by the water depth-gradient. We ran additional CCAs on the modern data sets from the Adriatic Sea and the Western Mediterranean Sea to explore species–environment relations in the individual data sets. We again used CCAs for both data sets because DCAs indicate a unimodal species response in both modern data sets (see section 4.2.1). For the Adriatic Sea data set, all species were included into the analysis and water depth, sand content (%), and organic carbon content (%) were used as environmental parameters. For the Western Mediterranean Sea data set, we only used species with a relative abundance of ≥5% in at least one sample. Here, we used water depth, chlorophyll a content (mg/m3) of the surface water, and substrate (fine-grained sediment <63 μm in per cent) as environmental parameters. Two samples from the Western Mediterranean Sea (samples 352–1 and 393–1) were excluded from the analysis because no grain size data were available for them. All environmental parameters were standardized, and we added binary co-variables to the CCA analysis of the data set from the Western Mediterranean Sea in order to minimize the differences between the three investigated areas [61] (Fig 1). All CCA analyses were performed with the R-package “vegan”.
3.4 Development of a regional transfer function for paleo-water depth estimates
For the development of a regional TF, we used WA-PLS. The TF performance was evaluated on the basis of the bootstrapped (1000 cycles) coefficient of determination (R2boot), allowing for a first estimation of the strength of the linear relationship between the observed and estimated water depths in the training data, and the root mean squared error of prediction (RMSEP), allowing for the evaluation of the overall predictive ability of the TF. To avoid model overfitting, a maximum of three components were extracted. The component used eventually was selected according to the lowest RMSEP value if the reduction in prediction error is significant and exceeded 5% for this component compared to the next lower component [42, 87]. Bootstrapping cross-validation (1000 cycles) was also used to evaluate the sample-specific errors of prediction in the fossil data set. We further applied the weighted-averaging (WA) approach to calculate optima and tolerances of the modern species with respect to water depth. For WA-PLS and WA calculations, we used the R-package “rioja”, version 0.9–9 [88].
To test the robustness of our regional TF, we removed randomly selected samples from the modern Adriatic Sea and Western Mediterranean Sea data sets (25% of the samples from each basin) and developed random TFs (n = 100) in “rioja” from the remaining modern samples which were then applied to the randomly removed samples.
To determine the individual species’ influence on the reconstructions in the Pefka E section, we applied jackknifing. For this, we reconstructed TFs leaving out one species at a time, and calculated the Kendall rank-order correlation [89] between the final TF model reconstruction and the jackknifed models, as well as the RMSEP of each reconstruction.
3.5 Statistically testing the accuracy of paleo-water depth estimates
We applied the goodness-of-fit statistics to evaluate whether the fossil samples from the Pefka E section have a good fit to the water depth in the modern data [57, 58, 90]. We used the 90th and 95th percentiles in the goodness-of-fit statistics as thresholds to differentiate between samples with good (<90), poor (90–95), and very poor (>95) fit to water depth and applied CCA for the ordination. For the calculations, we used the R-package “analogue”, version 0.16–3 [91].
To determine the robustness of our reconstructions, we compared them with reconstructions using an alternative TF method (i.e., MAT). We selected the seven closest analogues in the combined modern data set following the suggestions in Kemp and Telford [90], and used bootstrapping (1000 cycles) to calculate RMSEP and the cross-validated coefficient of determination (R2boot) between observed and estimated water depth in the combined modern data set. For all calculations, we used the R-package “analogue”.
We further tested the significance of our reconstructions by applying the randomTF test [59]. This test calculates the explanatory significance of the TF by comparing its proportion of explained variance with that of many (n = 999) alternative models trained on random environmental data, using an ordination technique (e.g., redundancy analysis (RDA); [92]). For the random TF test, we used the R-packages “palaeoSig”, version 1.1–3 [93], “rioja”, and “vegan”. To choose the correct ordination method, we determined whether the species have a linear or unimodal distribution in the fossil data set by using DCA in R package “vegan”.
All data needed to replicate our analyses are available from DRYAD under doi: 10.5061/dryad.349h0.
4 Results and discussion
4.1 Age model
At the Pefka E section, the FCO of H. balthica in the >125 μm fraction was identified at 917 cm section depth representing an age of 1.492 Ma [72] (Table 1). The combined biostratigraphic and cyclostratigraphic framework revealed a depositional age ranging from 1.656 Ma near the base to 1.278 Ma at the top of the LBF suggestion a Calabrian (early Pleistocene) age of the section (Fig 3). Thus, the 16.4 cm thick Pefka E section covers ~387 kyrs. Sedimentation rates range from 2.80 to 5.94 cm kyr−1 (average of 4.35 cm kyr−1) with generally higher values in the lower part and upper parts and lower values in the middle part of the section (Fig 4). Although we cannot exclude potential errors in peak-to-peak correlations, our age model is in general accordance with available biostratigraphic information on the appearance of calcareous nannofossils (Gephyrocapsa aperta, Rhabdosphaera claviner and Pseudoemiliana spp.) placing the LBF succession of Pefka E into early Pleistocene nannofossil zones MNN19b–d (E. Martini, pers. comm, 2005; [70, 71]). However, the planktonic foraminiferal biostratigraphy revealed considerably older ages of 1.79 to 2.09 Ma for the sediments of the boundary interval between the here investigated LBF and the underlying Kolymbia limestone [71] (Fig 2). This discrepancy may be attributed to reworking of older sediment particles into (both) the Kolymbia limestone and basal LBF leading to overestimated ages for this part of the section.
Results of spectral analysis and Gaussian band pass filtering of the relative abundance of epifaunal species (in the >125 μm fraction) in the Pefka E section and the Northern Hemisphere summer insolation [73]. Maxima of the filtered 102 cm component of relative epifaunal abundance were correlated with the filtered 23 kyr component of Northern Hemisphere summer insolation. The correlation is referenced to the first common occurrence (FCO) of Hyalinea balthica, which is dated at 1.492 Ma [72].
Age-depth plot of the Pefka E section with calculated sedimentation rates. Age tie points are derived from the first common occurrence (FCO) of Hyalinea balthica and graphical correlation of filtered relative epifaunal abundance and Northern Hemisphere summer insolation. See also Table 1.
4.2 Development of the regional transfer function
4.2.1 Response of the species in the modern data set to water depth.
The accuracy of a transfer function for paleo-water depth estimates strongly depends on whether or not modern species have a close relation to water depth. To investigate the response of foraminifera in the combined modern data set along the water-depth gradient, we first calculated the length of the gradient along the first axis by using DCA. We measured a gradient length of 4.6 standard deviation (SD) units, which indicates that most species should have a unimodal distribution (Table B in S3 File). However, to determine the individual species response in more detail and along the water depth gradient, we applied our reaction norm-test. The goodness of fit statistics indicates that 57.1% of the modern species have a significant (p<0.05) non-linear response to water depth (Table C in S3 File). By using the GAM, 71.4% of the species have a non-linear rather than a linear response to water depth. Piecewise robust regression, using three bins, suggests that only 41.4% of the species have a non-linear response. Binning (three mean bins) indicates 85.7% of the species as having a non-linear response to the water depth gradient, including all of the most dominant species.
To determine the amount of variance in the combined modern data set is explained by water depth, we applied CCA, with water depth as the sole environmental variable. The results show that 11.2% of the variance in the combined modern data set is explained by the constrained (first) axis (i.e., water depth), and 16.5% by the second, unconstrained, axis (Table B in S3 File). The water depth gradient is significant with p<0.001 (1000 permutations), but the ratio of the eigenvalue of the constraint axis (λ1) to the unconstrained axis (λ2) is 0.68, indicating that water depth is not the sole ecological gradient in the modern data set [20]. Among the most dominant fossil species in the Pefka E section, and hence important for the water-depth estimates (see section 4.3.1), Gyroidinoides altiformis, Uvigerina peregrina, Planulina ariminensis, Cibicidoides pachyderma s.l., Sphaeroidina bulloides, and to some extent Hyalinea balthica show a close relation to the first axis and thus the water depth gradient. In contrast, other dominant fossil species (Cibicides pseudoungerianus, Globocassidulina subglobosa, Cassidulina obtusa, Cassidulina carinata s.l., Lobatula lobatula s.l., Melonis barleeanum and Bulimina marginata s.l.) have a close relation to the second axis (Fig 5). Therefore, the CCA indicates that at least one additional important abiotic parameter influences the modern species distribution. This interpretation is in accordance with earlier observations, where substrate type and food availability (as reflected by total sedimentary organic carbon contents (OC) in the Adriatic Sea and by chlorophyll a values reflecting surface production in the Western Mediterranean Sea) have been identified as further important environmental factors [37, 61]. A CCA of the modern data set from the Adriatic Sea shows water depth, ranging between 20 and 1216 m, to be closely related to the first CCA axis explaining 17.5% of the variance in the data set (Table B in S3 File). Sand content and OC are related to the second CCA axis, which explains only 7.4% of the variance in the data set (Fig D and Table B in S3 File). The CCA of the modern data set from the Western Mediterranean Sea shows that water depth, ranging between 20 and 235 m, but also substrate type are the most important parameters for the species distribution. Both are closely related to the first axis, explaining 24.2% of the variance in the data set (Fig E and Table B in S3 File). In the combined modern data set, with water depths ranging between 20 and 1216 m, most of the species which plot on the negative side of the first CA axis are epifaunal species and have their optimum at relatively shallow water depths (Fig 5). This habitat preference suggests the relevance of other environmental factors than water depth, including a close relation to substrate type and bottom water currents, and to a lesser extent to organic carbon fluxes.
In bold the species that are dominant in the fossil data set are marked while the eigenvalues of the constrained (CCA) and unconstrained (CA1) axes are given in parentheses.
4.2.2 Influence of different trophic conditions on the water-depth distribution of benthic foraminifera in the modern Mediterranean Sea.
The bathymetric zonations of benthic foraminifera in the modern Mediterranean Sea reveal a characteristic shift from west to east, which can be related to the trophic contrasts between the meso- to eutrophic western basin and the oligotrophic eastern basin [39, 40]. Accordingly, there is a west to east shallowing of the upper or lower depth ranges of most oligotrophic taxa and a decrease or absence of eutrophic (i.e., shallow to deep infaunal) taxa [40]. These observations indicate that organic carbon flux to the seafloor, rather than water depth, is the driving factor for the bathymetric zonation and hence, the western and eastern basins may not be directly comparable due to their different trophic conditions. This raises the question whether or not the environmental conditions in the early Pleistocene were comparable to the modern oligotrophic conditions in the eastern basin. We found shallow to deep infaunal species throughout the section, indicating enhanced organic carbon fluxes to the sea floor. Specifically, we observed between 7% and 80% (on average 34%) infaunal species of which between 4% to 77% (on average 30%) are considered as shallow infaunal (species of the genera Bulimina, Brizalina, Cassidulina, and Uvigerina) and between 0% and 19% (on average 3%) as intermediate infaunal (i.e., M. barleeanum) in the Pefka E section (Fig F in S3 File). Shallow infaunal species, dominated by C. carinata s.l. and U. peregrina, are found in high relative abundance throughout the section with a slight decline in the middle part. Intermediate infaunal species are present in low percentages throughout but show a slightly higher abundance in the middle and upper parts of the section.
Some of the dominant fossil species, i.e., U. peregrina, C. carinata s.l., and B. marginata s.s., are supposed to be adapted to fresh organic matter pulses [94–102]. These species have only rarely been observed in the modern Eastern Mediterranean Sea while deep infaunal species are nearly absent [40] except in areas directly influenced by riverine nutrient input and locally enhanced organic matter fluxes [103]. The fossil foraminiferal distributions, i.e., the high relative abundance of infaunal species, suggest at least mesotrophic conditions in the Eastern Mediterranean Sea during the early Pleistocene. This is also supported by Plio- and early Pleistocene foraminiferal assemblages from the Kallithea section in NE Rhodes (Fig 1B; [18]), and by the occurrence of the non-zooxanthellate cold-water corals Lophelia pertusa and Madrepora oculata in the Lardos section (Fig 1B, [66]). Thus, trophic conditions in the early Pleistocene environments near the island of Rhodes may be better comparable to those of the modern Western Mediterranean Sea, and particularly to those of the modern Adriatic Sea then to the modern environments around Rhodes as also concluded by Rasmussen and Thomsen [18]). More specifically, the modern Adriatic Sea and the Western Mediterranean Sea comprise on average 24% and 8% infaunal species, respectively. Although trophic conditions in the early Pleistocene Eastern Mediterranean Sea and the modern Adriatic Sea (where 77% of the modern data set come from) may to some extent be comparable, we cannot completely exclude bathymetric offsets in the foraminiferal zonations between the modern Western Mediterranean and Adriatic Seas, and between the modern and Pleistocene environments.
4.2.3 Transfer function performance and reliability.
We ran an initial TF model on the modern data to identify outliers, that is, samples with residuals larger than 2.5 times of the standard deviation of the residuals (5 samples; 2.5%). These samples were then removed from further analysis [104]. For our final TF model, we observed a significant (p<0.001) improvement of the RMSEP of 54% from the first component (based on WA with inverse deshrinking) to the second component, and a further significant (p<0.001) improvement of 10% from the second to the third component (Table 2). We consequently selected the third component, having a cross-validated coefficient of variation R2 of 0.95 between the observed and estimated water depths in the combined modern data set and an RMSEP of 49 m, which corresponds to 4% of the observed water-depth range (20–1216 m) in the modern data set. This value is comparable to transfer functions from other intertidal and subtidal settings [11, 34, 35, 87], indicating that our model performs well. However, it should be noted that the RMSEP may be over-optimistic because of the unevenly covered water-depth gradient (i.e., due to the small number of samples from deeper water depths) in the modern data set. To address this issue, we split our data set into two segments (20–200 m with 171 samples and 200–1216 m with 28 samples) and calculated segment-wise RMSEPs as suggested by Telford and Birks [54]. These calculations give RMSEPs of 9 m for the shallower segment but of 130 m for the deeper segment. This shows that the overall predictive accuracy is, in fact, overoptimistic and decreases with increasing water depth. This can be explained by a combination of the low number of samples from larger water depths, and the wider tolerance of deep-water species relative to water depth (i.e., G. altiformis, U. peregrina, Gyroidina orbicularis, P. ariminensis, Bulimina costata, Pseudoclavulina crustata, Uvigerina mediterranea, Pyrgoella irregularis and C. obtusa; Fig G in S3 File). The comparison of the observed and estimated water depths in the combined modern data set shows a good fit between both, with residuals ranging from −144 to 132 m water depth (Fig 6).
Shown are the observed versus estimated water depths (A) and their residuals (B) in the combined modern data set (different colors and symbols are given for the two modern data sets from the Adriatic and Western Mediterranean Seas). See also Table 2.
In TF models for intertidal environments, commonly the first (WA) component is selected because higher components provide no or only slight model improvements as a result of the dominant influence of marsh elevation (and co-correlated environmental variables) [105–107]. In contrast, TF models for shelf environments showed that higher components commonly yield model improvements due to the influence of more than one environmental variable on the species distribution [11, 34, 35]. The improvement in RMSEP from the first to the second and third components may be the result of the uneven sampled gradient and/or the influence of other environmental parameters. The first WA-PLS component is equivalent to the first axis of a CCA with a single environmental variable [104], and CCA results show that only 11.2% of the variance in the modern data set can be explained by the first axis but 16.5% by the first unconstrained axis (Table B in S3 File). This indicates that additional environmental variables affect the foraminiferal distribution in the combined modern data set. Consequently, the model improvement for the second and third components suggests that WA-PLS, which corrects for edge effects, is also able, to some extent, to exploit the structured pattern in the residuals for higher components after fitting a WA model [42, 104]. However, higher components should not be chosen solely based on the model improvement to avoid a model overfitting [87, 90]. In our case, the selection of the third component seems to be ecologically reasonable because of an obvious influence of additional environmental variables (such as substrate type and food availability) on the modern species distribution (Fig 5, Figs D and E in S3 File).
To test the reliability of our regional TF, we applied 100 randomizations, where the modern dataset was randomly separated into training- (75%) and test-data (25%). We note that there is a correlation between estimated and reconstructed water depths (τ = 0.783) based on a Kendall rank-order correlation. The correlations are better for water depths ranging between 20 m and 848 m, but for water depths >848 m all test TFs underestimate the measured water depths (Fig H in S3 File). The performance statistics is almost comparable to the primary TF (Table 2) with mean R2boot of 0.949 (SD of 0.02) and mean RMSEP of 51 m (SD of 4 m). The component used for the reconstruction was dynamically selected in each randomization according to the criteria described above, and the third component was used in 98 cases. These results imply that our regional TF seems to be robust and can reliably predict water depths as deep as ~850 m, but with a bias towards underestimation beneath that depth.
4.3 Application of the regional transfer function
4.3.1 Fossil foraminifera in the Pefka E section.
A total of 171 different taxa have been identified in the samples from the Pefka E section, with individual samples containing between 29 and 64 taxa. In the data set used for the application of the transfer function, a total of 39 and 28 taxa have an abundance of >1% and >3% in at least three samples, respectively. The most important fossil species (with a relative abundance of ≥10% in at least one sample) in this data set comprise C. pachyderma s.l. (1–51%), C. pseudoungerianus (0–48%), U. peregrina (0–46%), C. carinata s.l. (0–38%), B. marginata s.l. (0–33%), G. subglobosa (0–23%), G. altiformis (0–23%), M. barleeanum (0–20%), C. obtusa (0–19%), L. lobatula s.l. (0–19%), S. bulloides (0–18%), H. balthica (0–15%), and P. ariminensis (0–13%) (Fig 7). In general, we observed fluctuating relative abundances for the most dominant species in the Pefka E section, but in more detail, C. pachyderma s.l. and C. pseudoungerianus are dominant throughout the section with the former increasing in relative abundance towards the top and the latter being slightly more abundant in the lower and middle parts of the section. Uvigerina peregrina has a higher relative abundance in the upper part compared to the lower part of the section. Cassidulina carinata s.l. has a higher relative abundance in the lower part of the section but it also occurs frequently in the upper part of the section. Cassidulina obtusa exhibits two abundance maxima in the middle and upper parts of the section, respectively. Bulimina marginata s.l. has a generally low relative abundance except for two samples in the middle part of the section. Gyroidinoides altiformis is present in the lower part but is more abundant in the middle part and has not been observed in the upper part of the section. Globocassidulina subglobosa is present throughout the section but is more abundant in the lowermost and uppermost parts. Hyalinea balthica firstly occurs in the middle part of the section and is then permanently present until the top of the section. Melonis barleeanum, P. ariminensis, and S. bulloides are present throughout the section. Lobatula lobatula s.l. is present in the lower and middle parts of the section but exhibits a higher abundance in the uppermost part of the section.
Relative abundance of the most important fossil species (≥10% in at least one sample) in the Pefka E section versus age.
According to their optimum water depths in the combined modern data set, calculated by weighted averaging, the most important species G. altiformis (797 m), U. peregrina (745 m), and P. ariminensis (560 m) can be considered as indicators for deep water, C. obtusa (333 m), S. bulloides (325 m), C. pachyderma s.l. (288 m), and H. balthica (208 m) as indicators for intermediate water depths, and C. carinata s.l. (155 m), M. barleeanum (152 m), L. lobatula s.l. (106 m), B. marginata s.l. (104 m), G. subglobosa (98 m), and C. pseudoungerianus (88 m) as indicators for shallow water depths (Fig G in S3 File). Although we probably do not cover the complete water depth range of some of the most important species in the fossil record, our interpretations are supported by a number of observations from the Mediterranean Sea [16, 38–40, 76, 102, 108, 109].
4.3.2 Determining the robustness and accuracy of the paleo-water depth estimates in the Pefka E section.
The accuracy of paleo-environmental reconstructions depends not only on the species response to the target environmental variable but also on the availability of modern analogues for fossil faunas and the fit of the fossil species to water depth. To explore whether or not our fossil samples have a good fit to water depth, we applied the goodness-of-fit statistics [57, 58]. The goodness-of-fit test indicates that only 7 fossil samples (5%) have a good fit, 39 samples (26%) have a poor fit, but the majority of the fossil samples have a very poor (103 samples; 68.0%) or extremely poor (>99, 3 samples; 2%) fit to water depth (Fig 8). These results are intriguing, and are investigated in more detail in the following. A coverage plot of the maximum relative abundance of individual species in the modern and fossil data sets illustrates that a number of species have a higher relative abundance in the fossil compared to the modern data set (Fig I in S3 File). Among these, the most dominant (and thus most important) species in the Pefka E section have a higher maximum relative abundance (between 38% and 51%) but a lower maximum relative abundance (of 16% to 30%) in the modern data set. We therefore assume that the poor to very poor fit to water depth in the Pefka E fossil data set is the consequence of a lack of close analogues in the modern data, indicating that local paleo-environmental conditions (i.e., productivity) during the early Pleistocene perhaps were less comparable to the modern conditions than previously thought (see 4.2.2). A number of tests with simulated data showed, however, that WA-PLS performs well under mild analogue situations [110, 111]. Also, Hutson [112] tested how well different transfer function techniques work under non-analogue situations and found a technique based on weighted averaging to be most accurate, because it does not extrapolate as other techniques do.
Squared residual length of the fossil samples from the Pefka E section versus depth is shown. The 90th, 95th, and 99th percentile and the maximum residual length are indicated with lines.
To determine the individual influence of the species on the paleo-water depth estimates, we ran additional TFs by excluding species by jackknifing and compared the reconstructions of the derived TFs with our final model by calculating Kendall’s rank-order correlation and RMSEP. These results show that the exclusion of U. peregrina resulted in the lowest correlation (τ = 0.08) and highest RMSEP (278 m) in the Pefka E section, indicating that this species has the strongest influence on the estimates (Fig J in S3 File). Further species with a strong influence on the reconstructions (higher RMSEP) comprise C. pseudoungerianus (118 m; τ = 0.836), C. pachyderma s.l. (104 m; τ = 0.888), C. carinata s.l. (61 m; τ = 0.847), and G. altiformis (61 m; τ = 0.846) (Fig J in S3 File).
The species with a higher relative abundance in the fossil record compared to the modern data sets and/or with a strong influence on the estimations have their optimum modern water depths at 88 m (C. pseudoungerianus), 155 m (C. carinata s.l.), 288 m (C. pachyderma s.l.), 745 m (U. peregrina), and 797 m (G. altiformis) (Fig G in S3 File). Accordingly, the TF may overestimate paleo-water depths for samples strongly dominated by U. peregrina and G. altiformis and may underestimate paleo-water depths for samples strongly dominated by C. pseudoungerianus and C. carinata s.l. As a consequence, we screened the three samples with an extremely poor fit to water depth (squared residual lengths larger than the 99th percentile). Two of the samples are strongly dominated by both U. peregrina (45% and 23%) and C. pachyderma s.l. (35% and 50%), respectively. The third sample has a high relative abundance of S. bulloides (18%). We finally removed these three samples from the reconstructions.
To determine the robustness of TF reconstructions, it is recommended to compare different TF techniques [90, 104]. We therefore compared our paleo-water depth reconstructions based on the WA-PLS method to those using MAT. We found a cross-validated correlation between observed and predicted water depths in the combined modern data set of R2 = 0.98 and a RMSEP of 55 m for MAT, indicating a performance comparable to WA-PLS (compare Table 2). The reconstructions in the Pefka E section are very similar for both models (r = 0.77; Fig K in S3 File). This indicates that our reconstructions seem to be reasonable.
Recent studies [59, 113] highlighted the necessity to test the significance of paleo-environmental reconstructions when using microfossil-based TFs. We therefore applied the randomTF test [59] to test the significance of our reconstruction. We used RDA for ordination because of the relative short gradient length (2.15 SD units) along the first axis, indicating a more linear rather than unimodal species distribution in the fossil data set (Table B in S3 File). We ran the test to the data set excluding the three fossil samples with extremely poor fit to water depth as described above. The result shows that our paleo-water depth reconstructions are significant at the 95% confidence level (Fig 9), confirming that our reconstructions seem to be robust [59].
Histogram of the proportion of variance in the fossil data set explained by 999 transfer functions trained with random data, and the proportion of variance explained by the transfer function constructed by us (solid black line). Dotted red line shows the 95%-explained limit and dashed line the maximum variance in the fossil data set.
4.4 Paleo-water depth trends in the Pefka E section and comparison with other records
To determine the reliability of our paleo-water depth reconstructions, we here want to compare them to other records available from the island of Rhodes. In general, our paleo-water depth reconstructions show a long-term transgressive trend, followed by a regressive phase in the upper part of the section, which may reflect tectonic movements of the island of Rhodes during the early Pleistocene (Fig 10A). However, the reconstructions also show distinct and repeating short-scale fluctuations with amplitudes of up to around 800 m, particularly in the middle part of the section. These abrupt fluctuations cannot be explained by tectonic adjustments but may represent a bias linked to orbital-driven environmental changes and resulting bathymetric species shifts. The Neogene climate and hydrology of the Mediterranean Sea was strongly influenced by changes in precessional insolation, resulting in cyclic depositions of sapropels (black shales), particularly in the eastern basins [72, 114–119] (Fig 10D). The sapropels were formed during Northern Hemisphere summer insolation maxima/high-amplitude precession minima, resulting in an enhanced African monsoon [120], and as a consequence, in higher freshwater influx, particularly via the Nile River, and higher export productivity into the Mediterranean when compared to today (review in Rohling et al. [74]). To determine whether or not these orbital-driven changes in hydrology and paleo-productivity (Fig 10D) may have an influence on our reconstructions, we evaluated our reconstructions and orbital precession [73] for multicollinearity [121]. We allowed to find the optimal time lag between both records and found a slight but significant Pearson product-moment correlation between both records with r = 0.271 and p<0.001, suggesting that changes in the precession band indeed influence our reconstructions. This is not surprising because infaunal species are abundant and show marked fluctuations throughout the Pefka E section, suggesting a response to precession-related changes in paleo-productivity (Fig F in S3 File). For example, the most dominant species U. peregrina in the middle and upper parts of the Pefka E section has been found in high relative abundance in a middle to late Holocene sediment record, reflecting enhanced seasonal phytodetritus fluxes related to Nile river nutrient plumes [103]. As a consequence, we applied a moving average smoothing using a calculated average precession period of 20.5 kyrs to filter out the precession-driven component from our reconstructions (Fig 10A). The precession-corrected paleo-water depth estimates reveal two marked transgressive and regressive trends with water depth changes of 190 to 410 m between 1652 and 1546 ka, followed by relatively stable paleo-water depths of around 190 m between 1544 and 1531 ka, and two distinct transgressive phases with maximum paleo-water depths of 610 and 650 m, separated by one short-time regression phase with a minimum paleo-water depth of 290 m between 1531 and 1369 ka (Fig 10A). The upper part of the section is characterized by two regressive phases with minimum paleo-water depths of 417 and 249 m between 1368 and 1316 ka, separated by a short-term transgressive phase with a maximum paleo-water depth of 560 m. Paleo-water depths then slightly fluctuate between 250 and 340 m between 1316 and 1285 ka.
A. Original paleo-water depth reconstructions (excluding those samples with extremely poor fit to water depth) as thin line with sample specific error (error bars) and moving average of paleo-water depths with window size matching average precession period. For reference, global ice volume transferred into global sea level in meter relative to present [125] is given. B. Paleo-water depth estimates based on planktic/(benthic−stress markers) (P/(B−S)) ratios using the transfer function from Van der Zwaan et al. [123], with the thin green line showing the original reconstruction and the thick green line showing a five point moving average. C. Benthic foraminifera stable oxygen isotope ratios (δ18O). D. Comparison of the original paleo-water depth reconstructions (light blue) versus precession (dark blue) [73]; note inverse scaling of the precession curve. Symbols show Pleistocene sapropels identified in ODP Site 967 and ODP Hole 969D (Eastern Mediterranean Sea) [72].
By using the smoothed paleo-water depth curve as a basis, we compared our estimates to paleo-water depth estimates based on the relative abundance of planktonic foraminifera on the total foraminiferal fauna (P/B ratio) counted in the >125 μm fraction of the Pefka E section. Specifically, we used a modified ratio by excluding infaunal species (stress markers S) as described by van Hinsbergen et al. [122] and applied the formula given in Van der Zwaan et al. [123] (Fig 10B). The paleo-water depths (smoothed using a five-point moving average) based on the P/(B−S) ratios fluctuate between 220 and 780 m from 1652 to 1530 ka, suggesting a general increase in paleo-water depth, between 600 and 1080 m from 1527 to 1360 ka, and between 630 and 950 m from 1359 to 1297 ka. It has been shown that paleo-water depth estimates based on P/B ratios do not always give reasonable quantitative results due to specific local environmental conditions such as exposure to the open ocean, surface currents and trophic conditions [124]. However, for the Pefka E section the long-term transgressive and regressive trends of our TF based reconstructions are also reflected by the reconstructions based on the P/(B−S) ratios although absolute values of the P/(B−S) ratios are generally higher.
During the time of the deposition of the Pefka E section, orbital-scale eustatic sea-level changes can account for a maximum fluctuation of 85 m (Fig 10A; [125, 126]), suggesting that the long-term and corrected relative sea-level trend observed in the Pefka E section, with a maximum difference of about 470 m, is caused by substantial local to regional vertical tectonic movements, particularly between 1531 and 1316 ka. This interpretation is supported by the epibenthic δ18O record of this section that basically retraces the long-term relative sea-level trend (Fig 10C). The δ18O values range between 1.32 and 3.44 ‰ (on average 1.90 ‰) from 1620 to 1481 ka, and reveal an increase with maximum values of 3.72 and 4.42 ‰ (average 2.95 ‰) from 1478 to 1351 ka, followed by a decrease to values between 1.40 and 3.04 ‰ (average 1.95 ‰) from 1349 to 1271 ka. Shifts to higher water depths are accompanied by calcification under lower ambient bottom water temperatures. In the modern Eastern Mediterranean Sea, there is a decrease in temperature of ~0.5–0.7°C per 100 m [127] which would account for 0.8 to 1 ‰ in δ18O when paleo-water depth increased by ~400 m [128]. Additional effects of global ice-volume changes of around 1 ‰ for a maximum global sea-level change of 85 m during the deposition of the section [125] and potential salinity changes of Eastern Mediterranean surface and intermediate waters [129, 130] add to the temperature influence and likely account for the high-amplitude δ18O shifts of more than 2 ‰, particularly in the middle and upper parts of the section. However, it should be noted that we cannot exclude the potential influence of diagenetic effects on the stable isotope signal although the measured tests did not exhibit visible signs of dissolution and presence of cements.
For further regional evaluations of our reconstructed relative sea-level history for the Pefka E section, we compared the range of our corrected estimates with reported paleo-water depth ranges from other sediment sections on the island of Rhodes. According to Hanken et al. [63], the Kolymbia Formation reflects a deepening upward indicated by shallow-water benthic foraminiferal faunas (dominated by Elphidium spp. and Ammonia spp.) in the lower part and deeper water assemblages with presence of Uvigerina spp. and Bolivina spp. in the uppermost part of the formation. This deepening is also reflected by ostracod faunas, which revealed water depths of >250 m at the boundary between the Kolymbia Limestone and the LBF (Hanken et al. [63], and references therein). These values are in accordance with our estimates of 260 to 290 m for the lowermost part of the LBF. Moissette and Spjeldnæs [68] reported an increase in water depth across the LBF, mainly based on bryozoan zonations, with water depths ranging between 250 to 450 m and a maximum value of ~600 m for the Lindos Bay, Vasfi and Cape Vagia sections (Fig 1B). Finally, water depths ranging between 320 and 500 m based on benthic foraminiferal and ostracod assemblages for the LBF deposition at the Kallithea section in the NE part of Rhodes were reported [18, 131] (Fig 1B). These observations are in general accordance to our estimates, indicating a deepening with maximum water depths of 610 to 650 m in the middle to upper part of the section.
Although it is not the primary focus of our study, the reconstructed bathymetry of the Pefka E section presents some important implications for the regional development of the island: (i) When comparing the onset of the LBF deposition at Pefka E section with other localities [66], and especially the new chronostratigraphy of some sections provided by Quillévéré et al. [71], the onset of the deposition seems rather synchronous in the southern part of the island. Only at Kalithea the onset of LBF deposition potentially occurred subsequent to the phase of maximum transgression within the Pefka E section [132] (compare with Cornée et al. [133] for chronostratigraphic inconsistencies). This might imply a differential paleo-bathymetric development of northern Rhodes compared to central and southern Rhodes. (ii) The phase of maximum transgression within the Pefka E section between 1.5 and 1.35 Ma might imply that the it was either reached earlier than suggested by Titschack et al. [66] or Cornée et al. [133], or that it lasted much longer (until about 900 ka) and was interrupted by a short phase of uplift (represented in the upper part of the Pefka E section). The regressive phase during LBF deposition, reflecting sustained uplift of Rhodes, is preserved in the Lardos section (Fig 1B) where a shallowing from ~320 to <100 m water depth towards the top of the LBF has been observed [66].
The good agreement of the long-term water depth trend for the Pefka E section with available semi-quantitative reconstructions for other sediment sections of Rhodes confirms the general reliability of our reconstructions. Accordingly, we provide a first quantitative reconstruction of long-term relative sea-level changes during the deposition of the LBF. Our results reflect substantial local to regional vertical tectonic movements during the deposition of the LBF and provide information on the exact timing of neotectonic events for the island of Rhodes.
Conclusions
We developed a regional TF model based on benthic foraminiferal assemblages, covering a wider water depth range than any other model before. For this, we used a combined data set of modern shelf to slope samples from the Western Mediterranean and Adriatic Seas. We found a relation of the modern assemblages to water depth but also to other factors such as substrate and food availability. For the developed TF, we observed a predictive accuracy (RMSEP) of about ~50 m over a water depth range of ~1200 m which is comparable to other TF performances. However, it should be noted that segment-wise RMSEPs indicate a better predictive accuracy of 9 m for shallower water depths (0–200 m) but a worse predictive accuracy of 130 m for deeper water depths (200–1216 m). By applying a randomization test, we found our model to be robust for water depths ranging between 20 and 848 m, indicating that paleo-water depth reconstructions within this range can be reliable if the assemblages are not strongly influenced by other environmental parameters. We applied the model to the early Pleistocene Pefka E sediment section situated in the SE of the island of Rhodes and tested whether the fossil species have a good fit to water depth in the modern data set. We observed that most fossil samples have a poor to very poor fit to water depth which can be explained by the lack of close modern analogues. We screened the fossil samples, resulting in exclusion of samples with extremely poor fit to water depth, and determined whether or not the reconstructions are robust by comparing with an alternative TF model (i.e., based on the modern analogue technique). This comparison showed that both reconstructions are comparable. Although a randomization transfer function (randomTF) test revealed that our reconstructions are significant at the 95% confidence level, and hence appear generally trustworthy, we found our reconstructions to be strongly influenced by precession-related changes in paleo-productivity in the Eastern Mediterranean Sea and consequently filtered out this precession component. Our precession-corrected paleo-water depth estimates are in general accordance with reconstructions based on modified P/B ratios and epibenthic δ18O data of the same section. In addition, our reconstructions are in good agreement with semi-quantitative paleo-water depth estimates from contemporaneous sediment records of the island of Rhodes. The reconstructed water depth changes reflect the presence of substantial long-term regional tectonic movements of Rhodes during the early Pleistocene.
Supporting information
S1 File. Taxonomic list of species used in this study.
https://doi.org/10.1371/journal.pone.0188447.s001
(XLSX)
S3 File. Additional figures and tables.
Evaluation of the reaction norm test, ordination analyses, abundances of species in the samples, ecological preferences of the species, transfer function evaluation.
https://doi.org/10.1371/journal.pone.0188447.s003
(DOCX)
Acknowledgments
We acknowledge Maurice Ballein, Franziska Schmidtke, and Benedikt Walker for processing of sediment samples from the Pefka E section and for counting some of them. We gratefully acknowledge Erlend Martini for determining the nannoplankton association and assigning them to nannoplankton zones. We thank the reviewers Bruce W. Hayward (Geomarine Research Auckland), Maria Triantaphyllou (University of Athens), and Brent Wilson (University of the West Indies) for their critical revision of an earlier draft; and the handling editor Fabrizio Frontalini (Università degli Studi di Urbino). We finally thank Richard Telford (University of Bergen) for an additional and voluntary review, helping to further improve our manuscript.
References
- 1.
Imbrie J, Kipp NG. A new micropaleontological method for quantitative paleoclimatology: Application to a late Pleistocene Caribbean core. In: Turekian KK, editor. The Late Cenozoic Glacial Ages. New Haven, Connecticut: Yale University Press; 1971. p. 71–181.
- 2. Koc Karpuz N, Schrader H. Surface sediment diatom distribution and Holocene paleotemperature variations in the Greenland, Iceland and Norwegian Sea. Paleoceanography. 1990; 5(4): 557–80.
- 3. Pisias NG, Mix AC. Spatial and temporal oceanographic variability of the eastern equatorial Pacific during the late Pleistocene: evidence from Radiolaria microfossils. Paleoceanography. 1997; 12(3): 381–93.
- 4. Pflaumann U, Jian Z. Modern distribution patterns of planktonic foraminifera in the South China Sea and western Pacific: a new transfer technique to estimate regional sea-surface temperatures. Mar Geol. 1999; 156: 41–83.
- 5. Peyron O, de Vernal A. Application of artificial neural networks (ANN) to high-latitude dinocyst assemblages for the reconstruction of past sea-surface conditions in Arctic and sub–Arctic seas. J. Quaternary Sci. 2001; 16(7): 699–709.
- 6. Nakagawa T, Tarasov PE, Nishida K, Gotanda K, Yasuda Y. Quantitative pollen-based climate reconstruction in central Japan: application to surface and Late Quaternary spectra. Quat Sci Rev. 2002; 21(18–19): 2099–113.
- 7. Kucera M, Weinelt M, Kiefer T, Pflaumann U, Hayes A, Weinelt M, et al. Reconstruction of sea-surface temperatures from assemblages of planktonic foraminifera: multi-technique approach based on geographically constrained calibration data sets and its application to glacial Atlantic and Pacific Oceans. Quat Sci Rev. 2005; 24(7–9): 951–98.
- 8. Charman DJ, Blundell A. A new European testate amoebae transfer function for palaeohydrological reconstruction on ombrotrophic peatlands. J Quaternary Sci. 2007; 22(3): 209–21.
- 9. Gehrels WR, Callard SL, Moss PT, Marshall WA, Blaauw M, Hunter J, et al. Nineteenth and twentieth century sea–level changes in Tasmania and New Zealand. Earth Planet Sci Lett. 2012; 315–316: 94–102.
- 10. Kotthoff U, Greenwood DR, McCarthy FMG, Müller-Navarra K, Prader S, Hesselbo SP. Late Eocene to middle Miocene (33 to 13 million years ago) vegetation and climate development on the North American Atlantic Coastal Plain (IODP Expedition 313, Site M0027). Clim Past. 2014; 10: 1523–39.
- 11. Avnaim-Katav S, Milker Y, Schmiedl G, Sivan D, Hyams-Kaphzan O, Sandler A, et al. Impact of eustatic and tectonic processes on the southeastern Mediterranean shelf during the last one million years: Quantitative reconstructions using a foraminiferal transfer function. Mar Geol. 2016; 376: 26–38.
- 12. Hohenegger J. Estimation of environmental paleogradient values based on presence/absence data: a case study using benthic foraminifera for paleodepth estimation. Palaeogeogr Palaeoclimatol Palaeoecol. 2005; 217: 115–30.
- 13. Pérez-Asensio JN, Aguirre J, Schmiedl G, Civis J. Messinian paleoenvironmental evolution in the lower Guadalquivir Basin (SW Spain) based on benthic foraminifera. Palaeogeogr Palaeoclimatol Palaeoecol. 2012; 326–328: 135–51.
- 14. Hayward BW. Foraminifera-based estimates of paleobathymetry using Modern Analogue Technique, and the subsidence history of the early Miocene Waitemata Basin. New Zeal J Geol Geop. 2004; 47: 749–67.
- 15. Hayward BW, Triggs CM. Using multi-foraminiferal-proxies to resolve the paleogeographic history of a lower Miocene, subduction-related, sedimentary basin (Waitemata Basin, New Zealand). J Foramin Res. 2016; 46(3): 285–313.
- 16.
Wright R. Neogene paleobathymetry of the Mediterranean based on benthic Foraminifers from DSDP Leg 42A. In: Hsü KL, Montadert L, editors. Initial Rep Deep Sea Drill Proj. XLII: Scripps Inst. Oceanogr; 1975. p. 837–46.
- 17. Spezzaferri S, Tamburini F. Paleodepth variations on the Eratosthenes Seamount (Eastern Mediterranean): sea level changes or subsidence? eEarth Discuss. 2007; 2: 115–32.
- 18. Rasmussen TL, Thomsen E. Foraminifera and paleoenvironment of the Plio-Pleistocene Kallithea Bay section, Rhodes, Greece: evidence for cyclic sedimentation and shallow-water sapropels. Cushman Foundation for Foraminiferal Research. 2005; Special Publication 39: 15–51.
- 19. Hutson WH. The Agulhas current during the late Pleistocene: Analysis of modern faunal analogs. Science. 1980; 207(4426): 64–6. pmid:17730815
- 20. Juggins S. Quantitative reconstructions in palaeolimnology: new paradigm or sick science? Quaternary Sci Rev. 2013; 64: 20–32.
- 21. Horton BP. The distribution of contemporary intertidal foraminifera at Cowpen Marsh, Tees Estuary, UK: implications for studies of Holocene sea-level changes. Palaeogeogr Palaeoclimatol Palaeoecol. 1999; 149: 127–49.
- 22. Gehrels WR. Using foraminiferal transfer functions to produce high-resolution sea-level records from salt-marsh deposits, Maine, USA. Holocene. 2000; 10(3): 367–76.
- 23. Gehrels WR, Marshall WA, Gehrels MJ, Larsen G, Kirby JR, Eiríksson J, et al. Rapid sea-level rise in the North Atlantic Ocean since the first half of the nineteenth century. Holocene. 2006; 16(7): 949–65.
- 24. Horton BP, Edwards RJ. Quantifying Holocene sea-level change using intertidal foraminifera: Lessons from the British Isles. Cushman Foundation for Foraminiferal Research. 2006; Special publication No. 40: 1–97.
- 25. Kemp AC, Horton BP, Culver SJ, Corbett DR, van de Plassche O, Gehrels WR, et al. Timing and magnitude of recent accelerated sea-level rise (North Carolina, United States). Geology. 2009; 37(11): 1035–8.
- 26. Leorri E, Gehrels WR, Horton BP, Fatela F, Cearreta A. Distribution of foraminifera in salt marshes along the Atlantic coast of SW Europe: Tools to reconstruct past sea-level variations. Quat Int. 2010; 221: 105–55.
- 27. Kemp AC, Hawkes AD, Donnelly JP, Vane CH, Horton BP, Hill TD, et al. Relative sea-level change in Connecticut (USA) during the last 2200 yrs. Earth Planet Sci Lett. 2015; 428: 217–29.
- 28. Kemp AC, Kegel JJ, Culver SJ, Barber DC, Mallinson DJ, Leorri E, et al. Extended late Holocene relative sea-level histories for North Carolina, USA. Quat Sci Rev. 2017; 160: 13–30.
- 29. Hawkes AD, Horton BP, Nelson AR, Hill DF. The application of intertidal foraminifera to reconstruct coastal subsidence during the giant Cascadia earthquake of AD 1700 in Oregon, USA. Quat Int. 2010; 122(1–2): 116–40.
- 30. Milker Y, Horton BP, Vane CH, Engelhart SE, Nelson AR, Witter RC, et al. Annual and seasonal distribution of intertidal foraminifera and stable carbon isotope geochemistry, Bandon Marsh, Oregon, USA. J Foramin Res. 2015; 45(2): 146–55.
- 31. Milker Y, Horton BP, Nelson AR, Engelhart SE, Witter RC. Variability of intertidal foraminiferal assemblages in a salt marsh, Oregon, USA. Mar Micropaleontol. 2015; 118: 1–16.
- 32. Scott DB, Medioli FS. Vertical zonations of marsh foraminifera as accurate indicators of former sea-levels. Nature. 1978; 272: 528–31.
- 33. Morigi C, Jorissen FJ, Fraticelli S, Horton BP, Principi M, Sabbatini A, et al. Benthic foraminiferal evidence for the formation of the Holocene mud-belt and bathymetrical evolution in the central Adriatic Sea. Mar Micropaleontol. 2005; 57: 25–49.
- 34. Rossi V, Horton BP. The application of subtidal foraminifera-based transfer function to reconstruct Holocene paleobathymetry of the Po Delta, northern Adriatic Sea. J Foramin Res. 2009; 39(3): 180–90.
- 35. Milker Y, Schmiedl G, Betzler C. Paleobathymetric history of the Western Mediterranean Sea shelf during the latest glacial period and the Holocene: Quantitative reconstructions based on foraminiferal transfer functions. Palaeogeogr Palaeoclimatol Palaeoecol. 2011; 307: 324–38. doi: j.palaeo.2011.05.031.
- 36. Milker Y, Wilken M, Schumann J, Sakuna D, Feldens P, Schwarzer K, et al. Sediment transport on the inner shelf off Khao Lak (Andaman Sea, Thailand) during the 2004 Indian Ocean tsunami and former storm events: evidence from foraminiferal transfer functions. Nat Hazards Earth Sys. 2013; 13: 3113–28.
- 37. Jorissen FJ. The distribution of benthic foraminifera in the Adriatic Sea. Mar Micropaleontol. 1987; 12: 21–48.
- 38. De Stigter HC, Jorissen FJ, Van der Zwaan GJ. Bathymetric distribution and microhabitat partitioning of live (Rose Bengal stained) benthic foraminifera along a shelf to bathyal transect in the southern Adriatic Sea. J Foramin Res. 1998; 28(1): 40–65.
- 39. De Rijk S, Troelstra SR, Rohling EJ. Benthic foraminiferal distribution in the Mediterranean Sea. J Foramin Res. 1999; 29(2): 93–103.
- 40. De Rijk S, Jorissen FJ, Rohling EJ, Troelstra SR. Organic flux control on bathymetric zonation of Mediterranean benthic foraminifera. Mar Micropaleontol. 2000; 40: 151–66.
- 41. Avnaim-Katav S, Hyams-Kaphzan O, Milker Y, Almogi-Labin A. Bathymetric zonation of modern shelf benthic foraminifera in the Levantine Basin, eastern Mediterranean Sea. J Sea Res. 2015; 99: 97–106.
- 42. Ter Braak CJF, Juggins S. Weighted averaging partial least squares regression (WA-PLS): an improved method for reconstructing environmental variables from species assemblages. Hydrobiologia. 1993; 269/270: 485–502.
- 43. Hill MO, Gauch HG. Detrended Correspondence Analysis: An improved ordination technique. Vegetatio. 1980; 42(1/3): 47–58.
- 44. Ter Braak CJF. Canonical Correspondence Analysis: a new eigenvector technique for multivariate direct gradient analysis. Ecology. 1986; 67(5): 1167–79.
- 45. Austin MP. On non-linear species response models in ordination. Vegetatio. 1976; 33(1): 33–41.
- 46. Austin MP. A silent clash of paradigms: some inconsistencies in community ecology. Oikos. 1999; 86(1): 170–8.
- 47. Oksanen J, Minchin PR. Continuum theory revisited: what shape are species responses along ecological gradients? Ecol Modell. 2002; 157(2–3): 119–29.
- 48. Austin M. Species distribution models and ecological theory: A critical assessment and some possible new approaches. Ecol Modell. 2007; 200(1): 1–19.
- 49. Graham JH, Duda JJ. The humpbacked species richness-curve: A contingent rule for community ecology. Int J Ecol. 2011; 2011(868426): 15 p.
- 50. Yee TW, Mitchell ND. Generalized additive models in plant ecology. J Veg Sci. 1991; 2(5): 587–602.
- 51. Huntley B, Green RE, Collingham YC, Hill JK, Willis SG, Bartlein PJ, et al. The performance of models relating species geographical distributions to climate is independent of trophic level. Ecol Lett. 2004; 7(5): 417–26.
- 52. Wamelink GWW, Goedhart PW, Van Dobben HF, Berendse F, Díaz S. Plant species as predictors of soil pH: Replacing expert judgement with measurements. J Veg Sci. 2005; 16(4): 461–70.
- 53. Elith J , H. Graham C, P. Anderson R, Dudík M, Ferrier S, Guisan A, et al. Novel methods improve prediction of species’ distributions from occurrence data. Ecography. 2006; 29(2): 129–51.
- 54. Telford RJ, Birks HJB. Effect of uneven sampling along an environmental gradient on transfer-function performance. J Paleolimnol. 2011; 46(1): 99–106.
- 55. Telford RJ, Birks HJB. The secret assumption of transfer functions: problems with spatial autocorrelation in evaluating model performance. Quat Sci Rev. 2005; 24: 2173–9.
- 56. Telford RJ, Birks HJB. Evaluation of transfer functions in spatially structured environments. Quat Sci Rev. 2009; 28: 1309–16.
- 57. Birks HJB. Numerical tools in palaeolimnology—Progress, potentialities, and problems. J Paleolimnol. 1998; 20: 307–32.
- 58.
Simpson GL, Hall RI. Human Impacts: Applications of numerical methods to evaluate surface-water acidification and eutrophication. In: Birks BHJ, Lotter FA, Juggins S, Smol PJ, editors. Tracking Environmental Change Using Lake Sediments: Data Handling and Numerical Techniques. Dordrecht: Springer Netherlands; 2012. p. 579–614.
- 59. Telford RJ, Birks HJB. A novel method for assessing the statistical significance of quantitative reconstructions inferred from biotic assemblages. Quat Sci Rev. 2011; 30(9–10): 1272–8.
- 60. Jorissen FJ. Benthic foraminifera from the Adriatic Sea: principles of phenotypic variation. Utrecht Micropaleontological Bulletins. 1988; 37: 1–174.
- 61. Milker Y, Schmiedl G, Betzler C, Römer M, Jaramillo-Vogel D, Siccha M. Distribution of Recent benthic foraminifera in neritic carbonate environments of the Western Mediterranean Sea. Mar Micropalaeontol. 2009; 73: 207–25.
- 62. Meulenkamp E, de Mulder EFJ, Van de Weerd A. Sedimentary history and paleogeography of the late Cenozoic of the island of Rhodes. Zeitschrift der Deutschen Geologischen Gesellschaft. 1972; 123: 541–53.
- 63. Hanken N-M, Bromley RG, Miller J. Plio-Pleistocene sedimentation in coastal grabens, north-east Rhodes, Greece. Geol J. 1996; 31(4): 393–418. <393::AID-GJ712>3.0.CO;2-H.
- 64.
Schlitzer R. Ocean Data View. http://odv.awi.de; 2017.
- 65. Mutti E, Orombelli G, Pozzi R. Geological studies on the Dodecanese Islands (Aegean Sea). IX. Geological map of the island of Rhodes (Greece), explanatory notes. Annales de Géologie des Pays Helléniques. 1970; 22: 77–226.
- 66. Titschack J, Joseph N, Fietzke J, Freiwald A, Bromley RG. Record of a tectonically-controlled regression captured by changes in carbonate skeletal associations on a structured island shelf (mid-Pleistocene, Rhodes, Greece). Sediment Geol. 2013; 283: 15–33.
- 67. Titschack J, Bromley RG, Freiwald A. Plio-Pleistocene cliff-bound, wedge-shaped, warm-temperate carbonate deposits from Rhodes (Greece): Sedimentology and facies. Sediment Geol. 2005; 180(1–2): 29–56.
- 68. Moissette P, Spjeldnæs N. Plio–Pleistocene deep-water bryozoans from Rhodes, Greece. Palaeontology. 1995; 38: 771–99.
- 69. Hansen KS. Development of a prograding carbonate wedge during sea level fall: Lower Pleistocene of Rhodes, Greece. Sedimentology. 1999; 46: 559–76.
- 70. Frydas D. Die Pliozän/Pleistozän-Grenze auf der Insel Rhodos (Griechenland). Münstersche Forschungen zur Geologie und Paläontologie. 1994; 76: 331–44.
- 71. Quillévéré F, Cornée J-J, Moissette P, López-Otálvaro GE, van Baak C, Münch P, et al. Chronostratigraphy of uplifted Quaternary hemipelagic deposits from the Dodecanese island of Rhodes (Greece). Quat Res. 2016; 86(1): 79–94.
- 72.
Lourens LJ, Hilgen FJ, Raffi I. 15. Base of large Gephyrocapsa and astronomical calibration of early Pleistocene sapropels in Site 967 and Hole 969D: solving the chronology of the Vrica section (Calabria, Italy). In: Robertson AHF, Emeis K-C, Richter C, Camerlenghi A, editors. Proceeding of the Ocean Drilling Program, Scientific Results. 1998; 160: 191–7.
- 73. Laskar J, Robutel P, Joutel F, Gastineau M, Correia ACM, Levrard B. A long term numerical solution for the insolation quantities of the Earth. Astron Astrophys. 2004; 428(1): 261–85.
- 74. Rohling EJ, Marino G, Grant KM. Mediterranean climate and oceanography, and the periodic development of anoxic events (sapropels). Earth-Sci Rev. 2015; 143: 62–97.
- 75. Paillard D, Labeyrie L, Yiou P. Macintosh program performs time-series analysis. EOS, Transactions, American Geophysical Union. 1996; 77: 379.
- 76. Sgarrella F, Moncharmont Zei M. Benthic foraminifera in the Gulf of Naples (Italy): systematics and autoecology. Bollettino della Societa Paleontologica Italiana. 1993; 32(2): 145–264.
- 77. Cimerman F, Langer MR. Mediterranean foraminifera. Ljubljana: Slowenska Akademija Znanosti in Umetnosti. Academia Sciertiarum et Artium Slovenca Cl 4 Hist. Nat., 30; 1991. 118 p.
- 78.
Jones RW. The Challenger Foraminifera, 1. Edition. Oxford University Press Inc, New York. 1994. 149 p.
- 79. Rasmussen TL. Systematic paleontology and ecology of benthic foraminifera from the Plio-Pleistocene Kallithea Bay section, Rhodes, Greece. Cushman Foundation Special Publication. 2005; 39: 53–157.
- 80. Milker Y, Schmiedl G. A taxonomic guide to modern benthic shelf foraminifera of the western Mediterranean Sea. Palaeontol Electron. 2012; 15(2,16A): 1–134.
- 81.
Aitchison J. The Statistical Analysis of Compositional Data. London: Chapman & Hall; 1986. 416 p.
- 82.
R-Development-Core-Team. R: a Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna. Available at: http://www.R-project.org/. 2016.
- 83.
Oksanen J, Blanchet FG, Kindt R, Legendre P, Minchin PR, O'Hara RB, et al. vegan: Community Ecology Package. R package version 2.3–3, http://CRAN.R-project.org/package=vegan. 2016.
- 84.
McDonald JH. Handbook of Biological Statistics. Baltimore: Sparky House Publishing; 2009. 313 p.
- 85.
Huber PJ. Robust Statistics. New York, Chichester, Brisbane, Toronto, Singapore: John Wiley & Sons, Ltd.; 1981. 308 p.
- 86. Akaike H. A new look at the statistical model identification. IEEE Trans Autom Control. 1974; 19(6): 716–23.
- 87. Barlow NLM, Shennan I, Long AJ, Gehrels WR, Saher MH, Woodroffe SA, et al. Salt marshes as late Holocene tide gauges. Glob Planet Change. 2013; 106: 90–110.
- 88.
Juggins S. rioja: Analysis of Quaternary Science Data. R package version 09–5, http://CRAN.R-project.org/package=rioja. 2015.
- 89. Kendall MG. A new measurement of rank correlation. Biometrika. 1938; 30: 81–93.
- 90.
Kemp AC, Telford RJ. Transfer functions. In: Shennan I, Long AJ, Horton BP, editors. Handbook of Sea-Level Research. Chichester: John Wiley and Sons; 2015. p. 470–99.
- 91. Simpson GL. Analogue methods in palaeoecology: Using the analogue package. J Stat Softw. 2007; 22(2): 1–29.
- 92. van den Wollenberg AL. Redundancy analysis–an alternative for canonical correlation analysis. Psychometrika. 1977; 42(2): 207–19.
- 93.
Telford RJ. palaeoSig: Significance tests for palaeoenvironmental reconstructions. R package, version 1.1–3, http://CRAN.R-project.org/package=palaeoSig. 2015.
- 94. Bartels-Jonsdottir HB, Knudsen KL, Abrantes F, Lebreiro S, Eiriksson J. Climate variability during the last 2000 years in the Tagus Prodelta, western Iberian Margin: Benthic foraminifera and stable isotopes. Mar Micropaleontol. 2006; 59: 83–103.
- 95. Schmiedl G, Mackensen A. Late Quaternary paleoproductivity and deep water circulation in the eastern South Atlantic Ocean: Evidence from benthic foraminifera. Palaeogeogr Palaeoclimatol Palaeoecol. 1997; 130: 43–80.
- 96. de Stigter HC, van der Zwaan GJ, Langone L. Differential rates of benthic foraminiferal test production in surface and subsurface sediment habitats in the southern Adriatic Sea. Palaeogeogr Palaeoclimatol Palaeoecol. 1999; 149(1–4): 67–88.
- 97. Schmiedl G, de Bovee F, Buscail R, Charriere B, Hemleben C, Medernach L, et al. Trophic control of benthic foraminiferal abundance and microhabitat in the bathyal Gulf of Lions, western Mediterranean Sea. Mar Micropaleontol. 2000; 40: 167–88.
- 98. Fontanier C, Jorissen FJ, Chaillou G, David C, Anschutz P, Lafon V. Seasonal and interannual variability of benthic foraminiferal faunas at 550 m depth in the Bay of Biscay. Deep Sea Res Part I: Oceanographic Research Papers. 2003; 50(4): 457–94.
- 99. Duchemin G, Jorissen FJ, Le Loc´h F, Andrieux-Loyer F, Hily C, Thouzeau G. Seasonal variability of living benthic foraminifera from the outer continental shelf of the Bay of Biscay. J Sea Res. 2008; 59: 297–319.
- 100. Goineau A, Fontanier C, Jorissen FJ, Lansard B, Buscail R, Mouret A, et al. Live (stained) benthic foraminifera from the Rhône prodelta (Gulf of Lion, NW Mediterranean): Environmental controls on a river-dominated shelf. J Sea Res. 2011; 65(1): 58–75.
- 101. Goineau A, Fontanier C, Mojtahid M, Fanget AS, Bassetti MA, Berné S, et al. Live–dead comparison of benthic foraminiferal faunas from the Rhône prodelta (Gulf of Lions, NW Mediterranean): Development of a proxy for palaeoenvironmental reconstructions. Mar Micropaleontol. 2015; 119: 17–33.
- 102. Dimiza MD, Triantaphyllou MV, Koukousioura O, Hallock P, Simboura N, Karageorgis AP, et al. The Foram Stress Index: A new tool for environmental assessment of soft-bottom environments using benthic foraminifera. A case study from the Saronikos Gulf, Greece, Eastern Mediterranean. Ecol Indic. 2016; 60: 611–21.
- 103. Schmiedl G, Kuhnt T, Ehrmann W, Emeis K-C, Hamann Y, Kotthoff U, et al. Climatic forcing of eastern Mediterranean deep-water formation and benthic ecosystems during the past 22000 years. Quat Sci Rev. 2010; 29: 3006–20.
- 104.
Juggins S, Birks HJB. Quantitative environmental reconstructions from biological data. In: Birks HJB, Lotter AF, Juggins S, Smol JP, editors. Tracking Environmental Change Using Lake Sediments–Data Handling and Numerical Techniques. Developments in Paleoenvironmental Research. Netherlands: Springer Science Business Media B.V; 2012. p. 431–94.
- 105. Hawkes AD, Horton BP, Nelson AR, Vane CH, Sawai Y. Coastal subsidence in Oregon, USA, during the giant Cascadia earthquake of AD 1700. Quat Sci Rev. 2011; 30: 364–76.
- 106. Engelhart SE, Horton BP, Nelson AR, Hawkes AD, Witter RC, Wang K, et al. Testing the use of microfossils to reconstruct great earthquakes at Cascadia. Geology. 2013; 41(10): 1067–70.
- 107. Milker Y, Nelson AR, Horton BP, Engelhart SE, Bradley L-A, Witter RC. Differences in coastal subsidence in southern Oregon (USA) during at least six prehistoric megathrust earthquakes. Quat Sci Rev. 2016; 142: 143–63.
- 108. Parker FL. Eastern Mediterranean Foraminifera. Reports of the Swedish Deep-Sea Expedition, Göteburg, Sweden. 1958; 8(4): 217–83.
- 109. Dimiza MD, Koukousioura O, Triantaphyllou MV, Dermitzakis MD. Benthic foraminiferal assemblages from coastal environments of the Aegean Sea (Greece). Rev Micropaleontol. 2016; 59: 19–32.
- 110.
Ter Braak CJF, Juggins S, Birks HJB, Van der Voet H. Weighted Averaging Least Squares regression (WA-PLS): definition and comparison with other methods for species-environment calibration. In: Patil GP, Rao CR, editors. Multivariate environmental statistics. Amsterdam: Elsevier Science Publishers B.V. (North Holland); 1993. p. 525–60.
- 111. Ter Braak CJF. Non-linear methods for multivariate statistical calibration and their use in palaeoecology: a comparison of inverse (k-nearest neighbours, partial least squares and weighted averaging partial least squares) and classical approaches. Chemometr Intell Lab. 1995; 28: 165–80.
- 112. Hutson WH. Transfer functions under no-analog conditions: Experiments with Indian Ocean planktonic Foraminifera. Quat Res. 1977; 8(3): 355–67.
- 113. Payne RJ, Babeshko KV, van Bellen S, Blackford JJ, Booth RK, Charman DJ, et al. Significance testing testate amoeba water table reconstructions. Quat Sci Rev. 2016; 138: 131–5.
- 114. Hilgen FJ. Astronomical calibration of Gauss to Matuyama sapropels in the Mediterranean and implication for the Geomagnetic Polarity Time Scale. Earth Planet Sci Lett. 1991; 104(2): 226–44.
- 115.
Kroon D, Alexander I, Little M, Lourens LJ, Matthewson A, Robertson AHF, et al. Oxygen isotope and sapropel stratigraphy in the eastern Mediterranean during the last 3.2 million years. In: Robertson AHF, Emeis K-C, Richter C, Camerlenghi A, editors. Proceedings of the Ocean Drilling Program, Scientific Results. 160; 1998. p. 181–9.
- 116. Tuenter E, Weber SL, Hilgen FJ, Lourens LJ. The response of the African summer monsoon to remote and local forcing due to precession and obliquity. Glob Planet Change. 2003; 36: 219–35.
- 117. Ziegler M, Tuenter E, Lourens LJ. The precession phase of the boreal summer monsoon as viewed from the eastern Mediterranean (ODP Site 968). Quat Sci Rev. 2010; 29(11–12): 1481–90.
- 118. Athanasiou M, Triantaphyllou MV, Dimiza MD, Gogou A, Τheodorou G. Zanclean/Piacenzian transition on Cyprus (SE Mediterranean): calcareous nannofossil evidence of sapropel formation. Geo–Mar Lett. 2015; 35(5): 367–85.
- 119. Athanasiou M, Bouloubassi I, Gogou A, Klein V, Dimiza MD, Parinos C, et al. Sea surface temperatures and environmental conditions during the “warm Pliocene” interval (4.1–3.2 Ma) in the Eastern Mediterranean (Cyprus). Glob Planet Change. 2017; 150: 46–57.
- 120. Rossignol-Strick M, Nesteroff W, Olive P, Vergnaud-Grazzini C. After the deluge: Mediterranean stagnation and sapropel formation. Nature. 1982; 295: 105–10.
- 121. Dormann CF, Elith J, Bacher S, Buchmann C, Carl G, Carré G, et al. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography. 2013; 36(1): 27–46.
- 122. Van Hinsbergen DJJ, Kouwenhoven TJ, van der Zwaan GJ. Paleobathymetry in the backstripping procedure: Correction for oxygenation effects on depth estimates. Palaeogeogr Palaeoclimatol Palaeoecol. 2005; 221: 245–65.
- 123. Van der Zwaan GJ, Jorissen FJ, de Stigter HC. The depth dependency of planktonic/benthic foraminiferal ratios: Constraints and applications. Mar Geol. 1990; 95: 1–16.
- 124.
Milker Y. Western Mediterranean shelf foraminifera: Recent distribution, Holocene sea-level reconstructions, and paleoceanographic implications. Hamburg. PhD thesis, University of Hamburg. 2010. Available from: http://ediss.sub.uni-hamburg.de/volltexte/2010/4709/
- 125. Bintanja R, van de Wal RSW. North American ice-sheet dynamics and the onset of 100,000-year glacial cycles. Nature. 2008; 454(7206): 869–72. pmid:18704083
- 126. Rohling EJ, Foster GL, Grant KM, Marino G, Roberts AP, Tamisiea ME, et al. Sea-level and deep-sea-temperature variability over the past 5.3 million years. Nature. 2014; 508: 477–82. pmid:24739960
- 127. Tanhua T, Hainbucher D, Schröder K, Cardin V, Álvarez M, Civitarese G. The Mediterranean Sea system: a review and an introduction to the special issue. Ocean Sci. 2013; 9: 789–803.
- 128. Wefer G, Berger WH. Isotope paleontology: growth and composition of extant calcareous species. Mar Geol. 1991; 100(1): 207–48.
- 129.
Rohling EJ, Cooke S. Stable oxygen and carbon isotopes in foraminiferal carbonate shells. In: Sen Gupta BK, editor. Modern Foraminifera. London: Kluwer Academic Press; 1999. p. 239–58.
- 130. Rohling EJ, Bigg GR. Paleosalinity and δ18O: A critical assessment. J Geophys Res: Oceans. 1998; 103(C1): 1307–18.
- 131. Hastrup A, Thomsen E. Paleoenvironmental interpretation of the Plio-Pleistocene Kallithea Bay section, Rhodes, Greece based on ostracods. Cushman Foundation Special Publication. 2005; 39: 159–91.
- 132. Thomsen E, Rasmussen TL, Hastrup A. Calcareous nannofossil, ostracode and foraminifera biostratigraphy of Plio-Pleistocene deposits, Rhodes (Greece), with a correlation to the Vrica section (Italy). J Micropalaeontol. 2001; 20(2): 143–54.
- 133. Cornée J-J, Münch P, Quillévéré F, Moissette P, Vasiliev I, Krijgsman W, et al. Timing of Late Pliocene to Middle Pleistocene tectonic events in Rhodes (Greece) inferred from magneto-biostratigraphy and 40Ar/39Ar dating of a volcaniclastic layer. Earth Planet Sci Lett. 2006; 250(1): 281–91.