Disentangling environmental drivers of subarctic dinocyst assemblage compositional change during the Holocene

Analysis of compositional changes in fossil organic-walled dinoflagellate cyst (dinocysts) assemblages is a widely used tool for the quantitative reconstruction of past environmental parameters. This approach assumes that the assemblage composition is significantly related to the reconstructed environmental parameter, which requires an independent correlation between the assemblage and the variable, meaning that the variable explains a dimension of the assemblage variance that is not also explained by other parameters. Theoretically, dinocyst assemblages can be used to reconstruct multiple environmental variables. However, there is evidence that primary and subordinate drivers for assemblage compositions regionally differ, and it remains unclear whether a significant relationship to specific parameters in the present ocean always implies that this relationship is significant in fossil assemblages, questioning if past changes in these multiple parameters can be reconstructed simultaneously from fossil assemblages. Here, we analysed a local subset of the Northern Hemisphere dinocyst calibration dataset (n=1968), including samples from the Baffin Bay area (n=421), and evaluated the benefits of a local versus a more regional or global calibration for the environmental reconstruction of Baffin Bay oceanography during the Holocene. We determined the dimensionality of the dinocyst ecological response and identified environmental drivers in the Baffin Bay area for the modern dataset. We analysed four existing Holocene records along a north–south transect in the area and evaluated the statistical significance of downcore reconstructions by applying the local and global datasets with different techniques: the modern analogue technique (MAT), the weighted average partial least square (WA-PLS) and maximum likelihood (ML). We evaluated reconstructions tested as significant in light of the existing state of knowledge about West Greenland's Holocene palaeoceanography. Our analyses imply that present-day and Holocene dinocyst assemblages in the Baffin Bay area are primarily driven by salinity changes; other parameters were less important in driving assemblage compositions, and their contribution differed among the studied records. We show that the objectively occurring and temporally coherent shifts in dinocyst assemblages in the Holocene of Baffin Bay can be robustly interpreted only by transfer functions that are locally calibrated. Transfer functions based on the broad Northern Hemisphere calibration yielded many insignificant environmental reconstructions. At the same time, we show that even in the local calibration, not all parameters that appear to significantly affect dinocyst assemblages in the calibration dataset can be meaningfully reconstructed in the fossil record. A thorough evaluation of the calibration dataset and the significance of downcore applications is necessary to reveal the region-specific information contained in dinocyst assemblage composition.

1 Introduction

During the last decades, the growing interest in environmental responses to changing climate produced many studies documenting climate variations based on reconstructions making use of marine microfossils (e.g. Muller et al., 1983; Sarnthein et al., 1988; Berger et al., 1989; Meyers, 1997; Rühlemann et al., 1999; Fischer et al., 2000). Next to analyses of chemical and isotopic signals locked in microfossil skeletons, quantitative analyses of the changing taxonomic composition of microfossil assemblages, using ecological models, have been common in palaeoceanography since Imbrie and Kipp (1971). This approach assumes that the composition of marine plankton and benthic communities reflects the environment in which they lived and that this information is stored within the sediment by their fossils. The hypothesis that an assemblage of taxa can be expressed as a function of (a) certain environmental variable(s) (ter Braak, 1987) allows for quantitative modelling applying transfer functions and analogue methods (from now on called transfer functions). Multiple microfossils like planktonic foraminifera (e.g. Kucera et al., 2005), diatoms (e.g. Sha et al., 2014) and dinoflagellate cysts (dinocysts) (e.g. de Vernal et al., 2005) have been used for environmental reconstruction applying multiple transfer function methods.

However, the key assumption of this approach requires that the assemblage composition is significantly related to the given environmental parameter to reconstruct both at present for the purpose of calibration and in the past for the purpose of reconstruction. In the presence of multiple environmental parameters affecting assemblage composition, it is challenging to disambiguate the effect of each parameter and determine which one(s) acted on the analysed fossil assemblages in the past.

There is evidence that the assemblages of zooplankton microfossils, such as planktonic foraminifera and radiolarians, being exclusively oceanic, are primarily controlled by temperature (Morey et al., 2005; Hernández-Almeida et al., 2017). Phytoplanktonic microfossil assemblages however, like diatoms and dinocysts, seem to be affected by multiple environmental factors (de Vernal et al., 2000, 2001; Dale et al., 2002; Holzwarth et al., 2007; Radi and de Vernal, 2008; Ribeiro and Amorim, 2008; Lopes et al., 2010; Hohmann et al., 2020). Apparently, phytoplankton communities have a complex relationship to the environment and independently reflect the attribute of hydrography and nutrient availability, rendering the possibility of multiple environmental factor reconstruction. It is tempting to explore whether such dependency on multiple factors can be used to reconstruct changes in all the driving factors in the past.

About 20 % of dinoflagellate species produce organic-walled cysts within their life cycles (Wall and Dale, 1968; Dale, 1983; Taylor and Pollingher, 1987; Head, 1996). As the walls of these cysts consist of a condensed, predominantly aromatic macromolecular structure (Bogus et al., 2014), they are relatively well preserved in seafloor sediments (e.g. Zonneveld et al., 2008), not being prone to dissolution after death, unlike mineralised fossils such as calcareous (e.g. foraminifera and nannoplankton) or siliceous (e.g. diatoms and radiolarians) microfossils (Koç et al., 1993; Matthiessen et al., 2001; Seidenkrantz et al., 2007; Schröder-Adams and van Rooyen, 2011; Zamelczyk et al., 2012). However, during exposure on the seafloor, oxidation of the dinocyst may result in selective organic matter degradation, which might affect the assemblage composition (Zonneveld et al., 2010). Nevertheless, if organic matter degradation is negligible, in settings with carbonate or opal dissolution, dinocyst-based transfer functions are often the best means to reconstruct past surface ocean properties.

Dinoflagellates form a major group of planktonic marine primary producers and include diverse autotrophic, mixotrophic and heterotrophic forms that all thrive in the upper layers of the oceans (e.g. de Vernal and Marret, 2007). Dinocyst studies of marine surface sediment samples revealed that dinocysts are excellent tracers of hydrological and biological conditions (Mudie, 1992; Rochon and Vernal, 1994; Matthiessen, 1995; de Vernal et al., 1997, 2000, 2001, 2013; Radi et al., 2001; Radi and de Vernal, 2004, 2008; Rochon et al., 2008). While dinocyst assemblages can theoretically be applied to quantify multiple environmental drivers, identification of the primary driving factor(s) for the dinocyst assemblage compositional changes is challenging. Relationships between the assemblages and environmental parameters are complex, pointing to regionally different relationships between assemblages and environmental conditions. For example, whereas temperature and salinity gradients are determinants in the Atlantic–Arctic region, primary productivity seems to explain a large portion of taxa variability in the Pacific (Radi and de Vernal, 2008; Hohmann et al., 2020). Studies suggest that the number and type of primary factors driving dinocyst assemblages vary regionally (Zonneveld and Siccha, 2016; Hohmann et al., 2020). Additionally, the evaluation of the performance of transfer functions seems to be method-specific (Hohmann et al., 2020) and biased by spatial autocorrelation (Telford and Birks, 2005, 2009).

