Subfossil tree trunks representing Scots pine wood (Pinus sylvestris L.) were collected from lake sediments in several sites at the Arctic timberline and at the Arctic Circle (Figures 1A,B; Supplementary Table S1). The trunks were pulled to the shore, sawn into disks, and returned to the lake. The samples were air-dried for 6 to 12 months and their cross-sectional surfaces were sanded down and scanned for digital images from which the series of tree-ring widths were measured using CooRecorder 7.7 software (Figure 2) (Rydval et al., 2014; Maxwell and Larsson, 2021). In addition to visual dating, the process of cross-dating was carried out using statistical tests (Holmes, 1983) using the tree-ring width data. Moreover, the surface reflectance of the wood was determined from the scanned images using the CooRecorder (Rydval et al., 2014; Maxwell and Larsson, 2021), and the minimum blue channel light intensity (BI) of the latewood portion of each annual ring was obtained. The BI data have been shown to provide a highly skilled surrogate for the MXD of P. sylvestris tree rings (Campbell et al., 2007). New MXD (maximum latewood density) data were produced for a smaller set of tree-ring samples from living and subfossil P. sylvestris trees using the standard protocols (Gunnarson et al., 2011). Regional Curve Standardization (RCS) (Briffa et al., 1992; Helama et al., 2017), in combination with signal-free implementation (Melvin and Briffa, 2008), was used for removing the non-climatic, age-related variations from the BI and MXD data. Age-related trends in tree-ring data were identified using the “Shah” and “Hugershoff” functions fitted to a cloud of tree-ring measurements (Björklund et al., 2013). The “Hugershoff” function (Bräker, 1981) is commonly used for modeling age-related trends in the MXD series (Briffa et al., 1992; Briffa et al., 1996). This function was found to be unsatisfactory as a model for the age-related changes in BI data, and the “Shah” function was used here for modeling the BI change as a function of tree biological age. The parameterizations of these functions were solved by the Levenberg–Marquardt algorithm (Moré, 1978). The series of (detrended) tree-ring indices were obtained by dividing BI and MXD measurement values by the respective values of the “Shah” (BI_x = a + b^x + cd^x, where x is the tree age and a, b, c, and d are coefficients from modeling) and “Hugershoff” (MXD_x = ax^b e ^−cx + d) functions. Cross-dated and detrended tree-ring indices were averaged into the mean BI and MXD chronologies. Statistical comparisons between the chronologies were carried out using Pearson’s product-moment correlation coefficient. Similar to a previous study (Esper et al., 2013), the chronologies were compared only over the well-replicated periods (replication, n ≥ 10 series per year). The earliest and latest dates fulfilling this criterion in our BI composite chronology (combining the sites at the Arctic timberline and at the Arctic Circle) are 5487 BCE and 1811 CE, respectively. The common period for the MXD-based records of Matskovsky and Helama (2014) [built previously using the MXD data from the same region by Schweingruber et al. (1988), Grudd (2008), Esper et al. (2013), and Melvin et al. (2013) (https://www.ncdc.noaa.gov/paleo-search/study/17264)] and of this study with n ≥ 10 series per year was achieved for the 1823–2010 CE period. The temperature-sensitivity of the BI data was tested by calculating Pearson’s correlations between the BI chronology and monthly and seasonal temperature records. Such temperature data representing the region of Lapland and overlapping with the BI records were available from the Tornio/Yli-Tornio district. These data cover the years from 1738 to 1835 CE, with gaps (number of years with data varied from 47 to 48, depending on the monthly variable) (Holopainen and Vesajoki, 2001; Helama et al., 2004b). Calculation and homogenization of daily temperature values were carried out using the established methodology (Moberg and Bergström, 1997), which accounted for changing measurement time of observations, as described previously by Vesajoki et al. (1995) and Holopainen and Vesajoki (2001). Later, Klingbjer and Moberg (2003) suggested excluding the observations of this temperature dataset from the year 1815 CE because of potential outliers. Here, we calculated the correlations first for a full set of temperature data (Holopainen and Vesajoki, 2001; Helama et al., 2004a) and then for data from which values of that particular year (1815 CE) have been removed. For comparison, similar correlations were calculated for the MXD-based records of Matskovsky and Helama (2014).

The data of the new MXD record were transformed into a temperature scale first by 1) linear regression and then, as an alternative, 2) by setting the mean and standard deviation of the new data equal to the mean and standard deviation of the existing temperature reconstruction (Matskovsky and Helama, 2014) over the same period. Our third approach was 3) to set the mean and standard deviation of the regression-based reconstruction 1) equal to the mean and standard deviation of the existing temperature reconstruction (Matskovsky and Helama, 2014). The 117–1811 CE period was known to have a replication of n ≥ 10 and was used for the calibration of the BI chronology, over which period this record was transformed into annually-resolved temperature estimates by 1) regression, 2) scaling, and 3) regression after scaling. The data withheld from calibration were used to verify the reconstruction skill, and the MXD- and BI-based temperature reconstructions were compared with the existing temperature record (Matskovsky and Helama, 2014). Calibration and verification statistics (Tables 1, 2) calculated over the instrumental and preindustrial periods included the coefficient of determination (R_cal ²), variance explained (r_ver ²), the reduction of error (RE) (Fritts, 1976), and the coefficient of efficiency (CE) (Briffa et al., 1988). The scaled temperatures resulted in equal outcomes in the case of R_cal ² and r_ver ² for their regression-based counterparts, but negative RE and CE results indicated that scaling was not a proper approach for transforming the tree-ring values into paleoclimate estimates. Compared to scaling, multiple regression models benefit from using more than one explaining variable in the time domain. In practice, regression models explaining the temperature in a year (t) with tree-ring values in concurrent (t), previous (t−1), and forthcoming (t + 1) years can result in improved calibration and verification statistics (Briffa et al., 1988; Briffa et al., 1992; Helama et al., 2009). Improvements in R_cal ² and r_ver ² statistics also indicate reconstructed temperature amplitudes that are closer to those observed in the explained variable (Esper et al., 2005; Helama et al., 2009). In the case of our MXD and BI data, multiple regression models showed higher R_cal ² and r_ver ² than scaling or regression models with one explaining variable. Combined with scaling, the multiple regression models resulted in positive RE and CE outcomes and were chosen as our paleoclimate reconstructions. These reconstructions showed higher temperature amplitudes due to both higher R_cal ² and scaling (here, subsequent to regression) (see Esper et al., 2005). Confidence intervals were estimated for the reconstruction using the same calibration method. Regarding the MXD data, a set of 50 tree-ring values was randomly selected from the 1823–2010 CE period and used for the calibration and production of temperature reconstructions over the 117–2010 CE period. A similar approach was applied to BI data, but a set of 100 tree rings was randomly selected from the 117–1811 CE period and used for calibrating and producing temperature reconstructions over the full preindustrial period (5487 BCE –1811 CE). This process was repeated n = 1,000 times and the 97.5 th and 2.5 th percentiles from the n = 1,000 reconstructions (calculated separately for MXD and BI data) were adopted as the confidence intervals for each annual estimate of the reconstructions.

