The δ¹⁸O and δ²Hn tree-ring cellulose data were analyzed from six sites across Spain, the European Alps, Iceland, and Greenland (Figure 1A; Table 1). While the datasets from Iceland and Greenland are unpublished, the Alpine arc data comes from the publicly available “Alpine Holocene Triple Tree Ring Isotope Record” (AHTTRIR) (Arosio et al., 2022), which comprises 202 trees at the tree line in the Swiss, Austrian, and Italian Alps, covering 46.03–47.03°N and 7.55–15.25°W. This dataset spans the last 9,000 years up to 2015 CE with a 5-year resolution. Analyses indicate a very low correlation (r < 0.2) between tree-ring width and isotopes in the Alpine dataset, confirming that favorable growing conditions (wider rings) do not significantly influence the 5-year average δ¹⁸O and δ²Hn values (Arosio et al., 2021). Greenland samples, with exceptionally narrow rings (mean 0.123 mm), were sectioned into 10-year blocks, while Icelandic wood was cut into annual slices. All samples underwent cellulose extraction (Ziehmer et al., 2018) and isotope analysis (Loader et al., 2015). The network spans six sites from 37°N to 80°N (Figure 1A), featuring larch and non-larch conifers. Samples from Spain and Iceland cover four centuries, the Alps span nine millennia, and Greenland’s wood extends over four centuries but predates the ¹⁴C dating limit by over 50 millennia (Table 1).

For the combined online measurements of δ¹³C, δ ¹⁸O and δ ²Hn, the samples were subjected to pyrolysis at temperatures exceeding 1,450°C in a Flash HT elemental analyzer to convert the cellulose to carbon monoxide and hydrogen (Leuenberger and Filot, 2007) following the method described in Filot et al. (2006). This high temperature is critical for decomposing the cellulose quantitatively into hydrogen and carbon monoxide, with the residue (partial carbon of the cellulose) (Leuenberger and Filot, 2007) collected in a graphite crucible. The resulting gases were separated via a 5 Å molecular sieve and analysed using a Thermo Delta XP or an Elementar/Isoprime mass spectrometer for δ¹³C and δ¹⁸O on CO and δ²Hn on H₂.

δ ²Hn values are determined using an online equilibration technique described by Filot et al. (2006) and extended to a triple isotope method (COH) described by Loader et al. (2015). The specific procedure involves equilibrating the cellulose samples with water vapour of known δ ²H values (close to VSMOW, variable in the range of −2 to +2‰ over years depending on the aliquot of internal water standards, calibrated against VSMOW with a precision of 0.5‰) at 110°C for 10 min using a continuous flow of helium carrier gas at 80 mL/min. The rate of equilibration was checked with different calibration waters in the range (−429.1 (Dome C) to +1 (Meerwasser) ‰) to cover a wide range of δ ²H values as already used by Filot et al. (2006). The samples and standards are stored in a desiccator when not being analysed. This equilibration method was cross-calibrated with the nitration method using cellulose standards with known δ ²Hn (Filot et al., 2006).

The mass spectrometric analysis involved the determination of isotope ratios of hydrogen and carbon monoxide, which were sequentially detected and measured after elution. This process included working gas pulses for both hydrogen and carbon monoxide, yielding raw delta values, which were converted in per mil (‰) deviations relative to Vienna Standard Mean Ocean Water (VSMOW) for hydrogen and oxygen as well as carbon based on the calibration procedure described below.

Calibration and quality control were integral to our methodology. A two-point calibration of the raw delta values obtained against our working gas was performed for all isotope ratios (Table 2). For hydrogen, we used Merck cellulose (−34.94‰ δ²H_total, and −70.9 δ²H_n) and IAEA-CH7 polyethylene standards (−100.3‰ for both δ²H_total, and δ²H_n as no exchangeable hydrogen is available). This calibration approach allows us to determine the offset and contraction/stretching of the scale relative to the international VSMOW-SLAP scale. These scale correction terms are independent of whether there are available exchangeable hydrogen atoms. To check for drift in the reference gas and to evaluate precision over time, we additionally used Wei Ming 101 cellulose (from WEI MING Pharmaceutical MFG) and/or IAEA-C3 cellulose (form IAEA) and IAEA-C6 sucrose (form IAEA) standards (Rozanski et al., 1992; Coplen et al., 2006). This led to the evaluation of the machine performance and the implementation of the necessary corrections to a linear drift. Furthermore, the inclusion of IAEA-CH7 specifically enabled quality control of the non-exchangeable hydrogen measurements, while the cellulose standards accounted for the exchangeable portion. For δ¹⁸O and δ¹³C, we used Merck and IAEA-C3 for the calibration, the corresponding values are given in Table 2. The system’s stability allowed for consistent and reliable data within the standard precision limits of ±0.3‰ for oxygen and ±3.0‰ for hydrogen.