Hohmann et al. (2020) analysed the dimensionality of the dinocyst environmental response separately for assemblages from the Pacific and the Atlantic–Arctic regions, calibrated and evaluated multiple transfer function methods for predicting environmental variables, and estimated their performance in light of spatial autocorrelation. This study identified different primary drivers for two areas affecting dinocyst assemblage composition, highlighting the merit of local calibrations and the necessity to carry out variable selections of main driving factors for each study area separately. Defining a local calibration dataset has the advantage of reducing the complexity of the assemblage to environment composition. However, this comes at the cost of reducing the pool of compositional variability, leading to the risk of biasing the results. For example, in Holocene cores off West Greenland, some of the best compositional analogues have been found in the Gulf of St. Lawrence (e.g. Allan et al., 2018).

However, because factor(s) driving microfossil assemblages differ between regions, it is possible that the same fossil assemblages may find equally good compositional analogues representing different combinations of environmental parameters. This is not to say that a specific taxon responds differently to a given environmental parameter in different regions (the behaviour of a taxon within the environmental gradient is global), but as the importance of environmental parameters differs regionally due to a higher variance of certain parameters, the taxa assemblages react mostly to these variables. This may lead to situations where the assemblages are the same, although some parameters in the environment were different. Hence, while local calibrations limit the range of compositional variability available to derive the transfer function, they are more likely to disambiguate the main driving variable in a certain region, which may render reconstructions in that region more robust assuming the driving factors remain the same.

By calibrating locally, we can exclude areas where dinocysts are susceptible to lateral transport during sinking and sedimentation processes. There is evidence that a dinocyst could be displaced up to thousands of kilometres from the location at the surface where it was produced (Nooteboom et al., 2019). Long-distance cyst transport by bottom waters and sediment flows was also suggested based on sediment trap studies in the central and North Atlantic Ocean (Dale, 1992; Dale and Dale, 1992). This implies that in areas with strong ocean currents, the surface above the sedimentation area does not necessarily mirror the environmental conditions under which the plankton was produced, generating a nuisance in the calibration dataset.

The aim of this study is to investigate if a local dinocyst-based transfer function calibration on samples covering the Baffin Bay and surrounding regions allows a robust disambiguation of the factors that affected the assemblage variation in this region during the Holocene. Baffin Bay is a climate-sensitive region characterised by strong seasonal sea ice variability and meltwater discharge from the Greenland Ice Sheet (GIS). In view of the current rapid warming in the Arctic, it is relevant to understand how this marine basin reacted to climate changes during the Holocene. Indeed, to this end, many Holocene reconstructions based on fossil dinocyst assemblages were made and could serve as a basis for a comprehensive analysis of factors affecting the assemblage composition (e.g. Levac et al., 2001; Ouellet-Bernier et al., 2014; Allan et al., 2018; Caron et al., 2019).

We begin by carrying out an objective variable selection, guided by canonical ordination, and investigating the data structure of the resulting designed local transfer functions. We proceed by testing their ability to interpret assemblage changes recorded in Holocene sediments by their applications to four records covering a north–south transect throughout the Baffin Bay area. To this end, we used data from three existing records and added new data from a core located in southern Melville Bay (see Fig. 1 and Table 1). We specifically focus on Holocene sections, assuming that the magnitude of environmental change in the Holocene was relatively small for the local calibration to be sufficient to capture the range of past environments recorded in the cores. To estimate the significance of the reconstructions we apply a test designed by Telford and Birks (2011), evaluating if reconstructions explain more of the variation in the fossil taxa data than random environmental variables assigned to the same fossil data. To assess the merit of a local versus a larger regional calibration, we compared local transfer function results to regional transfer functions based on a Northern Hemisphere calibration dataset.

Table 1 List of cores used in this study, including core location and water depth (m) and the modern sea surface conditions – sea surface temperature (SST) in summer, sea surface salinity (SSS) in summer, months per year of sea ice cover (SeaIce) and annual productivity (PPannual) of organic carbon – for each location provided by the WOA (World Ocean Atlas; Locarnini et al., 2013; Zweng et al., 2013) and NSIDC (National Snow and Ice Data Center; Walsh et al., 2015)

Figure 1 Map of the study area. (a) Schematic modern surface water circulation in the Baffin Bay area and around Greenland (BIC – Baffin Island Current, WGC – West Greenland Current, EGC – East Greenland Current, IC – Irminger Current). (b) Location of surface sediment samples (red dots) from the calibration dataset and cores (orange dots; see Table 1) used in this study. On the right, modern environmental parameter gradients gridded using DIVA gridding with Ocean Data View (Schlitzer, 2018). Units are listed in Table 3. The 6-month sea ice edge in the sea ice gradient map is indicated by the thick black line.

2 Material and methods

2.1 Data

The calibration dataset includes a local subset of the data from the Northern Hemisphere dataset (n=1968) in de Vernal et al. (2020). The local subset comprises census counts of recent organic-walled dinoflagellate cyst assemblages from core-top sediments at 421 sites from Baffin Bay with Kane Basin, Norwegian Bay and Lancaster Sound (n=233); Davis Strait (n=103); Labrador Sea (n=30); and Hudson Bay with Hudson Strait and Foxe Basin (n=55) (Fig. 1b). In samples containing low dinocyst concentrations, minimum counts of 60 specimens were accepted. The core-top samples were collected in the uppermost 1 cm of box or gravity cores and represent a few tens of years to centuries, depending upon sedimentation rates and mixing due to bioturbation (Radi and de Vernal, 2008). Palynological sample processing followed the procedures for palynological preparation described in de Vernal et al. (2010), which exclude the use of oxidation techniques to avoid selective degradation of more sensitive cyst taxa such as those produced by Protoperidiniales (Marret, 1993). The Northern Hemisphere dataset from de Vernal et al.  (2020) also served as the regional reference dataset, which we use to compare the results of the local dataset to. With the exception of the samples from Price et al. (2016, 2018), which were added at a later stage, the Northern Hemisphere dataset and therefore all samples included in the local dataset tested negative for significant selective organic matter degradation (Zonneveld et al., 1997) in Hohmann et al. (2020). The taxonomy used corresponds to that presented in Rochon et al. (1999) and updated in van Nieuwenhove et al. (2020). The taxonomic resolution of the local dataset was reduced to 33 dinocyst taxa after grouping following Hohmann et al. (2020) (Table 2, list of taxa). The original dinocyst dataset holding taxon percentages is archived at PANGAEA (de Vernal et al., 2019); the local data matrices with reduced taxonomy used in this analysis are available in the Supplement (Dataset S1).

