Deglacial increase of seasonal temperature variability in the tropical ocean
Study siteThe Cariaco Basin, located on the continental shelf off Venezuela, is a large (about 160 km long and about 65 km wide) depression, composed of two approximately 1,400-m-deep sub-basins. It is partially isolated from the Caribbean Sea by a series of sills with depths of less than 150 m (ref. 47). This limits renewal of deep water in the basin and, paired with the high oxygen demand resulting from intense surface primary productivity, leads to anoxic waters below a depth of about 275 m at present47,48.The marked seasonality in the Cariaco Basin, combined with anoxic bottom waters that effectively prevent bioturbation, results in the accumulation of annually laminated (varved) sediments. As sediments are varved for the last deglaciation and the Holocene, and because of the sensitivity of the area to climate change, they are considered to be one of the most valuable high-resolution marine climate archives and have been successfully used to study climate variability in the tropics3,11,16,17,18. Varve thickness is about 1 mm or more during the YD–Holocene transition18.Core and age modelCore MD03-2621 was retrieved during IMAGES cruise XI (PICASSO) aboard R/V Marion Dufresne in 2003 (Laj and Shipboard Party 2004). Cariaco cores have been collected under the regulations of the Ocean Drilling Program and the IMAGES coring programme. In this study, data from depths between 480 and 540 cm below the seafloor are presented, encompassing the YD–Holocene transition. A detailed age model for core MD03-2621 was established by Deplazes et al.11 and is based on the cross-correlation of total reflectance to dated colour records from the Cariaco Basin49,50. For the studied interval, the original age model is based on a floating varve chronology anchored to tree ring data by matching 14C data49. The age model for core MD03-2621 was further fine-tuned by correlation of reflectance data to the NGRIP ice core δ18O record on the GICC05 age scale11. The transition from the YD to the Holocene is characterized by a decrease in the sedimentation rate from 1.4 to 0.5 mm year−1.To account for possible depth offsets during storage and subsampling, we matched sediment colour data expressed as greyscale (GS) to the reflectance data from Deplazes et al.11 with the software QAnalySeries51. To enable comparison with our record, ages in Lea et al.3 were corrected for the age difference between the sediment-colour-based midpoint of the YD–Holocene transition in their record (11.56 kyr b2k) and in data from Deplazes et al.11 (11.673 kyr b2k). The start and end of the change in reflectance were determined by the RAMPFIT software52.Sample preparationSamples for MSI of molecular proxies were prepared as described in Alfken et al.53: the original core was subsampled by LL channels, from which X-ray pictures (Hewlett-Packard Faxitron 43855A X-ray cabinet) and high-resolution digital images (smart-CIS 1600 Line Scanner) were obtained. The LL channels were then cut into 5-cm pieces, which were subsequently freeze-dried, embedded in a gelatin:carboxymethyl cellulose (4%:1%) mixture and thin-sectioned on a Microm HM 505 E cryomicrotome. From each piece, one 60-µm-thick and one 100-µm-thick, longitudinal slice (spanning the whole 5 cm piece) were prepared and affixed to indium-tin-oxide-coated glass slides (Bruker Daltonik, Bremen, Germany) for MSI and elemental mapping, respectively. Slices for MSI were further amended with a fullerite matrix54.For all slices, a high-resolution picture was taken on a M4 Tornado micro-X-ray fluorescence spectroscopy system (Bruker Nano Analytics). This picture was used as a reference to set up elemental mapping and MSI analysis, and also for the 2D comparison of elemental and proxy data to sediment colour. Sediment colour is expressed as GS value. To account for differences between single slices, ΔGS was calculated as the difference between a value and the median GS of each individual slice. Very low GS values corresponding to areas devoid of sediment, identified by a black background, were excluded from analysis.Elemental mappingElemental mapping of 100-µm-thick slices was performed on a M4 Tornado micro-X-ray fluorescence spectroscopy system (Bruker Nano Analytics) equipped with a micro-focused Rh source (50 kV, 600 µA) with a polycapillary optic. Measurements were conducted under vacuum, with a resolution of 50 µm, two scans per spot and a scan time of 5 ms per scan. Data were initially processed and visualized with M4 Tornado Software version 1.3. XY matrices of relevant elements and sediment colour were imported into MATLAB (R2016b) for further processing. To assess the correspondence between sediment colour and elemental composition, for each 5-cm piece, signal intensities of Ca, Fe, Ti and Si in single spots were binned according to ΔGS and average intensities were calculated for each bin (Extended Data Fig. 5). The bin size was 5 GS units.Molecular proxy analysis by MSIMSI was carried out on a 7T solariX XR Fourier transform ion cyclotron resonance mass spectrometer coupled to a matrix-assisted laser desorption/ionization source equipped with a Smartbeam II laser (Bruker Daltonik, Bremen, Germany). Analyses were performed in positive ionization mode selecting for a continuous accumulation of selective ions window of m/z 554 ± 12. Spectra were acquired with 25% data reduction to limit data size. Spatial resolution was obtained by rastering the ionizing laser across the sample in a defined rectangular area at a 100-µm spot distance. Considering laminae thickness in the millimetre range18, such raster resolution is suited for seasonally resolved SST reconstruction. Settings for laser power, frequency and number of shots were adjusted for optimal signal intensities before each measurement; typical values were 250 shots with 200 Hz frequency and 60% laser power. External mass calibration was performed in electrospray ionization mode with sodium trifluoroacetate (Sigma-Aldrich). Each spectrum was also calibrated after data acquisition by an internal lock mass calibration using the Na+ adduct of pyropheophorbide a (m/z 557.2523), a chlorophyll a derivative generally present in relatively young marine sediments. Around 20,000 individual spots were thereby obtained for every 5-cm slice, each spot containing information on the abundance of diunsaturated and triunsaturated C37 alkenones needed to calculate the ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) SST proxy.We provide a two-pronged approach to decode SST proxy information: (1) a downcore ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) profile is obtained by pooling alkenone data from coeval horizons, and results in SST reconstructions with annual resolution, and (2) 2D images of alkenone distribution are examined in conjunction with maps of sediment colour and elemental distribution to filter single-spot alkenone data for season of deposition.SST reconstruction with yearly resolutionFor the downcore profile, MSI data were referenced to the X-ray image by the identification of three teaching points per 5-cm piece. Afterwards, the X-ray image was corrected for tilting of laminae in the LL channels. This was done by identification of single laminae in the X-ray image and selection of a minimum of four tie points per lamina. A detailed description can be found in Alfken et al.9. After applying the corresponding age model, downcore profiles were established with 1-year resolution: the intensity of the two alkenone species relevant to the ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) proxy (C37:2 and C37:3) were recorded for each individual laser spot and filtered for a signal-to-noise threshold of 3. Only spots in which both compounds were detected were further considered. Intensity values were then summed over the depth corresponding to 1 year. By pooling proxy data into 1-year horizons, the effect of changing sedimentation rate and, thereby, changing downcore resolution is minimized. If at least ten spots presenting both compounds were available for a single horizon, data quality criteria were satisfied54 and a ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) value was calculated as defined by Prahl and Wakeham22:$${{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}=frac{{text{C}}_{37:2}}{{text{C}}_{37:2}+{text{C}}_{37:3}}$$
(1)
To apply the gas chromatography (GC)-based calibrations for the ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) proxy, MSI-based data were converted to GC equivalents. Therefore, after MSI, sediment slices were extracted for conventional proxy analysis. Sediment was scraped off the slide and extracted following a modified Bligh and Dyer procedure55,56. Extracts were evaporated under a stream of nitrogen, re-dissolved in n-hexane and analysed on a Thermo Finnigan Trace GC-FID equipped with a Restek Rxi-5ms capillary column (30 m × 0.25 mm ID). For each 5-cm piece, a ratio between the ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) values obtained by GC flame ionization detector analysis and MSI of the whole piece was calculated. The average ratio of all pieces for which GC-based values could be obtained was 1.194, with a standard deviation of 0.021.$${{rm{U}}}_{37,{rm{G}}{rm{C}}-{rm{F}}{rm{I}}{rm{D}}}^{{{rm{K}}}^{{prime} }}=1.194times {{rm{U}}}_{37,{rm{M}}{rm{S}}{rm{I}}}^{{{rm{K}}}^{{prime} }}$$
(2)
This ratio was used to calculate GC-equivalent ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) values, which were then translated into SST using the BAYSPLINE calibration57. The average standard error of the BAYSPLINE model is 0.049 ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) units (corresponding to 1.4 °C) for samples with SST below 23.4 °C, but increases at higher values (to up to 4.4 °C)57. This is explained by the fact that sensitivity of the ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) to SST (that is, the slope of the regression) declines at higher values. In the current dataset, the 95% confidence interval is, on average, ±3.6 °C. The analytical precision of MSI-based SST reconstructions for the ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}), using at least ten data points, according to Alfken et al.9, is about 0.3 °C. Sources of uncertainty are summarized in Extended Data Fig. 10a.For frequency analysis, a continuous, annually spaced record was constructed by linearly interpolating 49 missing values. The record was subsequently detrended. Spectral analysis was performed with the REDFIT module58 using a Hanning window (oversample 2, segments 2). Continuous wavelet transforms were applied to investigate changes in cyclicity over time, using the Morlet wavelet with code provided by Torrence and Compo59 for MATLAB. All steps, except for the wavelet analysis, were performed with the PAST software60.For the assessment of the interannual variability, the SST record was band-pass-filtered for periods between 2 and 8 years. The record is based on 1-year binned data; seasonality is thereby nullified and the highest frequency to be evaluated (Nyquist frequency) corresponds to a period of 2 years. Variability of this time series was quantified by calculating the standard deviation of the band-pass-filtered ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) signal in 25-year intervals. To account for the potential impact of analytical precision on the observed signal (Methods, section titled ‘The effect of changing sedimentation rate on reconstructed interannual SST variability during the YD–Holocene transition’), the variability experiment from Alfken et al.9 was revisited. A sediment extract had been sprayed on an ITO slide and analysed by MSI. We then randomly selected n spots and obtained a ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) value for the summed intensities of these spots. Precision was calculated as the standard deviation of five replicate experiments for n = 1, 5, 10, 20, 30, 40, 50 and 60. Decreasing analytical variability with increasing number of observations was fitted to a curve (R2 = 0.838) described by the equation$${rm{Analytical}},{rm{variability}}=0.0741times {rm{number}},{rm{of}},{{rm{spots}}}^{-0.558}$$
(3)
On the basis of this equation, analytical variability for each horizon could be calculated on the basis of the number of values included (Extended Data Fig. 10b). The mean variability for each 25-year window was then subtracted from the observed variability in the band-pass-filtered signal and the resulting proxy values were translated to SST following the equation by Müller et al.61. Statistical significance of the change in corrected SST variability after 11.66 kyr b2k was assessed with a t-test.Assessment of SST seasonalityFor the assessment of SST seasonality, alkenone intensities from individual spots were binned according to ΔGS, with a bin size of 1 unit. Spots were then separated into the categories upwelling season and non-upwelling season by identifying the threshold ΔGS value that maximized the difference between average SST in the bins above and below it. Furthermore, this value had to fulfill three conditions: (1) be higher (lighter) than the bins with the highest relative abundance of Ca, Ti and Fe, which is indicative of the dark sediments associated to non-upwelling season, (2) be lower (darker) than the bin with highest relative abundance for Si indicative of light sediment associated to the upwelling season and (3) the number of spots categorized as upwelling and non-upwelling had to account for at least 25% of total spots. If criteria 1 and 2 prevented criteria 3 from being fulfilled, a limit of 15% was set. After separating data into these two categories, data were processed separately as described above for the unfiltered dataset and a downcore temporal resolution of 5 years was applied. Seasonality was calculated as the difference between both records and thus represents the difference between 5-year average SST in the non-upwelling and upwelling seasons.Shift in seasonality was fitted to two different ramps with the RAMPFIT software52. An unconstrained approach and a constrained approach (in which the start and end points of the ramp were restricted to the intervals 11.725–11.8 kyr b2k and 11.6–11.675 kyr b2k) were applied. Negative values were excluded from this fitting. The resulting groups of data were compared by a Mann–Whitney rank test.SST seasonality in the modern Cariaco Basin was calculated for the years 1980 to 2020 based on the HadISST dataset62 by dividing monthly data from each year into two groups and searching for the largest difference between the average temperatures of both groups. Each group had to include at least three consecutive months. In 36 out of 41 years, the warm season was defined from May to November or from July to November.Decadal-scale to centennial-scale SST changes during the YD–Holocene transition and in the early HoloceneAnnually reconstructed SST (average SST = 24.3 °C) remains relatively stable during the YD–Holocene transition. At around 11.4 kyr b2k, a warming trend is observed. Averaging all data before 11.39 kyr and after 11.37 kyr results in a warming from 23.9 ± 1.6 °C to 25.5 ± 1.4 °C. Trends identified by MSI are consistent with conventional ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) analyses performed in the present study and those previously reported by Herbert and Schuffert23 on Ocean Drilling Program core 165-1002C (Extended Data Fig. 1). These authors observed a slight warming several hundred years after the transition into the Holocene, between about 11.53 and 11.32 kyr b2k.Three prominent SST maxima are observed between about 11.50 and 11.45 kyr b2k. The average SST in these 50 years is 1.3 °C higher than in the 50 years before and after. These maxima are synchronous with the 11.4-ka cold event or PBO characterized by a negative excursion in δ18O and reduced snow accumulation rates in Greenland ice cores63 (Extended Data Fig. 2). The PBO coincides with the oldest of the Bond events, that is, pulses of ice rafting in the Northern Atlantic indicative of climatic deterioration64.A warm tropical response to the PBO would be supported by the lower-resolution foraminiferal SST record of Lea et al.3, which shows two data points of increased SST shortly after the end of the YD–Holocene transition. To enable direct comparison, ages in Lea et al.3 were corrected for the age difference between the sediment-colour-based YD termination midpoint in their record and in data from Deplazes et al.11. After this correction, these maxima correspond to 11.43 and 11.50 kyr b2k (Extended Data Fig. 2). Further, the SST maxima coincide with a short-lived change to lighter-coloured sediments. Hughen et al.19 described a correlation between brief North Atlantic cold events, such as the PBO, and changes in tropical primary productivity mediated by stronger upwelling that result in lighter sediments in the Cariaco Basin. Far-reaching effects of the PBO have previously been described in West Asia, with increased dust plumes being related to a southward shift of the westerlies65.The identification of the mechanisms behind a potential TNA response to the PBO is beyond the scope of this study. However, we wish to point out that high-resolution records are crucial to identify such events and to differentiate between underlying changes coinciding in time and, as in the present case, sharp signals that act on the same multidecadal timescales and can potentially be triggered by the same processes66.The effect of changing sedimentation rate on reconstructed interannual SST variability during the YD–Holocene transitionPooling proxy data into 1-year horizons establishes a constant sampling rate and thereby prevents potential effects of changing sedimentation rates. The onset of the Holocene in the Cariaco Basin sediments is characterized by a sharp decrease in sedimentation rates from 1.4 to 0.5 mm year−1 (refs. 11,19). Consequently, in the yearly pooled data, we observe a reduction in the number of values summed for each horizon (Extended Data Fig. 10b), as fewer laser spots fit into the thinner Holocene annual layers. At the same time, the mean intensity in each of these spots slightly increases, consistent with a relative increase of the contribution of haptophytes to primary production20.We have previously shown that the precision of MSI-based molecular proxy analysis is dependent on both the number of spots pooled per data point and the signal intensity in these spots54. All horizons used in the downcore record are above the established threshold of ten spots and proxy variability was shown to stabilize above this threshold9,54. However, as a decrease in the number of values per horizon might still result in lower analytical precision and contribute to higher signal variability, we corrected variability in the 2–8-year window with the estimated analytical variability (see equation (3)). With this correction, the magnitude of the described variability decreases across the record, but the trend towards higher interannual variability in the Holocene persists (Fig. 2c).Varve formation and alkenone deposition in the sediments of the Cariaco Basin during the YD–Holocene transitionComparison of elemental maps and sediment colour (Extended Data Fig. 5) shows a consistent pattern of lamination across the YD–Holocene transition that results from the seasonal interplay of precipitation, upwelling and dominant phytoplankton community composition. Darker laminae represent the rainy, non-upwelling (summer/fall) season and are enriched in Fe and Ti from terrigenous material and Ca sourced from biogenic CaCO3 produced by foraminifera or coccolithophores. Lighter laminae are characterized by high abundance of Si and correspond to the increased production of biogenic opal by diatoms during the upwelling (winter/spring) season67. This is in agreement with observations by Hughen et al.18, who described the laminae couplets in the Cariaco Basin as representing annual cycles, whereby light laminae are an indicator of high productivity associated with the winter/spring upwelling season and dark laminae are an indicator of summer/fall runoff and accumulation of terrigenous material. Deplazes et al.68 described a divergent origin of lamination for a deeper section of the YD, with light laminae being rich in calcareous and terrigenous elements characteristic for the summer season, whereas dark layers were enriched in Si and Br, indicative of diatoms and organic-walled primary producers characteristic for the more productive winter season. Such an alteration of the characteristic pattern of lamination is not observed in the late YD investigated here.This blueprint of seasonality was used to assess the seasonal behaviour of alkenones. Alkenones were deposited throughout the year, as evidenced by the fact that the number of spots containing detectable amounts of both alkenone species are not restricted to the upwelling or non-upwelling seasons but distributed across a relatively wide range of GS values to both sides of the median (Extended Data Fig. 6). Average alkenone signal intensity is higher in the non-upwelling season, pointing to a preference of alkenone producers for this season and/or to a stronger dilution of the signal in the upwelling season. In regards to the ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) SST proxy distribution in light versus dark layers, our observations are in agreement with the ability to capture the seasonal SST cycle with alkenones in sinking particles in the modern Cariaco Basin69.Effect of changing seasonality on YD and early Holocene SST records from the western TNAChanging seasonality can contribute to explaining contrasting lower-resolution SST records in the western TNA during the YD and the early Holocene. The strong warming during the YD–Holocene transition recorded in the foraminiferal Mg/Ca record of the Cariaco Basin (Lea et al.3; Extended Data Fig. 1) might be reflecting the more robust thermohaline stratification and increasingly warmer non-upwelling seasons, given the preference of Globigerinoides ruber for this season.Globigerinoides ruber (white), as used by Lea et al.3, is considered to be a dominant species in the tropics, with a relatively uniform annual distribution. However, in the modern Cariaco Basin, upwelling leads to a distinct foraminiferal community composition and seasonal turnover70, consistent with the notion of warm-water foraminifera narrowing their occurrence to the warmest season71. The relative abundance of G. ruber increases in the non-upwelling (warm) season but rarely exceeds 15%, whereas the upwelling season is clearly dominated by Globigerina bulloides72,73. Globigerinoides ruber fluxes are consistently lowest when upwelling is most vigorous, as expressed in annual minima in SST (Extended Data Fig. 9b). As upwelling during the YD and early Holocene was more intense than in the present70, the preference of G. ruber for the summer (non-upwelling) season might have been even more pronounced.The development of a stronger seasonality in the early Holocene would thus have led to a narrower temporal occurrence of G. ruber in the non-upwelling season, during which it would also be exposed to higher SST. The average SST difference between seasons obtained in our analysis can be converted into annual SST amplitude by assuming a sinusoidal curve. By doing so, we observe an increase in the seasonal amplitude of 1.5 to 1.9 °C (depending on the ramp fitted), which is similar to the warming described by Lea et al.3.This interpretation is in agreement with Bova et al.46, who observed that most Holocene climate reconstructions are biased towards the boreal summer/fall and reflect the evolution of seasonal rather than annual temperatures. As discussed above, this is probably not true for the ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) index in the Cariaco Basin, as alkenones are deposited throughout the year. The suggested weakening of summer stratification during the YD (as compared with the Holocene) might, however, explain why the lower-resolution ({{rm{U}}}_{37}^{{{rm{K}}}^{{prime} }}) records from the semi-enclosed Cariaco Basin show no or weaker warming23 than other, open-ocean, tropical YD records4, where the interplay of upwelling, freshwater input and stratification are less relevant to the SST signal. More