Tree-ring δ ¹⁸O and δ ²H data from three Spanish sites (Figure 1A; Table 1) were generated during the ISONET project (Treydte et al., 2006) and followed the methodology for δ ¹⁸O detailed in Treydte (Treydte et al., 2014). The δ ²Hn measurements were conducted using the technique described above. All samples in this study were measured for non-exchangeable, C-bound hydrogen from tree-ring cellulose (Filot et al., 2006; Loader et al., 2015).

All examined trees were absolutely dated with only one exception, the tree from Greenland, which was estimated to be older than 50,000 years. The fragmented remains of this tree were discovered within a glacier moraine at an elevation of 331 m asl, located more than 90 km away from the open sea.

The δ ²Hn and δ ¹⁸O of the Alpine dataset were corrected for the effect of ice volume with a gradient obtained for δ ¹⁸O (Fleitmann et al., 2009; Affolter et al., 2019) and converted to δ ²H by multiplying it by a factor of eight, resulting in a gradient of −0.064‰ per meter of sea-level rise (Rozanski et al., 1982); as previously done for speleothem records from Switzerland (Affolter et al., 2019). The effect of ice volume correction on the Alpine database for δ ²Hn and δ ¹⁸O values is shown in Supplementary Figure S5. The other datasets do not require a correction as they are from the late Holocene period, where the correction is ineffective (Supplementary Figure S5). The Greenland sample has not been corrected for temperature because the growth temperature is unknown.

The oxygen values of larch and non-larch trees are known to have similar values if they come from the same region and time, while those of hydrogen are known to be lower in larch than non-larch trees (Arosio et al., 2020b). Thus, we corrected all larch values based on the difference between the larch and non-larch species in the Alpine database (Figure 2, Supplementary Figures S2, S4-1, S5). The correction factor is −48.1‰.

Monthly temperature and precipitation data were extracted from the Climatic Research Unit’s CRU TS4.03 dataset (Harris et al., 2020) for the period 1905-2002 at each study site.

Summer (JJA) and winter (DJA) seasonal means deviation for monthly precipitation hydrogen (δ ²Hp) and oxygen (δ ¹⁸Op) stable isotope values were calculated per site by extrapolating the monthly values from the gridded data defined by Bowen et al. (2005), using each site’s latitude, longitude, and altitude. The modelled data are based on precipitation measurements from 1960 to the present (Bowen et al., 2005).

To assess the variable’s dependency, a linear regression has been calculated (Figures 1, 2, 3; Supplementary Figures S1–S5) and associated with the P value and the Statistical power. The Statistical power is defined as one minus the probability of type II error (Greenland et al., 2016). The choice to use the average and standard deviation to represent the site values is based on the normal distribution of the isotopic values (Supplementary Figures S3A, B).

To account for the uncertainty linked to the average of each region for δ ²Hn and δ ¹⁸O, we conducted a Monte Carlo simulation with 10,000 iterations whereby the linear model was fitted to δ ²Hn and δ ¹⁸O incorporating the uncertainty by randomly adding the standard deviation of both isotopes.

The inferred mean linear regression is 8.03 for larch and 6.64 for the non-larch group (Supplementary Figure S1.1). The same analysis was carried out for all regions together but with the correction for larch and not including the Greenland sample, with a mean linear regression of 8.19 (Supplementary Figures S1.2). To show the uncertainty generated from the Monte Carlo simulation, we calculated the distribution density of the simulated slopes and intercepts (Supplementary Figures S1.2B, C), simulated values and all linear interpolations (Supplementary Figure S1.2C).

The relationships between the isotope data (δ ²Hn and δ ¹⁸O) and geographical/climatic variables were analysed using two approaches: 1) linear regression between each isotope and each individual variable (Supplementary Figure S2, Supplementary Figure S4.1) and 2) multivariate regression analysis incorporating all variables simultaneously, excluding the Greenland sample (Supplementary Table S1).

All the data analysis was performed in R version 4.3.2.