We chose to carry out a canonical-ordination-based parameter selection (Lopes et al., 2010) to characterise the relationship between environmental variables and dinocyst assemblage composition. We defined a set of potential environmental controlling variables for dinocyst distributions in the Baffin Bay region: surface layer (here represented by conditions at 0 m depth) temperature and salinity in summer and sea ice cover and productivity in spring and summer as well as for the whole year. We did not consider temperature, salinity and productivity parameters for the cold seasons since most of the area is covered by sea ice, and there is hardly any productivity and hence no dinocyst production. Apart from environmental factors, some life-cycle factors (e.g. Fensome et al., 1993) such as water depth and distance from shore seem to affect dinocyst assemblages. But as they represent a separate independent dimension of variance (Hohmann et al., 2020), we do not include them in the local calibration dataset. The initial variable selection is listed in Table 3 and includes some of the variables considered in Hohmann et al. (2020) (cf. Hohmann et al., 2020, for more detail). Their spatial distribution is shown in Fig. 1.

Table 2 List of dinocyst taxa in the local calibration dataset, including abbreviations and highest relative abundances (minimum is always zero). Asterisks indicate heterotrophic taxa.

Table 3 Environmental variables used in the canonical-ordination-based parameter selection and their ranges within the local calibration dataset.

Figure 2 Flowchart describing the approach and sequence of analyses in this study.

We used four sediment sequences to test the ability of the locally calibrated transfer functions to interpret assemblage changes during the Holocene. These sediment cores cover a north–south transect throughout the Baffin Bay region (Fig. 1b, Table 1), including intervals within the last 10 ka.

We analysed dinocyst assemblages in core GeoB19927-3 (73^∘35,26^′ N, 58^∘05,66^′ W), which was taken at 932 m water depth by gravity coring during cruise MSM44 in 2015 (Dorschel et al., 2015). This core is from southern Melville Bay, an area where the relatively warm high-salinity West Greenland Current (WGC) interacts with the cold polar water of the Baffin Current (BC) originating from the Arctic Ocean and meltwater from glaciers of the north-western Greenland region. It is located optimally and offers a high resolution and continuous sedimentation (Dorschel et al., 2015). A high-resolution radiocarbon-based chronology and organic biomarker proxies for sea ice and primary productivity have already been presented by Saini et al. (2020). The core consists of 1147 cm of sediment. The working half of the core was sampled and divided into four sets of samples. One set was freeze-dried for dinocyst analysis. The other three sets were stored for other multi-proxy analyses and chronology (i.e. biomarker, foraminifera and provenance studies). In reference to summer temperature anomalies during the Holocene for the Northern Hemisphere (Marcott et al., 2013; Kaufman et al., 2020), we expected a possible trend towards the last 2000 years and sampled at a higher resolution within the top 270 cm. We sub-sampled every centimetre within the top 25 cm, every 5 cm within the 25 to 280 cm interval and every 10 cm within the 280 to 760 cm interval, covering the last 7677 years and including 124 samples (see Table 1). Sample processing followed the procedure for palynological preparation described in de Vernal et al. (2010), which consists of repeated HCl and HF treatment. The taxonomy of dinocysts used here was based on Rochon et al. (1999) and de Vernal et al. (2020), considering all taxonomic categories that were resolved in the calibration dataset. To obtain statistically reliable assemblage counts, at least 300 dinocyst specimens were counted per sample when possible. For samples with low dinocyst abundance as many specimens as possible were enumerated. For the sample with the lowest dinocyst abundance, we counted 89 specimens. The raw assemblage dataset for GeoB19927-3 will be available on time via the PANGAEA database.

In addition, we re-evaluated three existing sediment records that have already been studied in detail and have been analysed for their palynological content and dinocyst assemblages. HU91-039-008P (Levac et al., 2001) was sampled every 20 cm from the top to 820 cm, which corresponds to 1549–6756 years BP and includes 41 samples. MSM343310 (Allan et al., 2018) includes 194 samples taken around every 8 cm from the top to 925 cm (135–3578 years BP), while HU2008-029-004 (Gibb et al., 2015) includes 30 samples from every 4 cm from the top to 116.5 cm (0–9936 years BP) (see Table 1).

2.2 Data analysis

To explore the local relationship between dinocyst assemblages and environmental variables, we followed the procedure schematically illustrated in Fig. 2 and used statistical ordination (McGarigal et al., 2000). We performed detrended correspondence analysis (DCA) (Birks, 1995) to evaluate the gradient length, expressed as standard deviation (SD) units in the first DCA axis to decide whether to apply linear (SD < 3) or unimodal ordination (SD > 4, SD > 3 to <4; both methods possible, but unimodal ordination is advised) (Lepš and Šmilauer, 2003). For DCA and further analyses, we used the vegan package (Oksanen et al., 2019) in R (R Core Team, 2017). The dataset was found to have an SD of 3.0 (Fig. S1 in the Supplement) using log-transformed (log (x+1)) dinocyst assemblage per mil data. As the SD of our dataset is biased towards the linear threshold, we proceeded with linear techniques.

In the first step, we selected the main drivers affecting compositional changes in the dinocyst assemblages. We applied redundancy analysis (RDA), a constrained ordination (McGarigal et al., 2000) for linear data, to quantify the strength of the relationship between species variances and the tested forcing parameters. We applied the Hellinger transformation (Rao, 1995) to the dinocyst taxa ratios (33 taxa with abundances in percent). This transformation was proposed by Legendre and Gallagher (2001) for assemblage data that contain many double-zero abundances, as it allows the use of it without considering the common absence of species as a resemblance between assemblages. The Hellinger transformation also gives less weight to species with low counts and many zeros, focussing more on the composition of common species, as we assume that local datasets unlike global ones do not include indicators or inter-regional endemic species. It also ensures the preservation of Euclidean distance in a linear ordination when unimodal data are used as our species data SD did not strongly indicate linear data. Environmental data were standardised as described above to ensure the compatibility of variables for comparisons.

RDA was performed after ter Braak (1986) using the platform R (R Core Team, 2017). For the variable selection, we followed the procedure described in Hohmann et al. (2020), which includes a variable selection considering variance inflation factors (VIFs) ≤2 with the objective to extricate independent variables that show no collinearity between each other and explain a significant amount of the variance in the species data. Following this procedure, we obtained a set of three independent forcing variables, meaning that each of these three parameters explains a separate dimension of variance in the taxa data (for more explanation see Script S1 in the Supplement). These three independent variables, which were revealed to be summer sea surface temperature (SSTsummer), summer sea surface salinity (SSSsummer) and spring productivity (PPspring) (see Sect. 3.2), were used for a final RDA and further analyses.