Linear regression was used to estimate the slope of the temperature trends. This statistic, with the standard error value for the slope, was calculated for the temperature estimates from our BI data over the full preindustrial period (5487 BCE–1811 CE) and for sedimentary summer temperature records (Sundqvist et al., 2014) over the corresponding interval between 7,440 and 140 cal BP. These records represent various types of sedimentary proxy data (e.g., pollen, chironomids, and diatoms) and across the circumpolar Arctic sites (https://www.ncdc.noaa.gov/paleo-search/study/15444) (Sundqvist et al., 2014). Compared with tree-ring records, the data in these series were of much coarser resolution (temporal resolution of c. 120 years for Fennoscandian records, see Table 1 in Sundqvist et al., 2014), dated mainly by means of ¹⁴C ages and resulting age-depth models and, as a result, not equidistant. Only quantitative reconstructions indicative of summer (warm season) temperatures transformed into the centigrade scale (°C) and given in years converted to the calendar timeline (expressed as cal BP ages i.e., years before 1950 CE) were included. If the site had more than one record with the same proxy type, same climate variable, and same geochronology (age estimates), the available temperature values were averaged before trend estimation to avoid overrepresentation of that one proxy source. The one-sample Wilcoxon signed-rank test (Conover, 1980) was used for comparing the median of a set of temperature slopes (obtained for the sedimentary proxy data) against the slope obtained for the tree-ring temperature reconstruction. The 95% confidence intervals for the medians were obtained from bootstrapping, repeating the resampling n = 9,999 times (Hammer et al., 2001). Our tree-ring–based temperature record was compared with those originating from other Fennoscandian proxy data [excerpted from the database of Sundqvist et al. (2014)]. These data show mean temperatures negatively correlating with site latitude and standard deviations of c. 0.4°C (latitude-independent), which is similar to the standard deviation of our tree-ring–based temperatures over the preindustrial period (see following section). For visual comparisons, the individual Fennoscandian sedimentary proxy records were normalized to a mean of zero and a standard deviation of c. 0.4°C over this period. Then, the transformed data of the three proxy types (pollen, chironomids, and diatoms) were regionally averaged for 100-year periods (2010–1911 CE, 2000–2001 CE, 1990–1891 CE.) and spline functions acting as “low-pass” filters (Cook and Peters, 1981) representing a 500-year rigidity were applied to those three records for illustrating their temperature variations on centennial and millennial timescales. The BI (5487 BCE–1811 CE) and MXD (1823–2010 CE) records were combined into one temperature record by placing them sequentially. This record was averaged for 100-year periods and smoothed to correspond to the sedimentary proxy records for the purpose of this comparison. The sedimentary records (pollen, chironomids, and diatoms) were correlated with tree-ring records, and Monte Carlo tests (Macias-Fauria et al., 2012) were used for determining the exact p-values.

Our tree-ring data show non-climatic (age-related) values for both the MXD and BI data (Figures 1C–E). Trees from the Arctic timberline show higher BI values than those from the Arctic Circle as is demonstrated by age-related growth trends with differing levels of BI values (Figure 1F). This is consistent with the expectation to find lower wood density and thus higher BI (due to their negative correlation) under colder (timberline) climates. Using state-of-the-art dendroclimatic methods (Briffa et al., 1992; Briffa et al., 1996; Melvin and Briffa 2008; Melvin et al., 2013; Helama et al., 2017), these trends were removed from the tree-ring data from the Arctic Circle and timberline by their respective age-models. The resulting BI chronologies portrayed variations indicative of climate variability on inter-annual to longer scales (Figure 3). The chronologies from the Arctic timberline and Arctic Circle were generally similar in their growth characteristics. Overall, these BI chronologies correlate with r = 0.547 with each other over the well-replicated periods, thus demonstrating a common signal in their mean records with consistent latewood growth patterns (see Table 3 for correlations between the BI chronologies over multiple periods). The mean composite BI record (all trees included) compares also well with the existing MXD record (Figures 4A,B) as they correlate with r = −0.586 over their common interval 117–1811 CE. The mean BI data are also found to correlate with the instrumental summer (June–August) temperatures recorded in Lapland, with coefficients comparable (Figure 5A) and even slightly higher (Figure 5B) than those observed between the MXD data and June–August temperatures, despite the much lower replication of BI data over this instrumental period (Figure 5C).

Transforming the new MXD and BI records into temperature data using statistical calibrations over the instrumental (Figure 4C) and preindustrial periods results in statistically significant correlations and positive verification statistics with data withheld from regression-based fitting (Tables 1, 2). Produced by linear scaling, the reconstructed values show higher temperature amplitudes than regression (Figure 6A), which accords with the paleoclimate literature on this method (Esper et al., 2005), but the scaling-based temperature values do not verify satisfactorily over the withheld data (Tables 1, 2). When scaling is performed after multiple regression, the verification statistics remain positive with preserved (i.e., higher) temperature amplitudes. Here, we also compare the temperature amplitudes with existing data from sedimentary proxy records (Sundqvist et al., 2014) over the same preindustrial period. These records show a mean temporal resolution of c. 120 years (Sundqvist et al., 2014) and a standard deviation of c. 0.4°C during the preindustrial period (Figure 6B). These values are nearly identical to those based on our tree-ring record, which exhibits temperature variability with a standard deviation of 0.4°C on similar 100–150 year time scales (Figure 6A). Combining scaling with multiple regressions resulted in improved calibration and verification statistics (in comparison to simple scaling and univariate regression), including higher R_cal ² and r_ver ² and, accordingly, in higher reconstructed temperature amplitudes, which were compatible with the temperature evidence of sedimentary proxy data. These examinations further validate our BI data to constitute a reconstruction with temperature low-frequency amplitudes comparable to those of sedimentary records.

Our new high-latitude, high-resolution temperature reconstruction from BI data shows a negative temperature trend over mid- and late-Holocene times (Figures 4D,E). Computed over the well-replicated periods, this trend suggests that summer temperature has decreased by 0.13°C (±0.01°C) per millennium between 5486 BCE and 1811 CE (Figure 4F) and it parallels the coeval decrease observed in insolation values demonstrating the coterminous change in orbital forcing (Figure 4G). Over the Meghalayan Stage and Common Era, the trends are stronger showing a cooling of 0.18°C (±0.01°C) and 0.69°C (±0.04°C) per millennium, respectively (see Figure 4F), thus outweighing the previous MXD-based cooling estimate of 0.31°C (±0.03°C) per millennium for the same region over the Common Era (Esper et al., 2013). Apparently, these estimations are sensitive to shorter, sub-Milankovitch climate variations (e.g., low temperatures over the latest CE centuries due to the “Little Ice Age”) that are superimposed on the predictable long-term change. In any case, similar estimates of the temperature slope are obtained for a pollen record (Seppä and Birks 2001) central to our tree-ring sites showing a cooling of 0.11°C (±0.02°C), 0.21°C (±0.04°C), and 0.64°C (±0.12°C) over the full preindustrial period, Meghalayan Stage, and Common Era, respectively (see Figure 4H).

Differing temperature values are, however, observed when comparing the mid-Holocene temperature reconstructions from BI and MXD data. Despite the high correlation between the two reconstructions at year-to-year scales over this period (r = 0.749), the mean level of the BI-based temperatures (14.5°C) exceeds the MXD estimates (12.9°C) by 1.6°C as observed over their common interval 4260–4071 BCE (Figures 4D). Overall, the MXD-based temperatures appear anomalously low in the context of the-mid Holocene temperatures, as inferred from the BI data (for other sedimentary proxy records, see the following section). These findings suggest low-frequency decoupling of the BI and MXD data and incomparability of their paleoclimatic signals on orbital scales whereby the MXD (unlike BI) data remain insensitive to similar long-term temperature changes. Despite this, a high correlation between the mid-Holocene BI and MXD records, at year-to-year scales, provides strong evidence for their reliability as indicators of high-resolution, high-frequency climate variability.

Phytogeographical data show altitudinal and latitudinal retreats of the Pinus sylvestris timberline in the region of NW Finnish Lapland (for sites, see Figure 1B) since 2000 BCE, corresponding to 300 m downslope and 85 km toward the Equator (Figures 4I,J). Postglacial isostatic land uplift was found to explain no more than ca. 20 m of this change over the past 4,000 years (Helama et al., 2004b) and was not considered here. Using a temperature lapse rate of 0.5°C per 100 m (Laaksonen 1976), this change (300 m) translates into a summer temperature change of 1.5°C over the past 4,000 years. As a result, these estimates agree to show a cooling trend of 0.375°C per millennium at least since 2000 BCE, which markedly exceeds the trend seen in our BI reconstruction with a cooling of 0.18°C per millennium over the same interval (see Figure 4F). The trend obtained for our BI reconstruction agrees markedly with the slope of the trends obtained for all the circumpolar Arctic sedimentary proxy records over the full preindustrial period (Figure 7). More contrasting estimates of long-term cooling are, however, observed for the subsets of the Arctic records. The trends for pollen records and in Fennoscandian sites (where pollen data dominate the number of records) are generally stronger than for chironomid records and in other regions. As a result, the slope of the trend estimated for our BI-based temperatures is more gradual than the median of pollen and that of all Fennoscandian records but stronger than the median of all non-Fennoscandian and that of chironomid records. Likewise, the tests for comparing these medians against the slope of our BI temperatures (Supplementary Table S3) suggest comparability of the new BI-based cooling in a wider circumpolar Arctic context, excluding the pollen estimates.

Combined, the Fennoscandian proxy records have many similarities but also some unique features, with high and statistically significant correlations against tree-ring data, especially with pollen and diatom records (Figures 6C–F). Here, we neither aimed to check for their further inconsistencies nor attempted to evaluate their veracity but calculated the mean temperature records for different proxy types. In keeping with the Milankovitch scale forcing, these records show varying but generally declining summer temperatures over the mid- and late-Holocene times until the first millennium BCE, after which the tree-ring and especially the diatom and chironomid data indicate a transient rise in temperatures, followed by a subsequent temperature decline toward the coolest centuries of all records, in accordance with the “Little Ice Age” climate as frequently recorded in several types of proxy evidence around the circumpolar Arctic sites (Kaufman et al., 2009; Nicolle et al., 2018; Helama et al., 2021). Clearly, different types of proxy data come with their own characteristics. As evident from the circumpolar Arctic sites, pollen data appear to contain stronger trends of cooling than other proxy records (Figure 7). In the study region, discrepancies between the pollen and other sedimentary proxy records have recently been related to inaccuracies in the models transforming the pollen data into temperature estimates, to climatic disequilibrium of vegetation migration patterns, and to site-specific factors (Luoto et al., 2014b; Shala et al., 2017). Chironomid-based records, however, may be affected by limnic factors (Luoto et al., 2014a).