Our network is mostly focused over Europe and covers a large part of the Northern Hemisphere from the south of Spain to northern Greenland, across a range of time periods, and it shows that the slope of the linear interpolation of the regional averages of tree-ring δ ¹⁸O and δ ²H parallels that of the GMWL (Figure 1A) and the LMWL (Figure 2B). It is important to note that the GMWL has remained stable throughout the Holocene, providing a consistent reference for isotopic studies across this geological epoch (Rozanski, 1985). In contrast, we observe no relationship between δ ¹⁸O and δ ²Hn at the single site level (Figure 1B). An analysis of the Alpine dataset shows that the isotopes present a normal t-distribution (McCarroll and Loader, 2004) (Supplementary Figure S3), supporting the use of mean values and standard deviations for the representations of our datasets.

Larch was shown to have δ ²Hn values different from those of non-larch conifers, in agreement with the previous study (Arosio et al., 2020b). Thus, δ ¹⁸O and δ ²Hn data were initially analysed separately. The linear regression of the means of their δ ¹⁸O and δ ²Hn values was statistically significant (Figure 1A, p < 0.05, R² = 0.98 and slope = 8.13, power = 0.80, n = 4). This was confirmed by a Monte Carlo analysis, which includes the uncertainty of averages through the use of standard deviation. It produced a regression δ ¹⁸O vs. δ ²Hn with a similar slope (mean = 8.03, first quartile = 6.85 and third quartile = 9.01) (Supplementary Figure S1.1C). Importantly, the 8.13 slope value is close to that of the GMWL (Craig, 1961) (8.0), and to that of the LMWL (7.9 calculated with the modelled precipitation values) (Figure 2A). The linear regression of the δ ¹⁸O and δ ²Hn mean values of the non-larch was also significant but only marginally (p = 0.054, R² = 0.76 and slope = 7.46, power = 0.54, n = 5), confirmed by Monte Carlo simulations that produced a similar slope value (mean = 6.64, first quartile = 4.67 and third quartile = 8.23) (Supplementary Figure S1.1B). Thus, the linear regression analysis of the δ ¹⁸O-δ ²Hn mean values of all the sites was statistically significant (Figure 1A) and aligned with the GMWL (Figure 1A) and the LMWL (Figure 2). However, the intercepts of the larch (−358‰) and of non-larch (−275‰) linear regressions differed remarkably from that of the GMWL (+10) (Craig, 1961), because of the isotopic fractionation occurring in trees. Their different intercepts could be explained as the δ ¹⁸O in cellulose is ∼27‰ higher than that of leaf/precipitation water, while the δ ²Hn was shown to be −171‰ after photosynthetic and +158‰ at post-photosynthetic steps (Lehmann et al., 2022) and differ in larch and non-larch conifers (Arosio et al., 2020b).

The regression lines of larch and non-larch trees have very similar slopes, indicating that the difference between larch and non-larch values is roughly constant among different sites. The use of the two species together is possible but only by correcting the larch samples (Figure 2B, Supplementary Figure S1.2). After correction and considering all the regions together, the linear regression produced a slope value of 8.2 (p < 001, R² = 0.79, n = 8, power = 0.91, mean of 8.19 with a Monte Carlo analysis) (Figure 2B; Supplementary Figure S2). This value is similar to that of the regressions adjusted for the simulated isotopic precipitation values of the same sites (Figure 2A), which are 7.9 and 8.0 for the LMWL and the GMWL (Craig, 1961), respectively. Of note, we excluded the Greenland sample for the statistical analysis of Figure 2B, because it lacks an exact dating and its sea-level correction; moreover, its rather low values may affect the linear interpolation. Nevertheless, it remains an interesting verification sample.

Figure 2B shows that the mean values of all sites are well aligned on the regression line, except Spain1, which has higher δ ²Hn values than Spain 2. This difference cannot be explained by climatic or geographical effects (Table 1), or isotopic value of precipitation (Figure 2A), and we attributed the increased δ ²Hn values to physiological stress at Spain 1 that may increase autotrophic metabolism and affect isotopic fractionation (Cormier et al., 2018) (Figure 4).

After excluding the Greenland sample, we found that the relationship between the mean regional summer temperatures and the mean regional isotopic values for both δ ²H and δ ¹⁸O was significant (p < 0.01 R² = 0.85 and 0.83 for both δ ¹⁸O and δ ²Hn respectively) (Supplementary Figures S4-1A, B). In contrast, summer precipitation amounts do not seem to be a representative parameter of the regional tree-ring δ ²Hn and δ ¹⁸O, as the correlation is not significant (Supplementary Figures S4-1C, D). We conclude that temperature was the primary factor driving the variations observed in both δ ²H and δ ¹⁸O mean values across the different regions.