We determined the extent and type of autocorrelation in the independent dataset, as an accurate evaluation of transfer function performance is only possible in the case of spatial independence between sample sites in the dataset (Telford and Birks, 2009). In the case of spatial autocorrelation, the transfer function performance can be overestimated, which could result in inappropriate choices of transfer function models and parameters, possibly including environmental variables with no ecological relevance as “reconstructible”. Autocorrelation could occur due to spatial interdependence of proxy variables, which would render the main requirement of spatial independence for adequate transfer functions violated. It could also be due to an environmental similarity of nearby sample sites. A spatial autocorrelation due to environmentally similar sites, as found in datasets from the North Pacific as well as the North Atlantic and Arctic (Hohmann et al., 2020), may meet the requirement of transfer functions. For autocorrelation tests and further analyses, we additionally used R packages fields (Nychka et al., 2017), palaeoSig (Telford, 2015) and rioja (Juggins, 2017), implementing the routines developed in Telford and Birks (2009).

For the transfer functions remaining in the independent calibration dataset, we did calibration and cross-validation based on the leave-one-out (LOO) (cf. Efron and Gong, 1983) and h-block technique (Burman et al., 1994; Telford and Birks, 2005, 2009) for three distinctly different transfer function approaches: the maximum likelihood (ML), the weighted average partial least square (WA-PLS) and the modern analogue technique (MAT) (for further details about the applied approaches see Hohmann et al., 2020). For calibration and cross-validation tests we additionally used R packages bioindic (Guiot and Gally, 2014), gstat (Pebesma, 2004) and sp (Pebesma and Bivand, 2005) and followed the R code developed by Trachsel and Telford (2016) for h-block cross-validation and the estimation of h.

We reconstructed the three independent variables (SSTsummer, SSSsummer and PPspring) within four sediment sequences collected along a north–south transect through the Baffin Bay area. To determine changes in the environment visible through assemblage shifts downcore, we performed principal component analysis (PCA) on percentages of downcore dinocyst assemblages.

We tested the ability of the calibrated transfer functions to interpret the assemblage changes recorded in the sediments by applying the Telford and Birks method (2011). To estimate the significance of the reconstructions, this method evaluates if reconstructions explain more of the variance in the fossil data than most reconstructions derived from transfer functions trained on random environmental data. It also determines the best reconstruction when reconstructing several environmental variables and if there is sufficient information in the proxy data to support multiple independent reconstructions. We use redundancy analysis (RDA) to estimate the proportion of variance in the fossil data explained by a single reconstruction. Then, using the same taxa data from the calibration dataset, reconstructions are inferred from transfer functions trained on random environmental variables drawn from a uniform distribution. The proportion of the variance explained by these random reconstructions is then estimated. If the tested variable explains more of the taxa variance than 95 % of the random reconstructions, the reconstruction is deemed statistically significant (cf. Telford and Birks, 2011). Reconstructions that fail this test should be interpreted with caution.

Figure 3 RDA scores. (a, b) RDA ordination diagrams of dinocyst taxa (see Table 1 for abbreviations) vs. environmental variables (arrows) showing the correlation with the axes (direction of arrows) and the strength of the dinocyst vs. environmental variable relationship (length of arrows). Species in the centre region are without abbreviations. (c, d) Sample ordination for RDA constrained by independent environmental variables for axes 1, 2 and 3.

Reconstructions and significance tests within the four cores were performed by applying the calibrated local dataset and the regional Northern Hemisphere dataset, aiming to assess the advantages and disadvantages of using the local dataset.

3 Results

3.1 Redundancy analysis

The first step of the variable selection is based on RDA carried out for single variables. It revealed that, when considered individually, each of the six environmental variables appeared to explain a significant (p≤0.01) amount of variation in the assemblage composition in the dataset (Table 4). Together, they explain 32 % of the total inertia in the taxa dataset. We then evaluated the independence of the tested environmental variables through a selection, excluding variables with VIF >2. This successive forward selection uncovered three separate dimensions in the taxa variance driven by SSTsummer, SSSsummer and PPspring, respectively. The remaining variables were excluded from further analyses due to their collinearity with SSTsummer, SSSsummer or PPspring.

Table 4 Explained inertia (as a single variable) and variance inflation factors of the tested environmental variables in the RDA model. VIF (a) when all variables are considered, VIF (b) after manual forward selection for variables without dependence.

When only the independent environmental variables are used, the constrained part of the RDA explains about 24 % of the taxa variance (Table 5), with environmental variables having significant correlations with RDA axes (Table 6). The first axis explains 66 % of the constrained variance and is positively correlated with SSTsummer, while axes 2 and 3 explain 27 % and 7 % and are positively correlated with SSSsummer and PPspring, respectively (Fig. 3, Table 5). Regarding species–environment relationships, we find that the variance of I. minutum is strongly related to changes in SSTsummer, while the abundance of N. labyrinthus is related to SSSsummer. The cyst of P. dalei shows correlations with SSTsummer and PPspring. Phototrophic cysts are generally associated with higher summer temperatures and heterotrophic cysts with lower temperatures (Fig. 3a, b). Sample–environment relationships (Fig. 3c, d) reveal that all samples contribute to the variance in SSTsummer. The variance for PPspring is high for samples from Labrador Sea, Baffin Bay and Hudson Bay, while SSSsummer variance is highest in samples from Davis Strait and Hudson Bay.

Table 5 RDA results and results of permutation tests for the significance of RDA axes after manual forward selection (p≤0.01).

Table 6 Intraset correlation coefficients (r²) between independent variables and RDA axes and between the variables themselves.

3.2 Spatial autocorrelation

We tested spatial and environmental autocorrelation, considering only the independent variables and following the procedure described by Telford and Birks (2009). Inevitably, the performance (r²) deteriorates with increasing fraction of sites deleted (Fig. S2 in the Supplement). We find that the deterioration of performance is much worse with neighbourhood deletion than with random deletion, indicating the presence of spatial autocorrelation. When deleting environmentally similar sites, we find a pattern of performance loss that is slightly worse for neighbourhood deletion than for environmental deletion. Hence, removing geographically close sites causes a higher decline in performance than when the same fraction of environmentally similar sites is removed, indicating spatial autocorrelation over environmental autocorrelation. Variable performance drops below 0.55 when samples within a 200 km radius are deleted, implying that calibrated transfer functions need to be interpreted with caution.

3.3 Cross-validation

For transfer function development, LOO cross-validation results indicate the best performance for MAT, with WA-PLS RMSEP (root mean square error of prediction) values on average around 40 % higher and ML RMSEP values on average by around 100 % higher, depending upon the variable (Table 7). The differences in LOO RMSEPs between the calibration and the verification subsets for MAT show differences of 25 % for SSSsummer and 8 % for PPspring, indicating an inadequate spatial coverage for some areas within the regional subset for these parameters. For SSTsummer reconstructions the spatial coverage in the dataset seems to be sufficient. All three environmental parameters show similar performance within each reconstruction method. The h-block cross-validation shows MAT RMSEP increasing by 60 %, 60 % and 95 % for SSTsummer, SSSsummer and PPspring, respectively. RMSEPs of WA-PLS and ML remain largely unaffected by site deletion within the applied radius h for SSTsummer (1 %) and SSSsummer (7 %) (Table 7). For PPspring RMSEPs increase by about 25 %.

