Orbital forcing of tree-ring data

Tree-ring data and spatial coherence.

We collected core samples from living P. sylvestris trees growing at lakeshore and inland (that is ten or more metres distance from lakes) microsites, and disc samples from submerged logs in northern Finland and Sweden (Supplementary Table S1 and Fig. S2). MXD data were derived from high-resolution density profiles using X-ray radiographic techniques²⁸ (Fig. 1). Within and between-site coherence of the northern Scandinavian MXD network has been assessed using a total of nine data sets from living trees—of which three (Ket, Kir, Tor) are additionally subdivided into lakeshore and inland subsets—and 14 data sets from subfossil lake material. We calculated Pearson correlation coefficients among living-tree chronologies over the 1812–1978 common period (r_MXD = 0.72, r_TRW = 0.58; Supplementary Table S2 and Fig. S3), and over varying periods of overlap (AD 700–1600) between subfossil MXD chronologies (r_MXD = 0.71; Supplementary Table S3 and Fig. S4) to estimate data homogeneity throughout space and time. To ensure signal homogeneity, we considered MXD data from only lakeshore sites together with the subfossil material discovered from the lakes for the final reconstruction (N-scan). The record integrates 587 high-resolution P. sylvestris MXD measurement series.

Chronology development and assessment.

Various combinations of living-tree and subfossil MXD data were produced to assess the sensitivity of the long-term N-scan record to detrending methodology and microsite conditions. Application of negative exponential and RCS detrending techniques revealed substantial changes in retained low-frequency variance and sensitivity of twentieth-century trends to density differences between living-tree sites ranging from ~0.002 to 0.010 g cm⁻³ over the first 200 years of tree growth (Supplementary Fig. S5). Sensitivity of increased (decreased) twentieth-century chronology levels was assessed using MXD data from only lakeshore (inland) microsites in long-term RCS runs (Supplementary Fig. S6). N-scan characteristics were detailed by calculating 95% bootstrap confidence ranges, chronology age and replication curves, and interseries correlation and expressed population signal statistics (Supplementary Fig. S7). We analysed the RCS detrended N-scan data by classifying the measurement series into age classes ranging from 1–10 years to 201–210 years and calculating 100-year spline filters for each of these classes (Supplementary Fig. S8). This procedure provided insights into the coherence of long-term trends retained in (typically noisier) juvenile and (typically less replicated) adult tree-ring data. N-scan trend behaviour was additionally assessed by analysing the persistence of low-frequency variability in tree-ring parameters, indicating that only the MXD data preserved substantial variance on millennial timescales (Supplementary Fig. S9).

Proxy calibration and JJA temperature reconstruction.

The MXD climate signal was assessed using Pearson correlation coefficients between the lakeshore subsets Ket-L (r = 0.74), Kir-L (r = 0.75) and Tor-L (r = 0.74) and mean JJA temperatures recorded at the global historical climatology network stations Haparanda, Karasjok and Sodankyla over the 1876–2006 common period. Running correlations were applied to analyse the temporal characteristics of the signal revealing reduced coherence among the station records as well as between the station and proxy data centred in the 1910s (Supplementary Fig. S10). The long-term N-scan record integrating lakeshore and subfossil MXD data correlates at 0.77 (r² = 0.59) with regional JJA temperatures. We transferred this record into a JJA temperature reconstruction using ordinary least square regression with MXD as the independent variable. This approach provides conservative estimates—owing to the reduction of variability caused by unexplained variance²⁹—of pre-instrumental climate variability and derived long-term trends. Split-period calibration/verification statistics³⁰ with early and late r² (0.-57–0.61), reduction of error (0.57–0.59), coefficient of efficiency (0.50–0.54) and full period Durbin–Watson (1.75) statistics were applied to validate the reconstruction. N-scan confidence intervals were calculated considering the standard error (±2 × root mean square error) derived from verification against instrumental JJA temperatures over the early 1876–1941 period and a bootstrap confidence range derived from resampling the MXD data 1,000 times with replacement. A total of 2,000 MXD chronologies derived from randomly drawn subsets of the N-scan record were developed to test the influence of reduced sample replication, typical to earlier periods of the long-term reconstruction, on the calibration results. These tests revealed that the transfer model remains robust (r > 0.70) down to a replication of ten MXD measurement series (Supplementary Fig. S11). Extreme cool and warm summers (decades and centuries) since 138 BC are expressed as deviations from the 1951 to 1980 mean (Supplementary Table S4) and millennial scale JJA temperature trends estimated by calculating a linear ordinary least square regression over the 138 BC–AD 1900 period (Fig. 2). The robustness of the regression slope (−0.31 °C per 1,000 years) was tested by reducing the length of the regression period stepwise at both ends by 100 years to derive a 95% confidence range (±0.03 °C) of the millennial scale trend.

CGCM Holocene simulations.

Spatial patterns of JJA surface air temperatures derived from multimillennial ECHO-G (ref. 7) and ECHAM5–MPIOM (ref. 8) CGCM runs forced with and without long-term insolation changes were analysed to estimate low-frequency temperature trends throughout the Northern Hemisphere extratropics (Supplementary Fig. S12). ECHO-G is one of the coupled atmosphere–ocean models considered in the Intergovernmental Panel on Climate Change Fourth Assessment Report and was ranked among the best five models in simulating the mean patterns of surface atmospheric circulation and precipitation⁴. It integrates the atmospheric ECHAM4 model with a horizontal resolution of 3.75 × 3.75 degrees and 19 vertical levels, and the oceanic HOPE model with a horizontal resolution ranging from about 2.8 × 2.8 to 0.5 × 0.5 degrees towards the Equator and 20 vertical levels including a thermodynamic sea-ice model. To avoid artificial climate drift in the very long (7,000 years) climate simulations used here, a flux adjustment was applied to the atmosphere–ocean coupling. The second transient Holocene simulation⁸ consists of the spectral atmosphere model ECHAM5 run at truncation T31, corresponding to a horizontal resolution of a 3.75 × 3.75 longitude–latitude grid, with 19 vertical hybrid sigma pressure levels and the highest level at 10 hPa. It integrates the land-surface model JSBACH including a dynamic vegetation module, has been coupled to the ocean GCM MPIOM run with 40 vertical levels (30 levels within the top 2,000 m) and includes a zero-layer dynamic-thermodynamic sea-ice model with viscous-plastic rheology. No flux correction has been applied to this CGCM.