In addition, we found a strong significant correlation between the regional tree-ring δ ¹⁸O and δ ²Hn (larch corrected) and the modelled summer isotope precipitation values (Bowen et al., 2005; IAEA/WMOA, 2006; Bowen, 2012), while the correlation with winter isotope precipitation values was weaker, especially for δ ¹⁸O (Figure 3). The correlation between tree-ring δ ²Hn and precipitation δ ²Hp is lower than that between tree-ring δ ¹⁸O and δ ¹⁸Op (Figure 3). This difference is likely caused by the higher biological fractionation of δ ²Hn, which can dampen the seasonal climate signal. However, since hydrogen and oxygen come from the same water sources, the seasonal δ ²H should reflect the summer seasonality as δ ¹⁸O does. The correlations between tree-ring isotopes and the summer temperatures are similar(δ ¹⁸O) or higher (δ ²Hn) than those between model precipitation isotopes and temperatures (Supplementary Figures S4-1, S4-2).

The stronger summer δ ¹⁸O signal is probably due to the high altitude (>1,000 m) (Spain and the Alps) or the high latitude (Greenland and Iceland) of the sampling sites, where tree growth occurs only in the summer season (Moser et al., 2010) and depends on δ ¹⁸O of summer precipitation (Szymczak et al., 2020). This contrasts with lowland broadleaved trees, the leaf water of which reflects the winter δ ¹⁸O precipitation stored in the soil (Allen et al., 2019). Note that the Greenland sample was excluded from this analysis since it comes from another interglacial stage where temperatures and precipitation amounts are unknown.

The relationship of δ ²Hn and δ ¹⁸O with latitude showed significant negative regression slopes (p < 0.01 for both isotopes, R² = 0.88 for δ ¹⁸O and R² = 0.81 for δ ²Hn) (Figures 1A, 2A; Supplementary Figures SA, B). This agrees with the known temperature-controlled isotope fractionation in precipitation caused by poleward air mass transport that carries lower δ ²Hn and δ ¹⁸O values at higher latitude sites (Dansgaard, 1964). The altitude effect was not statistically significant (R² = 0.4, p = 0.08 for δ ¹⁸O; R² = 0.4, p = 0.09 for δ ²Hn) (Supplementary Figures S2C, D), possibly because of the limited altitudinal range of our data (Kern et al., 2014).

A multivariate regression analysis confirms that a key driver of spatial δ ¹⁸O and δ ²Hn variation is summer temperatures and that the δ ²Hn values are also significantly related to δ ²Hn summer precipitation, latitude, and altitude (Supplementary Table S1).

In the sections above, we considered the mean values of the different regions, but it is also important to analyse the isotopic variation in the individual sites. We focused on our largest database, the Alps1 sites. In the larch and non-larch trees of the dataset, the correlations between the two isotopes were statistically non-significant (Figure 1B). The low correlation (R² < 0.1) remained even after the normalisation of all trees (Supplementary Figure S3A), and also when considering the mean value of the different sites (R² < 0.1) (Supplementary Figure S3B). Thus, the low correlation is intrinsic to the trees and is not due to the merging of trees from the same region, with mean R² < 0.25 for all larch and R² < 0.2 for all cembra pine (Supplementary Figure S3C).

The Greenland sample was found at an altitude of 313 m a.s.l and 90 km from the ocean, indicating in situ growth during an interglacial phase when the climate in Greenland was warm. Now there are no trees nearby. Nevertheless, the Greenland sample falls remarkably close to the linear regression line for the region (shown by the red dot in Figure 2B).

Present work analyses six sites covering a large part of the Northern Hemisphere from the south of Spain to northern Greenland, across a range of time periods. Our data show that the δ ¹⁸O and δ ²Hn means of regions exhibit a covariance that follows the GMWL, indicating that the tree-ring isotopes capture the regional value of precipitation isotopes, although with a difference between larch and non-larch conifers (Figure 1A).

The consistent offsets of the larch and non-larch linear regressions from the GMWL across the entire dataset suggest a nearly constant difference in the isotopic fractionation of the two groups during the processes of cellulose formation (Figure 1A). This suggests that large-scale systematic physiological effects on isotope fractionation do not occur over the large spatial and temporal range of our data, although there are exceptions, such as the deviation observed at Spain 1 attributed to metabolic fractionation (Figure 2B).