An inspection of the distribution of residuals from cross-validation tests of the variables along their environmental gradient reveals uneven residual distributions for all independent variables. (Fig. S4 in the Supplement). For all transfer function methods, the spatial structure of the residuals is complex and large residual values emerge regionally. For ML large over- and underestimations suggest that the method is inaccurate for reconstructions. For MAT and WA-PLS, we observe a systematic underestimation in the Labrador Sea, especially for SSTsummer and PPspring and an overestimation in parts of Baffin Bay and Hudson Bay with especially large overestimations for PPspring in Hudson Bay. However, residuals for SSSsummer reconstructions regarding MAT and WA-PLS are quite small, suggesting solid transfer functions.

Table 7 Cross-validation performance. Performance for the local dataset applying LOO cross-validation and h-block cross-validation for MAT (c – calibration dataset, v – verification dataset) with the number of analogues used, for WA-PLS with the number of components used and for ML. The number of analogues and components that will perform best during reconstructions has been determined by calculating the RMSE of prediction as a function of the number of analogues used and a randomisation t test for testing the significance of cross-validated components used after van der Voet (1994), respectively; h values have been estimated by fitting a spherical variogram to detrended residuals of a WA model (Fig. S3 in the Supplement) following Telford and Birks (2009) as well as Trachsel and Telford (2016).

3.4 Transfer function significance in sediment core reconstructions

Having established the transfer functions, we applied them to the four selected sediment cores covering the north–south transect. First, for each record, we identified the interval characterised by the largest changes in assemblage composition, which we assumed to be related to the most important environmental change. PC analysis of downcore dinocyst assemblages revealed one to two assemblage changes within each of the four cores (Fig. 4). The first component (PC1) of the PCA identified three assemblage zones in the northernmost record HU91-039-008P, which account for 74.5 % of the total variance, while PC2 accounted for 17.0 %. In southern Melville Bay (GeoB19927-3) and Disko Bugt (MSM343310) PC1 identified two assemblage shifts, accounting for 76.6 % and 77.7 % of the total variance, while PC2 identified at least three shifts, accounting for 17.0 % and 18.3 %, respectively. Within the southernmost record from the Labrador Sea (HU2008-029-004) PC1 revealed four shifts (69.3 %) and PC2 three (21.3 %).

The reconstructions based on the local calibration dataset for the northern core HU91-039-008P (Fig. 4a) show modest changes, but SSTsummer resonates with the facies shifts. In core GeoB19927-3 (Fig. 4b), assemblage shifts correspond to variations in SSSsummer and PPspring reconstructions based on ML and WA-PLS. In core MSM343310 (Fig. 4c), the facies shift is caught in SSSsummer and SSTsummer reconstructions based on ML and WA-PLS. In core HU2008-029-004 (Fig. 4d), we identify shifts based on WA-PLS and ML reconstructions.

The different degree of correspondence between the environmental reconstructions inferred by the transfer functions and the principal community change trends in the dinocyst records is echoed in the significance tests made after the approach of Telford and Birks (2011) for downcore reconstructions. Here a significant reconstruction should explain more of the variance in the fossil dinocyst assemblage than a random transfer function calibrated using the same assemblages but random environmental data. These tests reveal that different transfer function methods and environmental variables produce significant reconstructions for each core. Significance test results of transfer function applications with the local dataset and the regional Northern Hemisphere dataset reveal the same pattern (Table 8).

When the Northern Hemisphere dataset is used for reconstructions of the four downcore assemblages (Fig. 4), we find that 17 % of all the tested transfer functions produce significant estimates (SSSsummer – 25 %, SSTsummer – 17 %, PPspring – 8 %) (Fig. 5, Table 8).

However, when applying the local calibration dataset, most transfer functions produce significant reconstructions for SSSsummer (75 %) and SSTsummer (50 %) and some for PPspring (8 %) (Fig. 5, Table 8). For the northern core HU91-039-008P, summer salinity estimates are significant with all three transfer function methods (Fig. 5a). In southern Melville Bay core GeoB19927-3, the WA-PLS transfer function has been tested as significant for PPspring, as well as SSSsummer and SSTsummer (Fig. 5b). For records in Disko Bugt (MSM343310), both MAT and WA-PLS produce significant estimates for SSSsummer, while WA-PLS shows additional significance for SSTsummer (Fig. 5c). Further south, the record from the Labrador Sea (HU2008-029-004) yields significant estimates for SSSsummer and SSTsummer based on MAT and ML (Fig. 5d).

Figure 4 Downcore reconstructions based on the local dataset (n=421) and the Northern Hemisphere dataset (n=1968) using three transfer function techniques. MAT, WA-PLS and ML are represented by blue, green and purple curves, respectively. The lighter-blue curves correspond to maximum and minimum possible values calculated from a set of five modern analogues. The ecostratigraphic zones are determined from PC analyses (dashed horizontal lines). The “|” signs (red) on the upper horizontal axis correspond to modern conditions. (a) HU91-039-008P, (b) GeoB19927-3, (c) MSM343310 and (d) HU2008-029-004.

Table 8 Significance for downcore reconstructions after the approach of Telford and Birks (2011) using the local dataset and the regional Northern Hemisphere dataset. Significant reconstructions are marked by a check mark and the non-significant ones by a cross.

Figure 5 Results of the significance tests for downcore reconstructions. (a) HU91-039-008P, (b) GeoB19927-3, (c) MSM343310 and (d) HU2008-029-004. For all significance tests, 999 random environmental variables were generated to produce the null distribution (grey bars). Black lines mark the variances explained by observed variables. Dotted red lines mark the 95th percentile of the random distribution (p≤0.05) above which reconstructions are deemed as significant. The proportion of the variance explained by the first axis of a principal component analysis (PCA) is also recorded (dotted black line), as this represents the maximum proportion of the variance in the fossil data that a reconstruction could possibly explain. For significance test analysis we applied the number of analogues (MAT) and components (WA-PLS) resulting in the highest significance with respect to all three variables.

4 Discussion

4.1 Methodology

The ecological models and transfer functions derived in this study rely on a few assumptions. Foremost is that the composition of dinocyst assemblages was primarily affected by the selected environmental parameters and that the set of parameters is comprehensive and includes all the relevant parameters. However, it is possible that a sedimentary assemblage reflects other processes than the considered parameters. For example, non-oceanographic factors such as those related to the dinoflagellate life cycle (Bravo and Figueroa, 2014) or post-depositional organic matter degradation (Zonneveld et al., 2019) could affect the dinocyst assemblage composition in the calibration dataset. The production of cysts is a part of the dinoflagellate's sexual reproduction. It protects the diploid cell, which may be dormant for a variable length of time. After dormancy the cell emerges from the cyst and re-establishes its metabolic activities within the pelagic ecosystem (Taylor and Pollingher, 1987). During encystment many cysts settle on the seafloor, requiring areas where the seafloor is not too deep to ensure the possibility of excystment and ascendance of the cells into surface waters. Hohmann et al. (2020) showed that water depth and distance from shore explain a significant amount of compositional variance in dinocyst assemblages (p<0.01). However, these variables were shown to describe a separate dimension of variability, orthogonal to the tested oceanographic and primary productivity factors. Also, when applied to a sedimentary record, factors such as the distance from shore and water depth should be mostly constant for the downcore record. The effect of organic matter degradation on the dinocyst assemblage composition can be detected independently of the ecological model because the ratio between cysts of autotrophic and heterotrophic taxa should systematically decrease with an increase in selective organic matter degradation, caused by their different wall polymer composition (Zonneveld et al., 1997). In the Northern Hemisphere calibration dataset, there is no evidence of a pervasive effect of degradation. In neither the Northern Hemisphere (Hohmann et al., 2020) nor the Baffin Bay local dataset (this study) has a systematic relationship between the residuals and the ratio of heterotrophic and autotrophic taxa been observed.

