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 2List of dinocyst taxa in the local calibration dataset, including abbreviations and highest relative abundances (minimum is always zero). Asterisks indicate heterotrophic taxa.

Download Print Version | Download XLSX

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

Download Print Version | Download XLSX

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

Download

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 3RDA 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.

Download

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.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 4Explained 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.

Download Print Version | Download XLSX

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 5RDA results and results of permutation tests for the significance of RDA axes after manual forward selection (p≤0.01).

Download Print Version | Download XLSX

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

Download Print Version | Download XLSX

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 7Cross-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).

Download Print Version | Download XLSX

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 4Downcore 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.

Download

Table 8Significance 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.

Download Print Version | Download XLSX

Figure 5Results 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.

Download

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 6Dinocyst 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).

Download

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 7Time 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.

Download

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.