From the regression of the means of δ ²Hn and δ ¹⁸O, using all the samples corrected, we obtained a slope value close to 8, which is very close to the slope of the GMWL (Craig, 1961). Similar recent study reported slightly different results with a δ ²Hn and δ ¹⁸O slope value of 10, higher than the expected 8.0 for the GMWL or 7.8 for the LMWL (Supplementary Figure S9 in the appendix of Vitali et al. (2022)). The work was based on the ISONET dataset (Vitali et al., 2022) that contains conifers and oak trees of high and low altitudes and different δ ²Hn measurement systems that may cause variability in tree-ring δ ¹⁸O and δ ²Hn (Kern et al., 2014; Allen et al., 2019). This may explain the difference from our study, which contains only high altitude/latitude conifers and uses the same δ ²Hn measurement system. The precipitation signal appears coherent across our network and reflects the regional summer isotopic precipitation signal (Figure 2).

The δ²H and δ¹⁸O values of local precipitation are primarily influenced by regional-scale processes (Rozanski et al., 1982), which include geographical and climatic factors such as temperature, precipitation, latitude, continentality, and altitude (Craig, 1961; Ingraham, 1998). Given that this study employs what is likely one of the largest datasets available for continental-scale isotopic analysis, it enables a deeper investigation into the geographical and climatic factors driving tree rings isotopic variation. Our findings demonstrate that regional tree-ring δ²Hn and δ¹⁸O mean values reflect these continental precipitation isotopic variations, with temperature and latitude emerging as dominant drivers. The seasonal investigation points to regional tree-ring δ²Hn and δ¹⁸O values as predominantly reflecting summer precipitation isotopic patterns, suggesting that temperature and latitude significantly influence spatial variation. The globally coherent nature of the water cycle during the Common Era (Konecky et al., 2023) further supports the capacity of tree-ring δ²Hn and δ¹⁸O measurements to capture geographic isotopic shifts across different time periods.

This finding may be important for archaeological studies (Bernabei et al., 2019), as well as for tracing the origin of Arctic driftwood (Hellmann et al., 2013). However, this is partially limited by the difference between larch and non-larch conifers and possibly other species. Another potential confounding factor is that at high latitudes and under certain climatic conditions, the active permafrost with low δ ²H and δ ¹⁸O values can contribute to the tree water supply (Kirdyanov et al., 2024).

The main limitation of this study is the small and variable sample sizes across regions due to the scarcity of δ²Hn measurements. Strengths include wide spatial dataset coverage, exclusive use of conifer species, and consistent δ¹⁸O and δ²Hn measurement systems. Notably, an ancient Greenland sample, excluded from statistical analyses due to its uniqueness, was used for verification.

We found an absence of correlation if a single site is analysed; this absence of correlation is intrinsic to the trees and is not due to the merging of trees from the same region (Figure 1B). This cannot be explained by the relatively small spatial range or variation in the LMWL values in the European Alps (between 7.84 and 8.29) (Flaim et al., 2013; Hürkamp et al., 2019) or by the different sampling sites or offsets between trees. More likely, it is due to temporal variations of δ ¹⁸O and δ ²Hn in response to climate and weather conditions, which include leaf transpiration, soil condition and rooting variation. These factors influence the fractionation of tree-ring δ ¹⁸O variations can retain temperature or hydroclimate signals utilized by paleoclimatic studies (McCarroll and Loader, 2004; Nagavciuc et al., 2019; Yang et al., 2021). However, they are masked in the continental analysis by averaging the regional tree-ring isotopic values. From this analysis, we conclude that the low correlation between δ ¹⁸O and δ ²Hn at the local level is due to the temporal variations of isotope fractionations rather than database heterogeneity.

The Greenland sample is not used for the correlation analysis as it is not dated. Nevertheless, its proximity to the Conifers MWL (shown by the red dot in Figure 2B) could indicate that it experienced similar hydrological conditions and biological δ ²Hn fractionations. However, the sample has lower δ ²Hn and δ ¹⁸O values compared to the present precipitation isotopic values (Figure 3), which can be interpreted as being due to a lower temperature than nowadays (Supplementary Figures S4-1A, B), or that the tree used melting ice water that has lower δ ²H and δ ¹⁸O values.