Since selective degradation plays at most a subordinate role in biasing dinocyst assemblages in this study, and life-cycle factors explain an independent dimension of the taxa variance, we consider the dinocyst assemblage compositions in the calibration datasets as being mostly driven by environmental factors. However, most of the taxa variance in the local calibration dataset cannot be explained even when all tested environmental factors (including water depth and distance from shore) are combined. First, the combination of taxa to achieve higher consistency may hide a finer ecologically relevant taxonomic level. Since studies revealed the existence of cryptic species in modern dinoflagellates (Montresor et al., 2003; Parkinson et al., 2016; Wang et al., 2019) and since some taxa are grouped voluntarily during analysis, finer relevant taxonomic levels are very likely. It is, however, currently unavoidable because routine analyses (de Vernal et al., 2020) do not allow a refinement in regard to taxonomic levels.

Further, cyst-producing dinoflagellates lose their mobility at cyst formation, from which time onwards they are subject to ocean current transport until they reach the ocean floor. Although we excluded sample side areas with strong oceanographic currents, small-size dinocyst particles might still be subject to transport by moderate currents when not incorporated in marine snow or faecal pellets (e.g. Turner, 2002, 2015) or bottom water currents and sediment flows (Dale, 1992; Dale and Dale, 1992). This may result in a lateral relocation and sedimentation from a cyst's formation environment and thus in an additional nuisance in the calibration dataset (Nooteboom et al., 2019). Biases in the reconstructions may also result from advected fossil dinocysts in the past when ocean currents were different. While this might not be eminent for stratigraphies in Baffin Bay during the Holocene, in general cores from such regions need to be treated with caution. Additional noise may also emerge from modern assemblages in the calibration dataset, which represent average fluxes over different time intervals. Assemblages may be affected by anthropogenic effects, such as local pollution or global warming, to a different degree. Assemblages from sites with low sediment accumulation may represent averages over millennia and contain a mixture of microfossils deposited during distinct oceanographic events of the Late Holocene (e.g. Wigley et al., 1981; Bradley and Jonest, 1993) or be affected by cysts deposited during conditions induced by modern anthropogenic global change.

Alternatively, it is possible that the tested parameters are not the most determinant ones or that the assemblage composition is to a certain degree driven by biotic interactions rather than by environmental factors. It is also possible that environmental factors, not yet considered, play an important role. These could include turbulence, micronutrients, stratification or any other parameters that might be difficult to quantify for data treatment purposes.

Regardless of the limitations that can be due to missing or unidentified driving parameters, we find three mutually independent ecological gradients driving dinocyst assemblages in the local calibration dataset. This confirms previous findings that phytoplankton (Lopes et al., 2010) and hence dinocyst assemblages reflect multiple environmental drivers (e.g. Radi and de Vernal, 2008; Ribeiro and Amorim, 2008; Hohmann et al., 2020).

Figure 6 Dinocyst assemblage changes according to the first two PC axes. Coloured dots represent PC transitions pointing to regime shifts around 1.3 ka (between 1.0 and 1.5 ka), 4.0 ka (between 3.1 and 4.3 ka) and 6.7 ka (shaded orange).

The existence of multiple environmental gradients affecting the dinocyst assemblage composition implies that it should be possible to extract information from the same fossil dinocyst assemblage on past variability for more than one environmental driver. Whilst this is correct theoretically, one must consider the (small) amount of variance constrained by the tested parameters against the (large) portion of the variance that remains unconstrained. For example, primary productivity, which represents the third gradient after SST and SSS, explains a maximum of 6 % of the constraint variability in the RDA (Tables 5 and 6). This parameter does not contribute to the first two dimensions of the RDA (see also Radi and de Vernal, 2008; de Vernal et al., 2020), and although the extracted relationship appears statistically significant, it may well be overwhelmed by factors that are responsible for the unconstrained portion of the variance in assemblage composition.

Additionally, we note that the presence of multiple environmental gradients in the local calibration dataset is associated with spatial structuring of the environment (Fig. S4 in the Supplement). This may result in the fact that some environmental drivers deemed significant in the RDA are only significant because they drive the taxa composition in a specific area with a strong gradient of the given factor. This is in line with the observed spatial autocorrelation of the tested variables (Fig. S2 in the Supplement) and explains the performance decrease in MAT in the h-block cross-validation (Table 7) as well as the strong spatial structure observed in the residuals of all methods (Fig. S4 in the Supplement). These observations are coherent with findings in Hohmann et al. (2020), where the observed multiple and mutually independent environmental factors driving the dinocyst assemblage composition in two larger regional datasets may result from the agglomeration of information from regions with unlike environmental (primary) drivers and endemic taxon–environment relationships. Hence, our results suggest localised driving mechanisms for dinocyst assemblage compositions in very specific areas. This implies that generating locally calibrated transfer functions will come at the cost of loss of generalisation and limited applicability to past oceanographic states.

Being based upon differing basic principles, the applied transfer function approaches seem to be affected differently by the detected local dataset's spatial structuring and the region-dependent primary drivers (Kucera et al., 2005). ML transfer functions appear to be unable to disclose and depict the taxon–environment relationship when multiple gradients affect the taxa assemblages (highest RMSEPs, lowest r², fewest significant transfer functions). However, ML and WA-PLS performance are less affected by the effect of spatial autocorrelation. During h-block cross-validation, their performance remains similar, while the performance of MAT transfer functions plummets, returning higher RMSEP and lower r² than ML and WA-PLS (Table 7). On the other hand, the performance of MAT based on LOO cross-validation is more optimistic than for ML and WA-PLS. While MAT finds the local structure in the assemblage data, ML and WA-PLS, which are based on a unimodal species–environment response model, find the general component of the variation in the assemblage data that are correlated with a specific environmental forcing. While all approaches have their advantages, ML and WA-PLS are more robust with regard to the detected spatial structure in the data than MAT (Telford and Birks, 2005). Cross-validation of the local dataset applying a calibration and verification dataset results in a 25 % lower RMSEP_v for the verification dataset compared to the calibration dataset for SSSsummer, while RMSEP_v values are <10 % lower for SSTsummer (2 %) and PPspring (8 %). This indicates that additional sample sites might benefit performance for SSSsummer transfer functions as with an adequate, hence representative, spatial coverage RMSEP should be similar for both datasets (Table 7).

Regardless of the statistical approach, a theoretical high-performance transfer function for one of the independent driving variables for assemblage composition in the local calibration dataset does not necessarily produce a reliable reconstruction applied to fossil downcore dinocyst assemblages. This is because we can never be sure that the main drivers for past assemblages were identical in the modern calibration dataset. To assess if the main past drivers resemble the main modern drivers, we compared the transfer function reconstructions to the major compositional trends in fossil dinocyst assemblages as revealed by multivariate analyses.

Across the four records throughout Baffin Bay, the main pattern of dinocyst assemblage change during the Holocene can be identified as two biofacies shifts occurring around 4.0 and 6.7 ka (Figs. 4 and 6) throughout the whole study area and one regime shift occurring in Baffin Bay and Disco Bay between 1.0 and 1.5 ka. Despite the presence of this clear pattern, most transfer functions based on the local and Northern Hemisphere datasets reconstruct modest changes close to uncertainty for SSTsummer, SSSsummer and PPspring (Fig. 4), likely reflecting the difficulty in attributing the taxa assemblage change to specific factors. Alternatively, one may argue that we have not retained the correct determinant parameters during the manual forward selection. For example, the Holocene records could reflect changes that were due to large shifts in sea ice coverage. Indeed, during the selection process, the analysis suggested that SeaIce does not explain a separate dimension of variance in the taxa assemblages within the local Baffin Bay dataset (Script S1 in the Supplement), as it revealed high collinearity with SSTsummer. This means that regionally SeaIce and SSTsummer explain a similar taxa variance, and while reconstructing one, we reconstruct a mixture of both parameters. It is also possible that the parameterisation of sea ice cover used here is not ideal to describe the relationship with dinocyst assemblages. In the context of Baffin Bay, SSSsummer, which is the parameter that appears to be best reconstructed, is likely a reflection of the regional hydrography closely related to the sea ice cover that matters for productivity. Therefore, we can only conclude that the environmental regime shifts around 1.3, 4.0 and 6.7 ka appear to be driven by changes in local hydrography that were associated with salinity changes, but considering the magnitude of the salinity changes implied by the transfer functions, it is unlikely that the dinocyst composition responded to salinity directly.

In contrast to the results of the local calibration, only a few of the reconstructions were identified as significant when applying the regional Northern Hemisphere dataset (SSTsummer = 16 %, SSSsummer = 25 %, PPspring = 8 %) (Table 8, Fig. 5). This suggests that transfer functions based on broad regional or even global calibration datasets are affected by nuisance variables, which makes it difficult to identify the driving mechanisms of dinocyst assemblages when the amplitude of environmental change is relatively small, as is the case during the Middle and Late Holocene. For longer time intervals with large amplitude changes, global calibrations may work better and would be even necessary as larger datasets contain a wider range of analogues. Hence, we suggest a thorough evaluation of the calibration dataset and the significance of fossil applications, ensuring that the driving mechanisms of dinocyst assemblages are caught through time and space.

Figure 7 Time series for significant SSSsummer reconstructions (Table 8) (local MAT and WA-PLS) from the four studied cores compared to different environmental records from the Baffin Bay region. (a) The Greenland Ice Sheet area (Larsen et al., 2015), (b) Laurentide Ice Sheet area (Dyke et al., 2003) (plots after Briner et al., 2016), (c) Holocene Greenland temperatures based on δ¹⁸O from the GISP (Greenland Ice Sheet Project) ice core (Alley et al., 2010), (d) sea ice cover in southern Melville Bay based on biomarkers (Saini et al., 2020) (MIZ – marginal ice zone), (e) sea ice cover in Upernavik based on dinocyst assemblages (three-point running mean) (Caron et al., 2019), (f) sea ice cover in the southern Baffin Bay based on dinocyst assemblages (five-point running mean) (Gibb et al., 2015; Jennings et al., 2014), (g) primary productivity in Upernavik based on dinocyst assemblages (three-point running mean) (Caron et al., 2019), (h) SSSsummer from core HU91-039-008P, (i) SSSsummer from core GeoB19927-3, (j) SSSsummer from core MSM343310 and (k) SSSsummer from core HU2008-029-004. Regime shifts (Fig. 6) are shaded in grey.

4.2 Salinity reconstructions and regime shifts in the context of the existing regional palaeoceanographic records

Results imply that in Baffin Bay most locally calibrated MAT and WA-PLS transfer functions produce robust Holocene reconstructions for summer salinity, being the main driver of dinocyst assemblage compositions while representing the regional hydrography closely related to sea ice cover and sea surface temperature driving the salinity changes. This poses the opportunity to interpret SSSsummer reconstructions within the context of the regional palaeoceanography.

Until 7.6 ka Baffin Bay was still a region of extended sea ice and low primary productivity (e.g. Ouellet-Bernier et al., 2014; Jennings et al., 2014; Gibb et al., 2015; Caron et al., 2019; Limoges et al., 2020; Saini et al., 2020) (Fig. 7d–g), while the Laurentide and the Greenland ice sheets continuously retreated (Fig. 7a and b) (Dyke et al., 2003; Larsen et al., 2015).

With the beginning of the Holocene Climate Optimum, sea ice cover in Baffin Bay retreated, and it became a region of marginal ice zones until 6.3 ka. During this time, the first observed regime shift occurred between 7.6 and 6.3 ka (Fig. 6), which is in accordance with the shift towards mild postglacial conditions in Baffin Bay proposed by Gibb et al. (2015) and Caron et al. (2019). In the north-western Labrador Sea (core HU2008-029-004) however, this postglacial transition already occurred around 11.9 ka (Gibb et al., 2015), implying that the observed regime shift also visible in both PCA axes of HU2008-029-004 cannot – like in Baffin Bay – be ascribed to the postglacial transition.

After the transition phase from 6.3 ka on, salinity reconstructions propose a continuous increase in the north-western Labrador Sea and a continuous decrease in salinity in southern Melville Bay (Fig. 7i and k). For the Labrador Sea this is in accordance with findings in Gibb et al. (2015), who ascribe the increase to a reduction in stratification in the upper waters caused by the reduction in meltwater input from the Laurentide ice sheet and the increased advection of North Atlantic water (de Vernal et al., 2001; Solignac et al., 2004). Hence, the observed regime shift in the north-western Labrador Sea between 7.6 and 6.3 ka is in line with the suggested regional shift towards modern water circulation south of Davis Strait (Gibb et al., 2015).

For southern Melville Bay, other records also propose a slow salinity decrease between 5.5 and 1.4 ka (Gibb et al., 2015; Caron et al., 2019), which might be attributed to the still-active process of postglacial meltwater discharge north of Davis Strait. This can be supported by a continuous decrease in sea ice cover in this region (Fig. 7d–f; Gibb et al., 2015; Caron et al., 2019; Saini et al., 2020) and the increase in primary productivity after the first observed regime shift (e.g. Caron et al., 2019; Limoges et al., 2020). At the south-western Greenland margin dinoflagellate blooms in summer–autumn have been associated with freshwater inputs due to meltwater runoff from the GIS, which triggers upwelling of nutrient-rich waters and high productivity (Juul-Pedersen et al., 2015; Krawczyk et al., 2015, 2018). Along most of the West Greenland margin, sea ice melt and meltwater runoff from the GIS result in seasonal stratification of surface waters, which is amplified in summer by solar heat, leading to relatively mild conditions (Juul-Pedersen et al., 2015; Tremblay et al., 2015) and also enhancing primary productivity. Pollen data indicate a high influx of long-distance pollen over north-western Greenland during this time, indicating a milder climate (Funder, 1978; Fredskild, 1985; Fredskild and Andrews, 1985).

In the northernmost Baffin Bay, our reconstructions show no general trend but a continuously oscillating salinity (Fig. 7h), which is supported by reconstructions in Levac et al. (2001). Like in Levac et al. (2001), our reconstructed values stay below 32, indicating lower values than today's modern instrumental values.

However, while a general trend towards lower salinity can be observed north of Davis Strait and south of the northernmost Baffin Bay after 5.5 ka, other studies (Gibb et al., 2015; Caron et al., 2019) also find a salinity peak around 3.7 ka, interrupting the general trend, which is not visible in our reconstruction, although, this salinity peak corresponds to the second regime shift indicated by the PCA of all four studied cores (Fig. 6). This change in dinocyst assemblage compositions is also mentioned in Caron et al. (2019) for southern Melville Bay, although other dinocyst records from this region did not observe a shift (Ouellet-Bernier et al., 2014; Gibb et al., 2015). In northernmost Baffin Bay, our reconstructions point to a period of slightly increased salinity during the regime shift (Fig. 7h). For the Baffin Bay area this shift may mark the final adaptation of assemblage compositions to full interglacial conditions, supported by still continuously decreasing sea ice cover (Fig. 7) and other palaeoclimate studies along the south-western Greenland margin, suggesting that optimum conditions were attained between 6.0 and 3.5 ka (Moros et al., 2006; Lloyd et al., 2007; Seidenkrantz et al., 2007; Ren et al., 2009; Andresen et al., 2010; Perner et al., 2013).

In northernmost Baffin Bay and north-western Labrador Sea however, the changes in assemblage compositions are not stable, changing again within the supposed period of the regime shift (Fig. 6). For northern Baffin Bay Levac et al. (2001) find a general deterioration of sea surface conditions after 3.6 ka, which was accompanied by extensive sea ice cover. These findings are supported by the second regime shift, implying changes in the dinocyst assemblage composition.

The Late Holocene is locally characterised by the neoglacial cooling trend (Weidick et al., 2012), which along the West Greenland coast is associated with subsurface cooling and glacier advances (e.g. Perner et al., 2013; Briner et al., 2016; Schweinsberg et al., 2017) over the last 4.0 ka (e.g. Knudsen et al., 2008; Seidenkrantz et al., 2008; Krawczyk et al., 2010; Andresen et al., 2010). This cooling trend is marked by climate oscillations, including the short warming phases of the Roman Warm Period (2250–1600 years BP) (Lamb, 2002) and the Medieval Climatic Optimum (1100–700 years BP) (Lamb, 1965; Mann, 2002). In Baffin Bay records, the trend's impact becomes visible from around 2.2 ka (Fig. 7). Regional cooling accompanied by decreasing summer solar insolation (Berger and Loutre, 1991), fostering regrowth of the GIS (Briner et al., 2010; Weidick and Bennike, 2007; Young et al., 2011), as well as cooling over the GIS (Dahl-Jensen et al., 1998) led to more extensive sea ice, which was accompanied by increasing sea surface salinity and a decline in primary productivity. Around 1.3 ka the southern Melville Bay core (GeoB19927-3) and the Disco Bay core (MSM343330) suggest a regime shift, supporting the assumption of a transition towards colder climate. This is also in accordance with proxy records from eastern Baffin Bay, suggesting glacier advances (Ouellet-Bernier et al., 2014; Gajewski, 2015) and decreased terrigenous detrital supply (Caron et al., 2020) in this time interval.

Overall, the reconstructions for summer salinity and the observed regime shifts fit into and support the current picture of West Greenland's palaeoceanography. Comparing salinity reconstructions for southern Melville Bay (GeoB19927-3) with reconstructions of sea ice cover from the same core derived from the biomarker proxy P_DIP₂₅ (Belt and Müller, 2013) for late-spring sea ice cover (Saini et al., 2020), similar changes can be observed. Increasing sea ice cover correlates with increasing salinity and vice versa. This additionally supports the significance of SSSsummer WA-PLS reconstructions in the Melville Bay core and supports the assumption that salinity and sea ice cover are tightly coupled environmental drivers in this area.

5 Conclusions

Understanding past natural variability in ecosystem drivers in Arctic regions requires reliable reconstructions. We calibrated local transfer functions from Baffin Bay samples for three independent parameters driving dinocyst assemblages in this area. We estimated how well these local functions were able to estimate past environmental parameters in the Baffin Bay region during the Holocene, comparing them to estimates derived from regional transfer functions from Northern Hemisphere samples and afterwards comparing the significant local transfer functions with the existing regional palaeoceanographic records.

Statistical analyses of surface sediment samples and transfer function application on downcore samples from a Holocene transect along the Baffin Bay indicate that salinity is the main driver of present-day and Holocene dinocyst assemblage composition in the region, with temperature being an important but minor driver.

Analyses of four Holocene records in a north–south transect along the Baffin Bay show three regime shifts that occurred throughout the region. The main environmental parameter to which the dinoflagellate communities seem to have reacted during regime shifts is salinity, but further minor drivers may have differed regionally. The reconstructed salinity records and regime shifts support the current state of knowledge of West Greenland's Holocene palaeoceanography.

We show that the relationship between dinocyst assemblage composition and multiple environmental drivers can be modelled by MAT and WA-PLS transfer functions. In our study, the identification of main driving variables and regime shifts was achieved through the application of a local calibration dataset, while through the application of the Northern Hemisphere dataset, we did not succeed in extracting this information from the taxa variability. While our results reinforce the assumption of counteracting nuisance in spatially large calibration datasets, suggesting the application of local datasets, regions with spatially or temporarily broader environmental gradients might supposedly benefit from larger calibration datasets regardless of a certain nuisance.

We suggest a thorough evaluation of transfer function performance and significance for downcore applications to disclose the specific drivers for present and fossil dinocyst assemblages in a studied core location. To this end, a local calibration dataset, which reduces several nuisances, is presumably beneficial when the amplitude of parameter change is relatively small, like during the Middle and Late Holocene.