Norway spruce shows stronger growth sensitivity and weaker intrinsic water-use efficiency response than Scots pine under increasing water limitation in southern Finland

Boreal forests, essential for carbon sequestration and multiple ecosystem services, are increasingly threatened by climate-induced water stress. This study investigates how increasing water limitation affects growth and intrinsic water-use efficiency (i_WUE) in Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies (L.) Karst) in southern Finland from 1990 to 2022. Tree-ring data on basal area increment (BAI) and carbon isotope discrimination (Δ¹³C) were collected from three sites per species, with Scots pine on mineral soils and Norway spruce on organic soils and under contrasting soil moisture conditions. The tree-ring records were combined with regional growth indices (GIs) from the Finnish National Forest Inventory (NFI) to enable broader inference across soil types. Our results indicate that water limitation intensified over the past decade, with pronounced growth declines after 2016 in Norway spruce, while Scots pine exhibited smaller reductions, evident at the site level from around 2009 but not resolved in the regional GI. Δ¹³C analyses indicated increased stomatal regulation in Scots pine and more limited physiological adjustment in Norway spruce after 2016, consistent with rising atmospheric moisture demand. Accordingly, i_WUE increased more strongly in Scots pine, whereas Norway spruce showed only modest changes despite a substantial growth decline. Interannual variability in growth and i_WUE was related both to vapour pressure deficit (VPD) and to the summer water balance, but their relative importance differed between species: VPD had the larger effect per standard deviation in Scots pine, whereas Norway spruce responded to atmospheric demand and to the summer water balance to a similar degree. In Norway spruce, growth after 2016 was significantly constrained by VPD, whereas Scots pine showed comparatively buffered growth responses and stronger increases in i_WUE. Soil-type effects were generally weak at the regional scale, but clearer patterns emerged along local moisture conditions. At the site level, both species tended to show greater negative post-change growth under drier conditions; however, Scots pine maintained higher mean growth on wetter sites, whereas Norway spruce showed higher mean growth under drier conditions. These contrasting growth–physiology responses indicate species-specific strategies in coping with increasing water limitation, with Norway spruce showing greater growth sensitivity and more limited physiological adjustment than Scots pine. Overall, boreal forest responses to drought appear to be structured primarily by species-specific physiological strategies, while being further modulated by local hydrological conditions. Accounting for both species differences and fine-scale variation in moisture availability will be critical for sustaining forest productivity and resilience under increasingly warm and water-limited climates in northern Europe, underscoring the need for site-adapted management that integrates species selection with local soil and moisture conditions.

1 Introduction

Boreal forests cover approximately one-third of the global forest area and play a crucial role in maintaining ecosystem services and acting as significant carbon sinks (Bonan, 2008; Bradshaw et al., 2009; Gauthier et al., 2015; Pan et al., 2011; Tagesson et al., 2020). However, rising temperatures and atmospheric water demand increasingly threaten forest health by reducing tree growth and weakening the forest carbon sink, increasing susceptibility to pests and diseases, and elevating mortality rates (Allen et al., 2010; Girardin et al., 2016; IPCC, 2023; Junttila et al., 2024; Linnakoski et al., 2017a, b; Martínez-García et al., 2024; Seidl et al., 2017; Zhang et al., 2024). Increasing water limitation has therefore emerged as an important driver of forest growth decline in the boreal region over recent decades (Gao et al., 2016; Gutierrez Lopez et al., 2021; Treydte et al., 2023). However, the effects of climate on species-specific growth and water-use regulation remain unclear.

Drought events have increased in frequency, duration, and intensity across Europe (Gutierrez Lopez et al., 2021; Knutzen et al., 2025; Lehtonen and Pirinen, 2019), with pronounced impacts emerging across the boreal region (Gao et al., 2016; Schär et al., 2004; Schuldt et al., 2020). These droughts, often coinciding with elevated air temperatures (T_a), increased vapour pressure deficit (VPD), and reduced soil moisture, are expected to become more frequent under continued warming (Gao et al., 2016; Swain et al., 2025). In southern Finland, characterized by cool and moist growing conditions, the recent rise in summer T_a has led to elevated atmospheric water demand, which amplifies vegetation water stress (Grossiord et al., 2020; Novick et al., 2024; Treydte et al., 2023). Although summer precipitation has shown a modest upward trend in parts of this boreal zone, these gains are minor relative to increasing evaporative demand (Grossiord et al., 2020; Irannezhad et al., 2014). Recent eddy-covariance evidence indicates that atmospheric drought has become more frequent across northern Europe and constrains ecosystem carbon uptake over a broad range of conditions, whereas strong soil-moisture limitation remains comparatively rare and spatially localized (Peltola et al., 2026) and the strongest reductions in carbon uptake occur when atmospheric and soil drought coincide.

In southern Finland, records from the 13th Finnish National Forest Inventory (NFI-13) have documented declines in tree growth, increased canopy dieback, and elevated mortality over the past decade, consistent with increasing water limitation. The two dominant tree species, Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies (L.) Karst.), which account for approximately 40 % and 38 % of the total stem volume, respectively, have both grown below expected trajectories, with Norway spruce exhibiting larger negative deviations (Mäkinen et al., 2022; Henttonen et al., 2024; Junttila et al., 2024; Korhonen et al., 2024; Fig. S1 in the Supplement). However, despite recent advances in identifying climatic drivers of growth decline (Laudon et al., 2024; Lee et al., 2026), important uncertainties remain regarding the extent to which species-specific physiological responses mediate these patterns. Addressing this gap requires integrated analyses of growth and physiology at both regional and site scales.

Species-specific traits shape forest responses to water limitation and warming, with consequences for forest composition and resilience (Anderegg et al., 2018). In boreal forests, Scots pine and Norway spruce exhibit distinct physiological and water-use strategies under varying water limitations, reflecting their respective strategies for coping with drought (Levesque et al., 2013; Oren et al., 1999). Norway spruce generally has a shallower and less plastic root system than Scots pine, increasing its sensitivity to short-term variation in water supply (Kalliokoski et al., 2008; Schmid and Kazda, 2002). Scots pine also generally has greater photosynthetic capacity (Thum et al., 2008) and stronger stomatal sensitivity to increasing VPD than Norway spruce (Lagergren and Lindroth, 2002). Under frequent high-VPD conditions, weaker stomatal down-regulation in spruce can lead to increased water loss, lower leaf water potentials, slower cell expansion, and heightened hydraulic risk when water supply is limited (Arend et al., 2021). These differences align with observed patterns of lower drought tolerance in spruce, often seen as more substantial growth reductions under water limitation (Motte et al., 2023; Ulrich and Grossiord, 2023). Resolving these species-specific sensitivities is critical for anticipating changes in boreal forest composition and guiding adaptive management.

Intrinsic water-use efficiency (i_WUE), defined as the ratio of net photosynthetic rate (A) to stomatal conductance (g_s), reflects a species' ability to balance carbon gain and water loss (Farquhar et al., 1989b; Seibt et al., 2008). Differences in i_WUE among species can signal variations in stomatal control, photosynthetic capacity, and water-use strategies, thereby affecting competitive interactions and influencing forest structure and dynamics (Ainsworth and Rogers, 2007; Anderegg et al., 2018). For instance, increases in i_WUE above the expected increase due to increasing atmospheric CO₂ generally indicate stronger stomatal regulation in response to water stress (Farquhar et al., 1982, 1989a; Kannenberg et al., 2021; Saurer et al., 2004; Szejner et al., 2018; Voelker et al., 2016). i_WUE is closely related to δ¹³C because photosynthetic carbon isotope discrimination (Δ¹³C) primarily reflects the ratio of internal-to-ambient CO₂ concentrations ($c_{\mathrm{i}}/c_{\mathrm{a}}$). Since A can be expressed as the product of g_s and the CO₂ gradient (c_a−c_i), i_WUE increases as $c_{\mathrm{i}}/c_{\mathrm{a}}$ decreases, leading to reduced Δ¹³C values. In some cases, tree species may exhibit long-term stable i_WUE, reflected by a consistent increase in Δ¹³C associated with a proportional increase in ($c_{\mathrm{i}}/c_{\mathrm{a}}$) under increasing c_a (Belmecheri et al., 2021; Guerrieri et al., 2019; Saurer et al., 2004; Voelker et al., 2016). Forest management such as thinning and fertilization can influence i_WUE by enhancing photosynthetic capacity and resource availability, often resulting in increased A and higher i_WUE (Brooks and Mitchell, 2011; Lehtonen et al., 2023a; Palvi et al., 2025). However, increasing i_WUE does not always correspond to increased tree growth or productivity, particularly in unmanaged forests where other limiting factors may prevail (Andreu-Hayles et al., 2011; Battipaglia et al., 2013; Kannenberg et al., 2022; Walker et al., 2020).

Tree responses to water limitation can differ between organic and mineral soils due to differences in water-holding capacity, aeration, and nutrient supply. Specifically, in organic peat soils, high water retention but low oxygen in the rooting zone constrains root function. Lowering of the water table (WT) through drainage can improve aeration and promote growth, but may increase drought exposure if the WT declines too deeply (Hökkä et al., 2025; Sikström and Hökkä, 2016). In peatlands, fine-root production is concentrated in the upper 20 cm, making trees particularly sensitive to WT declines during dry periods (He et al., 2023). This effect may be more significant in well-drained peatlands, where lower WT levels can more readily induce drought stress. In contrast, mineral soils typically drain faster and retain less water than organic soils, causing surface layers to dry rapidly during warm summers, though deeper rooting is possible where soil structure and oxygen availability permit (Heiskanen and Mäkitalo, 2002; Schenk and Jackson, 2002). Topography further modulates these patterns observed in mineral soils. In particular, lowlands tend to retain moisture and risk waterlogging, whereas uplands dry more quickly and face higher evaporative demand, creating microclimatic conditions that influence tree species performance (Kemppinen et al., 2023; McNichol et al., 2024). Understanding these tree-soil interactions is essential for assessing drought sensitivity and developing targeted forest management strategies.

Southern Finland is among the regions of the country most vulnerable to increasing water limitation (Lee et al., 2026; Lehtonen and Pirinen, 2019; Fig. 1). Here, we investigate how increasing water limitations influence growth and i_WUE of Scots pine and Norway spruce across mineral and organic soils under contrasting water-availability conditions in this region. To address this, we integrated site-level tree-ring analyses of basal area increment (BAI), Δ¹³C, and i_WUE with regional growth indices (GIs) derived from the Finnish NFI data spanning 1990–2022. We related these growth and physiological responses to T_a, VPD, and Standardized Precipitation-Evapotranspiration Index (SPEI; Vicente-Serrano et al., 2010) using statistical models. We hypothesize that: (i) Recent tree growth decline is primarily driven by increasing water limitation, with stronger effects in Norway spruce. (ii) Scots pine and Norway spruce exhibit divergent i_WUE responses to increasing atmospheric demand and water stress, reflecting species-specific strategies for coping with hydrological change. (iii) Tree growth and physiological responses differ between mineral and organic soils due to contrasting water availability, aeration, and nutrient dynamics.

Figure 1 Regional climate and multi-decadal trends in southern Finland. Data based on CRU TS data (version 4.08) for the period 1980–2023. Panels show (A) mean summer temperature (T_a) for June–August (JJA), (B) mean August standardized precipitation–evapotranspiration index (SPEI), (C) mean JJA vapor pressure deficit (VPD), (D–F) Mann-Kendall trends in JJA T_a (°C yr⁻¹), August SPEI (unitless yr⁻¹), and JJA VPD (kPa yr⁻¹), respectively. Circles in panels (A)–(C) show the location of the study sites, while the polygon outlines in panels (D)–(F) corresponds to the southern Finland region as defined by the Finnish NFI. Panels (G)–(I) show the regional averaged time series (1980–2023) of T_a, SPEI, and VPD, with 30-year spline curves and the shaded area shows two standard deviations around the period mean.

Our specific research questions were: (1) How has increasing water limitation affected regional growth trends of Scots pine and Norway spruce on mineral and organic soils in southern Finland over recent decades? (2) What are the species-specific physiological responses of Scots pine and Norway spruce to contrasting soil water availability conditions in this region? (3) How do T_a, VPD, and SPEI relate to growth (BAI) and physiological traits (Δ¹³C, i_WUE) in these species in southern Finland?

2 Methods

2.1 Study area

2.1.1 Land and forest cover

This study was conducted in Southern Finland (Fig. 1). This region includes 9 provinces (Uusimaa, Southwest Finland, Satakunta, Kanta-Häme, Pirkanmaa, Päijät-Häme, Kymenlaakso, South Karelia, and South Savo) covering an area of 75 222 km². According to data from the 13th Finnish National Forest Inventory (NFI-13, 2019–2023) (Korhonen et al., 2024), 53 073 km² (71 %) of the total region area is forestry land. This land is primarily privately owned and managed, with stands dominated by Scots pine (49 %) and Norway spruce (36 %). These tree species play a crucial role in both ecological functions and as a source of raw material for the forestry sector (Henttonen et al., 2024). Notably, on mineral soils, Scots pine dominated forests cover a larger area than those dominated by Norway spruce (20 558 versus 16 137 km²) (Korhonen et al., 2024). In contrast, on organic soils, Norway spruce dominated forests are more prevalent, covering 5298 km², compared to 4694 km² for Scots pine dominated forests. It is important to note that a significant proportion of the organic soils with Norway spruce and Scots pine have been drained (47 % and 85 %, respectively).

2.1.2 Regional climate

The 0.5° gridded Climatic Research Unit (CRU) Time-series (TS) data version 4.08 monthly data (Harris et al., 2020) for air temperature (T_a, °C), precipitation (P, mm), vapor pressure (VP_air, in kPa), and monthly average of potential evapotranspiration (PET, mm d⁻¹) over the period 1990–2023 was obtained using the Climate Explorer of the Royal Netherlands Meteorological Institute (KNMI) (Trouet and van Oldenborgh, 2013). Vapor pressure deficit (VPD, kPa) was estimated using Eq. (1).

$$\begin{array}{}\text{(1)}& \mathrm{VPD}=\left(\mathrm{0.6108}\times e^{\left(\left(\mathrm{17.27}\times T_{\mathrm{a}}\right)/\left(T_{\mathrm{a}}+\mathrm{237.3}\right)\right)}\right)-{\mathrm{VP}}_{\mathrm{air}}\end{array}$$

Where the expression in parentheses corresponds to saturation vapor pressure (VP_sat, kPa) estimated from T_a.

The PET was converted to mm per month to compute the water balance needed for the estimation of the SPEI. The water balance was computed as the difference between monthly precipitation and monthly PET. A 3-month SPEI was then computed from the water balance time series using the R package SPEI (Vicente-Serrano et al., 2010). The 3-month SPEI was selected because it captures short to medium term moisture variability relevant to tree physiological processes during the growing season. Data processing was performed in R software using the ncdf4 package for handling NetCDF data.

In this study, we selected T_a, VPD, and SPEI to represent thermal conditions, atmospheric water demand, and climatic water balance, respectively (Babst et al., 2019; Girardin et al., 2016). T_a was retained for descriptive and correlation analyses but excluded from the mixed models because of its high collinearity with VPD. The spatial distribution of T_a, SPEI, and VPD across southern Finland, along with their long-term temporal trends, reveals temperature increases of up to 0.05 °C yr⁻¹ in certain areas (Fig. 1D, G). Concurrently, the negative trend in SPEI indicates a decline in water availability (Fig. 1E, H). In addition, VPD exhibits a similar upward trend, reflecting increased atmospheric water demand (Fig. 1F, I).

2.2 Finnish National Forest Inventory data

GIs data were used to provide a broader regional context and validate site-specific observations. The GI describes the annual variation in tree growth that is not driven by time-dependent tree or stand factors, but by environmental variability, such as climate fluctuations (Henttonen, 2000). In practice, size- and age-related growth trends and other long-term developmental changes are removed by standardizing the tree growth time series. Consequently, fluctuations in the GI reflect year-to-year variation in temperature, precipitation, soil moisture, and other environmental factors rather than intrinsic growth trends associated with stand development. An annual GI value is typically expressed as a relative deviation from the long-term average of the annual GIs, which is set at a value of 100 (Henttonen, 2000).

In this study, we estimated GIs for the period 1990–2022 separately for Scots pine and Norway spruce growing on mineral and drained organic soils in Southern Finland. Specifically, the estimation of these GIs was based on sample plot data from five consecutive Finnish NFIs, covering the period from the 9th NFI (NFI-9, 1996–2003) to the 13th NFI (NFI-13, 2019–2023). The data selection focused on southern Finland, comprising the southernmost sampling area in NFI (Korhonen et al., 2021, 2024; Tomppo et al., 2011). Radial increments, measured from tree ring increment cores from NFI temporary sample plots, were used in the GI estimation. Tree-ring cores were collected at breast height (1.3 m above the ground), and the ring widths from these samples were measured in the laboratory. Data were limited to trees from the dominant crown layer. To avoid confounding effects of forest management, samples were restricted to sites with no felling in the past five years and no drainage operations on mineral soils within the past 30 years.

GIs were then estimated using mixed linear models (Henttonen, 2000), consistent with methods applied in the Finnish NFI (e.g., Korhonen et al., 2024). A detailed description of the statistical framework used for GI estimation, including the mixed-effects model formulation, uncertainty propagation, and breakpoint sensitivity analysis, is provided in Appendix A. The NFI-based GI dataset shows a gradual decline in sample size after approximately 2005, primarily reflecting changes in the design of the Finnish NFI, including the rotating sampling scheme and the availability of historical plot and tree-core data. Consequently, earlier time periods are represented by a larger cumulative number of observations than more recent years.

2.3 Tree-ring study sites

2.3.1 Site selection and plot categorization

Three study sites per tree species (Scots pine and Norway spruce) were established in southern Finland (Fig. 1) within representative forest stands. The Scots pine sites were identified as “Solböle”, “Vastamäki”, and “Ruotsinkylä”, while the Norway spruce sites were identified as “Janakkala”, “Ränskälänkorpi”, and “Lapinjärvi”. To validate the robustness and representativeness of our sampling design, we also tested tree growth correlations between the plot-level and regional NFI data (Fig. S2).

At each site, two plots (10 m radius) were established under contrasting water availability conditions, which were classified as “dry” and “wet”, resulting in a total of 12 plots across southern Finland. This moisture categorization was based on a combination of site-specific observations. Specifically, the Scots pine plots were in different topographical contexts on mineral soils with varying susceptibility to water limitation. Plots at each site were selected using the cartographic depth-to-water (DTW) index maps (2×2 m resolution) with 0.5 ha threshold (Murphy et al., 2007), where values below 100 indicate wet conditions and high water retention capacity. Thus, the “dry” and “wet” plots were situated on forest stands on upland and lowland soils, respectively. In the case of Norway spruce plots, they were in organic soils with distinct drainage systems. Plots at each site were selected according to the ditch spacing, with the narrower the distance between ditches, the deeper the water table level and the lower the soil moisture availability (Leppä et al., 2020a).

It should be noted that the selection of sites for Scots pine on mineral soil and Norway spruce on drained organic soil was based on the dominance of the species on these soil types in southern Finland (see Sect. 2.1).

2.3.2 Soil characteristics

Soil bulk density (BD, g cm⁻³) was measured in both Scots pine and Norway spruce plots at 0–10 and 10–20 cm depths on three randomly selected locations within the plot. In Scots pine plots, samples were collected using a stainless-steel cylinder (5.5 cm diameter, 4.0 cm height). In Norway spruce plots, samples were collected using a stainless-steel box sampler ($\mathrm{100}\times \mathrm{6.3}\times \mathrm{3.9}$ cm; length × width × height). All soil samples were dried to constant mass at 105 °C for 48 h, and BD was then determined by dividing the oven-dried soil mass by the original sample volume.

Soil samples were collected next to each BD sampling point at the same depths. In Scots pine plots, samples from the organic and mineral layers at the three locations were combined into a single composite sample per layer. In Norway spruce plots, all samples were pooled into a single composite sample per plot. Total carbon (C) and nitrogen (N) concentrations (%) were measured using a TruMac CN analyser (LECO Corporation, Saint-Joseph, MI, USA).

The volume of stones (V_s, %) was estimated in the 0–30 cm only in mineral soil layer at 32 sampling locations within each Scots pine plot, following 8 transects radiating from the plot centre along the cardinal directions with 4 points per transect. V_s was then assessed using a rod penetration method (Viro, 1952). Additionally, the depth of the organic layer (O_depth, cm) was measured at each of these sampling locations.

2.3.3 Climate at the site level

For each study site, we obtained the mean daily T_a and VP_air from the observational data interpolated to a 10×10 km grid covering Finland from the Finnish Meteorological Institute (FMI). The mean daily VPD was then estimated using Eq. (1), where VPD = VP_sat−VP_air. In addition, the 3-month SPEI was obtained from the 0.5° gridded regional SPEI dataset (see Sect. 2.1.2).

2.3.4 Tree increment core collection

Within each plot, we recorded the coordinates, diameter at breast height (DBH), total tree height (H_t, defined as the vertical distance from the base of the tree at ground level to the top of the highest living part of the crown), and tree crown base height (H_cb, defined as the vertical distance from the base of the tree at ground level to the lowest live branch) for every tree. The canopy composition was categorized by strata, quantifying the proportion of overstory, mid-story, and understory. For tree-ring analyses, 19–21 dominant and codominant trees per plot with DBH ≥ 20 cm were selected. When the required number of suitable trees was not available within the plot, additional individuals from a 5 m buffer zone surrounding the plot were included, excluding trees near drainage ditches in Norway spruce plots. From each selected tree, one 5 mm diameter increment core was extracted at DBH using an increment borer. Cores were oriented toward the pith to optimize analysis of annual tree-ring width and carbon isotope composition (δ¹³C). Trees were selected in apparently healthy condition: individuals showing crown dieback, stem damage, resin flow or bark-beetle entry holes were avoided at the time of coring. Damage status was not recorded systematically for each tree, so this was a field selection practice rather than a documented criterion.

2.3.5 Tree-ring preparation and analysis

Tree cores were prepared and analysed following standard dendrochronological procedures (Stokes and Smiley, 1996). Each tree core was immersed in water for at least 30 min to avoid breaking the cores, when mounted on wooden supports, and the core surface was prepared with a microtome (WSL Core Microtome; Swiss Federal Research Institute WSL, Switzerland) until the ring boundaries and wood cells were clearly visible under a binocular microscope. Then, each tree core was scanned at a resolution of 800 dpi using a colour scanner (Epson perfection v700 Photo). Finally, ring widths were measured with a precision of 0.001 mm using the WinDENDRO™ software (version Reg 2022a; Regent Instruments Inc., Quebec, Canada).

Cross-dating was performed to corroborate the calendar year assigned to each tree ring. Using standard techniques ensured the accurate assignment of calendar years to each growth ring (Stokes and Smiley, 1996). This process involved comparing ring width patterns among all samples to identify and correct any dating errors, thereby enhancing the reliability of the chronological series (Figs. S3–S14 and Tables S1–S12). All series were retained in the analyses, including those with lower correlations to the master chronology, as no evidence of cross-dating errors was detected; given the relatively short time series and the influence of natural forest dynamics, these cases are unlikely to affect the overall results.

After cross-dating, the ring-width time series were standardized and detrended to remove age-related growth trends and emphasize the environmental and climatic signals influencing growth variability (Cook, 1987). We utilized the R package dplR (Bunn, 2008, 2010) for data processing and analysis. Specifically, Regional Curve Standardization (RCS) was applied to the ring width data to construct site-specific chronologies that account for variations in growth related to tree age and size (Briffa and Melvin, 2011). This method involves creating a regional growth curve by averaging ring widths across all trees of the same cambial age, allowing us to preserve to some extent the low-frequency variability related to climatic factors and account for the age trends (Briffa and Melvin, 2011). Additionally, basal area increment (BAI, mm² per tree per year) was estimated using the annual diameter increment per tree, using the R package dplR (Bunn, 2008, 2010). Individual-tree BAI values were converted from mm² to m² and scaled using site-specific tree density from each plot, thereby expressing growth on an area basis (BAI, m² ha⁻¹ yr⁻¹).

2.3.6 Carbon isotope analysis (δ¹³C)

The δ¹³C analysis was performed on whole-ring samples to capture the integrated physiological signal of each annual ring. α-cellulose was extracted from the tree ring annual samples using standard chemical procedures adapted from Leavitt and Danzer (1993). From each plot, we selected five trees out of the 20 collected. The annual rings corresponding to each year were carefully sliced into small particles. To obtain sufficient and balanced masses for isotope analysis and reduce individual tree variability, the samples from the five trees were pooled for each annual ring, ensuring equal mass contribution from each tree (Belmecheri et al., 2022). Carbon isotope measurements were made on annual rings covering the period 1999–2022 at all 12 plots, giving 24 annual values per plot.

The pooled wood samples were placed on a Soxhlet with ethanol for 24 h to remove extractives. This step removes resins, oils, and other soluble substances that could interfere with the cellulose extraction. Following this, the samples underwent delignification by treatment with a sodium chlorite (NaClO₂) solution, buffered with acetic acid (CH₃COOH) to maintain a pH of around 4. This step was performed at 70 °C for 6 h to degrade lignin components and other secondary metabolites in the wood sample (Tarvo et al., 2010).

The delignified samples were then treated with a 17 % sodium hydroxide (NaOH) solution at room temperature for one hour to isolate α-cellulose by removing hemicelluloses (Green, 1963). After each chemical treatment, samples were thoroughly rinsed with deionized water until neutral pH was achieved. Finally, the purified α-cellulose was stabilized by treating it with a 10 % acetic acid solution to neutralize any remaining alkali, and then it was rinsed again with deionized water. The samples were freeze-dried and homogenized using an ultrasonic device prior to isotope analysis. α-cellulose was used in this study because it is a structural polymer closely derived from sugars produced during the growing season, and extracting cellulose allows for more effective homogenization when pooling samples from multiple trees, improving comparability and reliability (Belmecheri et al., 2022).

The δ¹³C of the α-cellulose samples were measured using an isotope ratio mass spectrometer (IRMS) at the Stable Isotope Laboratory of Luke (SILL), Natural Resources Institute Finland (Luke). Approximately 0.3 mg of each sample was weighed into tin capsules and combusted in an elemental analyser (EA) coupled to the IRMS. The results were normalized against IAEA-CH7 (−32.15 ‰) and in-house lactose and sucrose reference materials (−24.66 ‰ and −12.22 ‰, respectively), and are reported relative to the Vienna Pee Dee Belemnite (VPDB) standard. Repeat measurement of IAEA-CH3 material, run concurrently with the samples and normalized with the same scheme returns δ¹³C value of $-\mathrm{24.77}\mathrm{‰}\pm \mathrm{0.03}$ ‰ (mean ± standard deviation, n=18) in line with the certified value for this reference (−24.72 ‰).

2.3.7 Intrinsic water-use efficiency (i_WUE)

The δ¹³C can be used to estimate i_WUE at a variety of temporal resolution, from season (Tang et al., 2022) to multi-decadal and at spatial (leaf level to ecosystems) scales (Battipaglia et al., 2013; Gagen et al., 2011; Seibt et al., 2008). At the leaf level, the variability of δ¹³C in plant material represents the balance between net photosynthetic rate (A) and stomatal conductance to water vapor (g_s), which both processes discriminate against ¹³C. This isotopic effect and coupling is mainly influenced by environmental factors such as soil water availability, atmospheric water demand, photosynthetically active radiation, and atmospheric CO₂ concentrations (Belmecheri et al., 2021; Farquhar et al., 1982; Leavitt et al., 2011; Ma et al., 2023). In this study, we describe photosynthetic carbon discrimination (Δ¹³C) using a simplified mechanistic model that accounts for the primary factors affecting the discrimination against ¹³C during photosynthesis (Eqs. 2 and 3).

$$\begin{array}{}\text{(2)}& {\mathrm{\Delta }}^{\mathrm{13}}\mathrm{C}={\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{atmosphere}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{cellulose}}}{\mathrm{1}+{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{cellulose}/\mathrm{1000}}}}\end{array}$$

Where “δ¹³C_atmosphere” is the isotopic composition from the well mixed atmospheric CO₂ (Belmecheri and Lavergne, 2020) and “δ¹³C_cellulose” is the measured isotopic composition in cellulose from the tree rings, then Δ¹³C expression is:

$$\begin{array}{}\text{(3)}& {\mathrm{\Delta }}^{\mathrm{13}}\mathrm{C}=a+\left(b-a\right){\displaystyle \frac{c_{\mathrm{i}}}{c_{\mathrm{a}}}}\end{array}$$

Where “a” is the isotopic fractionation during CO₂ diffusion through stomata (∼4.4 ‰), “b” is the tissue-specific isotopic fractionation associated with net fractionation caused by Rubisco and PEPC (phosphoenolpyruvate carboxylase) and post-photosynthetic fractionation factors (Cernusak and Ubierna, 2022; Ubierna et al., 2022) integrated as (∼27 ‰), and “$c_{\mathrm{i}}/c_{\mathrm{a}}$” is the ratio of intercellular to ambient CO₂ concentrations (data from Belmecheri and Lavergne, 2020). The $c_{\mathrm{i}}/c_{\mathrm{a}}$ ratio is indicative of the relation between CO₂ supply to the leaves (via stomata) and CO₂ demand by photosynthesis in the chloroplast. At a given water evaporative demand determined by VPD, the $c_{\mathrm{i}}/c_{\mathrm{a}}$ ratio is linked to i_WUE as follows (Eqs. 4 and 5):

$$\begin{array}{}\text{(4)}& {\displaystyle }{c}_{\mathrm{i}}={\displaystyle \frac{{\mathrm{\Delta }}^{\mathrm{13}}\mathrm{C}-a}{b-a}}\times c_{\mathrm{a}}\text{(5)}& {\displaystyle }{i}_{\mathrm{WUE}}={\displaystyle \frac{A}{g_{\mathrm{s}}}}\approx {\displaystyle \frac{g_{{\mathrm{s}}_{\mathrm{c}}}\cdot (c_{\mathrm{a}}-c_{\mathrm{i}})}{{g_{\mathrm{s}}}_{\mathrm{w}}}}\approx {\displaystyle \frac{c_{\mathrm{a}}-c_{\mathrm{i}}}{\mathrm{1.6}}}\end{array}$$

Where “A” is the assimilation rate, $g_{{\mathrm{s}}_{\mathrm{c}}}\cdot (c_{\mathrm{a}}-c_{\mathrm{i}}$), and “$g_{{\mathrm{s}}_{\mathrm{c}}},{g_{\mathrm{s}}}_{\mathrm{w}}$” are the stomatal conductance of CO₂ and water (Ehleringer and Cerling, 1995; Farquhar et al., 1982). The i_WUE represents the efficiency of carbon gain per unit of potential water loss determined by g_s, and often varies with VPD, as higher VPD can lead to reduced stomatal conductance to minimize water loss (Medlyn et al., 2017).

2.4 Statistical methods

2.4.1 Post-hoc comparisons of forest stand and soil characteristics

Forest stand and soil characteristics were compared between “Dry” and “Wet” plots for each tree species. Pairwise comparisons for each variable were performed using Bonferroni-adjusted post-hoc tests to account for multiple comparisons, with significance set at p<0.05.

2.4.2 Detection of changes in the time series of regional- and site-level parameters

The time series derived from the NFI GI and site-level BAI, Δ¹³C, and i_WUE chronologies were modelled using segmented regression models (Muggeo, 2008). This segmented regression framework enables the estimation and detection of distinct temporal trends within different time series phases. Starting from an unsegmented linear trend model, Davies's test (Davies, 2002) was applied iteratively to evaluate the significance of additional structural breakpoints until no further significant breaks were detected. This approach was applied to each chronology, resulting in species- and variable-specific breakpoint estimates. All change-point analyses used the common 1999–2022 period defined by the isotope record, allowing direct comparison among GI, BAI, Δ¹³C, and i_WUE. Figure 2 presents the full GI record from 1990 to 2022.

Figure 2 National Forest Inventory (NFI) based tree growth patterns in southern Finland. Panels show NFI-based standardized growth index (GI) for Scots pine (A) and Norway spruce (B) from 1990 to 2022 on mineral (darkblue) and organic (lightblue) soils. The upper panels illustrate the annual variability in GI for each tree species and soil type, with thin lines depicting the mean annual values, while thick lines showing the 30-year smoothing splines to illustrate long-term trends. The lower panels show the number of trees used to derive each annual GI estimate for each tree species and soil type (see Appendix A for NFI uncertainty estimates).

Because breakpoints were not detected for every variable, common change years were used within each species to enable comparisons among soils and plot types: 2009 for Scots pine and 2016 for Norway spruce. Breakpoint results and their uncertainty are reported in Sect. 3.3 and Appendix A. The estimated marginal means treat the harmonised change years as fixed and should therefore be interpreted as conditional on those years. In contrast, breakpoint uncertainty was propagated through the segmented regressions and the mixed models in Eq. (6). Only Tables S13 and S14 use the NFI growth index and include NFI sampling uncertainty; Tables S15–S18 report site-level comparisons.

For each response and species, we fitted a model of the form response ∼ group × period × year, where group represented soil or plot type and year was centred within each period. Temporal slopes and estimated marginal means were compared among groups and periods using the emmeans package (Lenth and Piaskowski, 2017). Estimated marginal means are model-adjusted mean response values, whereas temporal slopes describe annual rates of change. Results are reported in Tables S13–S18. We distinguish between breakpoint-defined periods, based on species- and variable-specific change points estimated from segmented regression, and the common pre-2016/post-2016 periods used in linear mixed-effects models (LMMs) to facilitate direct comparisons of recent climate sensitivities between species. Accordingly, the terms “pre-change” and “post-change” refer specifically to breakpoint-defined periods, whereas “pre-2016” and “post-2016” refer to the common comparative temporal split.

2.4.3 Sensitivity of regional- and site-level parameters to environmental stressors

To evaluate the influence of the low and high frequency climate variability on inter-annual tree growth and physiology, we used both the detrended and non-detrended NFI-based GI, tree-ring variables (BAI and Δ¹³C), and the monthly climate indices (T_a, VPD, and SPEI). BAI and Δ¹³C were analysed because tree growth and physiological regulation of carbon gain relative to water loss may respond differently to climate variability and reflect different seasonal windows of sensitivity (Stockton Maxwell et al., 2020). Detrending was performed using the super smoother Friedman method implemented in the dplR package (Bunn, 2010). For each site and month, we then adjusted a simple linear model of the form of: (variable (detrended and not) ∼ Climate index (detrended and not)). We applied this simple linear model for all combinations of sites, species, and climate variables including previous fall and extracted the coefficient of determination (R²), slopes, and p-values from each model run.

Additionally, we used linear mixed-effects models (LMMs, Eq. 6) to examine the sensitivity of tree growth and physiological responses to SPEI and VPD. T_a was excluded from the models due to its high collinearity with VPD. The analysis was conducted on BAI and Δ¹³C chronologies, with monthly SPEI and VPD averaged across the May-June-July (MJJ) period, which was chosen because it best captures the climatic conditions influencing the physiological processes of the study species. Climate predictors and years were centred before modelling to facilitate the interpretation of changes in intercepts and differences in predictor effects, which correspond to interaction effects across periods. For comparative analyses, a common breakpoint year (2016), identified from the segmented regression analysis (Sect. 2.4.2), was used to define pre-2016 (1999–2015) and post-2016 (2016–2022) periods for both species. This categorical period variable was included in the LMMs to test whether climate sensitivity to SPEI and VPD differed between periods, potentially reflecting intensified water limitation after 2016. Before model fitting, climate predictors and year were mean-centred by subtracting their respective means. To enable comparison between SPEI and VPD, slopes were rescaled to represent a one-standard-deviation increase in each predictor. Predictor standard deviations were calculated separately for each species over the full analysis period using plot-years with complete SPEI and VPD data (n=144 per species). This linear rescaling does not alter the fitted relationships, test statistics, or p-values. Rescaled slopes are presented in Fig. 6 and Table S20. Models were fitted separately for Scots pine and Norway spruce using the following structure:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{y}_i& ={\mathit{\alpha }}_{\mathrm{period}\left(i\right)}+{\mathit{\beta }}_{\mathrm{period}\left(i\right)}\cdot {\mathrm{SPEI}}_i+{\mathit{\gamma }}_{\mathrm{period}\left(i\right)}\cdot {\mathrm{VPD}}_i\\ & +e_{\mathrm{site}\left(i\right)}+{\mathit{ϵ}}_i\end{array}\end{array}$$

Where y_i represents either BAI or Δ¹³C of the ith observation, e_site is a site-specific random intercept, and ϵ_i is a random residual error, which is assumed uncorrelated between observations. Because annual observations within a plot are temporally dependent, we modelled within-plot residuals with a first-order autoregressive process, AR (1), indexed by year. This allows residuals closer in time to be more correlated and improves precision of fixed-effect estimates under temporal dependence. Models were fitted using restricted maximum likelihood (REML). Parameters α_period, β_period, and γ_period, for period∈{pre-2016,post-2016}, as well as the between-site variances of random effect e and the AR(1) parameters were estimated using the lme function of the R package nlme (Pinheiro et al., 1999).

Because the split year in Eq. (6) was estimated rather than known, treating it as fixed would underestimate uncertainty in the period contrasts. We therefore propagated uncertainty in the Norway spruce breakpoint through the mixed models using the 1000 GI realizations generated for each soil type from the NFI uncertainty distribution in Appendix A (2000 draws in total). For each realization, the breakpoint nearest 2016 was used to define the pre- and post-change periods, and Eq. (6) was refitted separately for BAI and Δ¹³C. The fixed-effect estimates were then pooled across realizations using Rubin's rules (Rubin, 1987) with the small-sample degrees-of-freedom adjustment of Barnard and Rubin (1999). When no breakpoint was detected, the corresponding no-change model, BAI or Δ¹³C ∼ cSPEI + cVPD, was fitted with the same random-intercept and AR(1) structure. These realizations therefore contributed no period change to the pooled estimates, which represent model averages across breakpoint and no-breakpoint cases rather than estimates conditional on detecting a breakpoint. Breakpoints that would leave fewer than five years in either period were restricted to the nearest feasible year. We report the pooled estimates, the proportion of GI realizations in which a breakpoint was detected, and the proportion in which each model term was significant. The pooled p-value incorporates uncertainty across realizations, whereas the latter proportion indicates the stability of the inference to breakpoint uncertainty. Because the breakpoint year varies among realizations, the pre- and post-change periods also vary slightly, and the pooled variance therefore includes uncertainty in period definition as well as in the model coefficients.

To aid ecological interpretation, regression slopes were scaled relative to the observed variability of each response variable. BAI ranged from 0.14 to 1.28 m² ha⁻¹ yr⁻¹ (range = 1.14), and Δ¹³C from 16.5 ‰ to 22.3 ‰ (range = 5.8 ‰). Effect sizes were classified consistently for both variables as small (<10 % of the observed range), moderate (10 %–50 %), and large (>50 %), corresponding to slope thresholds of <0.11, 0.11–0.57, and >0.57 m² ha⁻¹ yr⁻¹ for BAI, and <0.58, 0.58–2.90, and >2.90 ‰ for Δ¹³C. These classes were used for interpretive purposes only and do not replace statistical inference.

3 Results

3.1 NFI-based tree growth patterns in southern Finland

Our analysis of NFI-based GI data for Scots pine and Norway spruce on both mineral and organic soils over 1990–2022 provides a regional overview of tree growth variability across southern Finland (Fig. 2). For Scots pine, GI values were broadly similar between soil types during the pre-change period (1990–2008), as defined by breakpoint analysis (see Sect. 3.3.). There were no significant differences between mineral and organic soils in either temporal trends (Δslope = 0.61±0.56 yr⁻¹, p=0.285) or mean growth levels (ΔEMM = $-\mathrm{3.0}\pm \mathrm{3.1}$, p=0.292) (Tables S13–S14). During the breakpoint-defined post-change period (2009–2022), GI declined on both soil types, with slightly stronger negative trends on mineral soils ($-\mathrm{1.46}\pm \mathrm{0.62}$ yr⁻¹) than on organic soils ($-\mathrm{1.06}\pm \mathrm{0.62}$ yr⁻¹), although this difference was not statistically significant (Δslope = $-\mathrm{0.40}\pm \mathrm{0.87}$ yr⁻¹, p=0.654) (Tables S13–S14). Mean GI was slightly higher on mineral soils during this period, but the difference was weak and not statistically supported (ΔEMM = 6.0±3.5, p=0.105), indicating limited divergence between soil types in both growth trajectories and average growth levels.

Propagating NFI sampling uncertainty through breakpoint detection establishes that no regional GI breakpoint is supported for Scots pine on either soil type, a breakpoint being detected in 0.9 % of 1000 GI realisations on mineral soils and 32 % on organic soils, and near 2019 rather than 2009 in the latter case (Appendix A). The 2009 partition used here is the breakpoint detected in site-level BAI at three of the six plots (Table S19), applied to the regional index as a common comparison window.

Norway spruce showed positive and similar GI trends across soil types during the breakpoint-defined pre-change period (1990–2015), with no evidence of soil-specific differences in either temporal dynamics (Δslope = $-\mathrm{0.17}\pm \mathrm{0.44}$ yr⁻¹, p=0.704) or mean growth levels (ΔEMM = 2.2±3.3, p=0.502) (Tables S13–S14). In the breakpoint-defined post-change period (2016–2022), GI declined sharply on both soil types, with more negative slopes on mineral soils ($-\mathrm{9.71}\pm \mathrm{2.20}$ yr⁻¹) than on organic soils ($-\mathrm{5.87}\pm \mathrm{2.20}$ yr⁻¹), although this difference was not statistically significant (Δslope = $-\mathrm{3.84}\pm \mathrm{3.11}$ yr⁻¹, p=0.245). Mean GI was somewhat higher on organic soils during this period, but soil differences remained statistically unsupported (ΔEMM = $-\mathrm{4.0}\pm \mathrm{6.2}$, p=0.527) (Tables S13–S14). This regional decline is well resolved by the NFI sample on mineral soils, where a breakpoint is detected in every one of 1000 GI realisations and the pooled post-breakpoint slope remains significantly negative; on organic soils the pooled slopes are of the same sign and magnitude but the sampling uncertainty is large enough that they are no longer distinguishable from zero (Appendix A, Table A2).

Overall, both species exhibited a marked reduction in GI following the breakpoint-defined post-change period, whereas soil-specific differences in both temporal trends and mean levels remained weak and statistically non-significant. Decreasing sample sizes toward the end of the study period (Fig. 2) may have increased uncertainty in late-period GI estimates, particularly because this reduction overlapped with the observed GI decline and could affect estimates of recent growth variation. Once this sampling uncertainty is propagated (Appendix A), the post-change decline in Norway spruce GI remains statistically supported on mineral soils and retains its sign and magnitude on organic soils, and in both cases its direction agrees with the independent site-level BAI analyses (see Sect. 3.2.).

3.2 Site-level patterns in tree growth and i_WUE

Sample plots in Scots pine sites on mineral soils exhibited contrasting topographical conditions. Upland dry plots showed significantly higher depth-to-water index (257±78 cm) and a greater proportion of rocks (55.3 %±3.0 %) than lowland wet plots (36±4 cm and 10.7 %±6.1 %, respectively). In contrast, plots in Norway spruce sites on organic soils were characterized by distinct drainage systems, with wet plots showing significantly wider ditch spacing (65±16 m) than dry plots (46±4 m). Despite these topographical and hydrological differences, plot-level measurements indicated no significant variation in forest stand or soil characteristics between dry and wet plots within each species (Tables 1 and 2).

Table 1 Forest stand characteristics for Scots pine and Norway spruce sites under “Dry” and “Wet” conditions. T_d: tree density, DBH: diameter at breast height, H_t: total tree height, H_bc: tree crown base height, BA: basal area, Age: mean stand age. Values are mean ± standard deviation. For each tree species and variable, a Bonferroni-adjusted post hoc test was applied to compare differences in means between plots (p<0.05). No significant differences were identified between the “Dry” and “Wet” plots for any variable across tree species. n=3 sites per tree species and plot type.

Table 2 Soil characteristics at Scots pine and Norway spruce sites under “Dry” and “Wet” conditions. BD: soil bulk density, C: carbon concentration, N: nitrogen concentration, $\mathrm{C}/\mathrm{N}$: $\mathrm{carbon}/\mathrm{nitrogen}$ ratio, O_depth: topsoil organic layer depth, V_s: volume of stones up to 30 cm depth. Values are mean ± standard deviation. Layer references: BD, C, N, $\mathrm{C}/\mathrm{N}$ correspond to topsoil or mineral layers as appropriate. O_depth and V_s are reported for Scots pine only; soils are organic with zero stones for Norway spruce. For each tree species and variable, a Bonferroni-adjusted post hoc test was applied to compare differences in means between plots (p<0.05). V_s differed significantly for Scots pine, while other variables showed no significant differences. n=3 sites per tree species and plot type.

^a,b Different superscript letters indicate a significant difference in mean stone volume between wet and dry Scots pine plots (Bonferroni-adjusted post hoc test, p<0.05).

From the breakpoint-defined post-change period onward (2009–2022 for Scots pine; 2016–2022 for Norway spruce), BAI showed contrasting trajectories between species and soil moisture conditions (Fig. 3A–B). In Scots pine, BAI exhibited a modest decline under dry conditions ($-\mathrm{0.0082}\pm \mathrm{0.0034}$ yr⁻¹) and near-stable dynamics under wet conditions ($-\mathrm{0.0003}\pm \mathrm{0.0034}$ yr⁻¹; Table S15), with consistently lower mean growth in dry compared to wet conditions (ΔEMM = −0.22 m² ha⁻¹ yr⁻¹, p<0.001; Table S16). In contrast, Norway spruce showed a pronounced breakpoint-defined post-change reduction in BAI under both moisture regimes (dry: $-\mathrm{0.069}\pm \mathrm{0.015}$ yr⁻¹; wet: $-\mathrm{0.061}\pm \mathrm{0.015}$ yr⁻¹; Table S15), accompanied by a significant shift in mean growth patterns, with higher BAI in dry than wet conditions (ΔEMM = 0.211±0.047, p=0.0007; Table S16), indicating a decoupling between temporal decline and average growth levels.

Figure 3 Trends in basal area increment (BAI) and intrinsic water use efficiency (i_WUE) for Scots pine and Norway spruce under “Dry” and “Wet” conditions during 1990–2022. Panels (A)–(B) show values for BAI, while panels (C)–(D) show values for i_WUE. Solid and orange lines display dry conditions, while dashed and blue lines represent wet conditions. Shaded areas indicate the 90 % confidence intervals around the smoothed lines using a 30-year smoothing spline.

Responses of i_WUE contrasted between species responses over the same breakpoint-defined post-change periods (Fig. 3C–D). Specifically, Scots pine exhibited a stronger increase in i_WUE under dry conditions (2.160±1.510 yr⁻¹), exceeding that under wet conditions (1.400±1.510 yr⁻¹; Table S17), with consistently higher mean i_WUE in dry relative to wet conditions (ΔEMM = 11.3±4.9 µmol mol⁻¹, p=0.039; Table S18). This pattern suggests enhanced physiological regulation of water use under increasing atmospheric demand and moisture limitation. In contrast, Norway spruce showed more moderate and less consistent i_WUE responses, with non-significant differences in temporal trends between dry and wet conditions (0.889±0.509 and 1.286±0.509 yr⁻¹, respectively; Table S17). Mean i_WUE values were slightly lower under dry conditions (69.6±1.2 µmol mol⁻¹) than under wet conditions (72.6±1.2 µmol mol⁻¹), although this difference was not statistically significant (ΔEMM = $-\mathrm{3.0}\pm \mathrm{1.7}$ µmol mol⁻¹, p=0.090; Table S18).

Overall, these results indicate contrasting post-change growth–physiology relationships between species. Norway spruce experienced a pronounced reduction in BAI under both moisture conditions, while showing only weak and variable changes in i_WUE. In contrast, Scots pine maintained relatively stable growth under wet conditions but exhibited increased i_WUE under dry conditions, suggesting a stronger physiological adjustment despite comparatively limited changes in growth. Temporal variability among individual sites is reflected in the spread of site-level trends (Figs. S15–S16), but the overall patterns remained consistent across the dataset.

3.3 Change-points in tree growth and i_WUE over the study period

Segmented regression analysis revealed species-specific differences in the timing of changes in growth and physiological trends. (Fig. 4; Table S19). In Norway spruce, significant breakpoints were detected in all four variables, although not at every plot: in GI on both soil types, BAI at five of six plots, Δ¹³C at three of six plots, and i_WUE at four of six plots. Most breakpoints occurred between 2015 and 2017, indicating a broadly coherent shift in growth and physiological responses. In contrast, Scots pine showed significant breakpoints only in BAI, at three of six plots in 2007, 2007, and 2011. No consistent breakpoints were detected for GI, Δ¹³C, or i_WUE. For GI in both species, reduced sample sizes in recent years may increase uncertainty in breakpoint timing and magnitude (Appendix A).

Figure 4 Change-points in tree growth and water use over the study period 1999-2022. Panels show the time series for (A–B) standardized Growth Index (GI) from the Finnish National Forest Inventory (NFI), (C–D) basal area increment (BAI), (E–F), carbon isotope discrimination (Δ¹³C), and (G–H) intrinsic water-use efficiency (i_WUE) for Scots pine (left panels) and Norway spruce (right panels). Within each panel, the grey lines represent time series in which no statistically significant detected change points were detected. In contrast, coloured lines indicate time series with one or more detected change points, while dashed lines show the corresponding segmented regression model fits. Below each time-series panel, density plots illustrate the temporal distribution of estimated change points across all time-series as probability density function (PDFs). The y-axis differs among plots because PDFs are normalized, whereas peak density values reflect the clustering of detected change points for each variable and species.

In Norway spruce, the breakpoints generally marked shifts from positive or stable trends to declining GI, BAI, and Δ¹³C (Fig. 4B, D, F; Table S19). The strongest and most consistent response was observed in BAI, indicating a substantial reduction in growth after 2016. Concurrent declines in GI and Δ¹³C further suggest that growth reductions were accompanied by shifts in physiological response.

In Scots pine, BAI breakpoints represented moderate reductions in growth rates rather than abrupt declines (Fig. 4C; Table S19). No breakpoint was detected in the observed GI series on either soil type, and the simulation analysis provided no support for a regional GI breakpoint. Thus the 2009 reference year represents a comparison point derived from site-level BAI rather than a regional change point.

Accordingly, 2016 was used for Norway spruce and 2009 for Scots pine in the within-species comparisons among soils and plot types reported in Tables S13–S18. For the subsequent LMM analysis of climate sensitivity, 2016 was used for both species to characterise post-change dynamics during 2016–2022 and enable interspecific comparison. The earlier Scots pine BAI breakpoint was not used in these models because it was not supported by the other response variables. For Norway spruce, breakpoint uncertainty was propagated by redrawing the split year from its empirical distribution in each realisation. The resulting estimates in Table 3 therefore incorporate uncertainty in breakpoint timing.

3.4 Environmental influence on tree growth and carbon isotope discrimination

Linear regression analysis (Fig. 5) illustrates how 3-month seasonal means of key climate variables, T_a, SPEI, and VPD, relate to tree-ring characteristics and NFI-based GIs. Among these, water availability (SPEI) and atmospheric demand (VPD) emerged as the strongest influential drivers. T_a showed similar temporal patterns to VPD, suggesting its effects are mediated through increased evaporative demand.

Figure 5 Heatmaps showing the strength and significance of relationships between climate variables and tree growth and physiological indicators across different sites and seasonal windows (3-month periods). Panels illustrate the relationships between T_a, SPEI, and VPD with Δ¹³C (A–C), BAI (D–F), and regionally NFI-based GI (G–I). In all panels, the x-axes represent the 3-month periods from September of previous year to November to the current growth year, while the y-axes show site-level chronologies for Scots pine and Norway spruce (A–F), as well as regional chronologies (G–I). Colour shading represents the coefficient of determination (R²) from simple linear regressions, ranging from white (low R²) to red (high R²). Positive (+) and negative (−) signs indicate the direction and statistical significance of the relationship (p<0.05). All analyses shown in this figure were conducted on non-detrended data.

At the site level, Δ¹³C showed the strongest relationships with climatic variables during July–August, when low SPEI and high VPD coincided with reduced carbon isotope discrimination (Fig. 5A–C). These patterns suggest that intensified summer water stress induces greater stomatal closure, particularly in Scots pine, leading to lower Δ¹³C values. T_a exerted a weaker but occasionally notable influence on Δ¹³C, due to its indirect role in increasing VPD. This indicates that high summer T_a can further exacerbate water stress, most noticeably in Scots pine sites. In contrast, BAI but especially NFI-based GI exhibited weaker correlations with moisture and temperature variables compared to Δ¹³C (Fig. 5D–I).

BAI in Scots pine shows minimal changes in intercepts across periods (Table 3). Growth sensitivities remain close to zero for SPEI (small effects) and are negative for VPD (moderate effects; Fig. 6A, Table 4). In contrast, Δ13C shows clearer and more consistent climate signals, with positive responses to SPEI and negative responses to VPD (Fig. 6B), with no significant differences between periods (Table 3). In terms of effect size, VPD exerts large negative effects on Δ13C, while SPEI effects remain small (Table 4). Together, these results indicate that although growth responses in Scots pine remain limited, the isotopic sensitivity to moisture stress has persisted in the post-2016 period.

Figure 6 Sensitivity of basal area increment (BAI) and carbon isotope discrimination (Δ¹³C) to recent changes in water limitation. Regression slopes and 95 % confidence intervals for the effects of mean May–July SPEI-3 and VPD on annual tree-ring BAI (A, C) and Δ¹³C (B, D) in Scots pine and Norway spruce, respectively, based on the LMMs described in Eq. (6). Because models were fitted separately for each species, slopes were rescaled to represent a one-standard-deviation increase in each predictor within the corresponding species-specific dataset. The predictor standard deviations were 0.759 for MJJ SPEI-3 and 0.100 kPa for MJJ VPD in Scots pine, and 0.741 and 0.089 kPa, respectively, in Norway spruce. Slopes therefore indicate the estimated change in BAI (m² ha⁻¹ yr⁻¹) or Δ¹³C (‰) associated with a one-standard-deviation increase in the corresponding climate predictor. Blue and orange points represent the pre-2016 (1999–2015) and post-2016 (2016–2022) periods, respectively. Error bars denote 95 % confidence intervals. Numerical estimates are reported in Table S20.

In Norway spruce, the baseline level of Δ¹³C differed between periods (Δα=0.455, SE 0.146, p=0.002), whereas the baseline level of BAI did not (Δα=0.033, SE 0.051, p=0.520; Fig. 6, Table 3). VPD emerges as a clear constraint on BAI after 2016 (Fig. 6C): the post-2016 slope is negative, moderate in effect size, and significantly different from zero (Table 4), whereas the pre-2016 slope is not. The change in slope between the two periods is not itself statistically supported (Table 3), so the evidence is for a significant VPD constraint on growth after 2016 rather than for a measured increase in the strength of that constraint (Table 3 and 4). The SPEI effect on Norway spruce BAI is negative in both periods and becomes stronger after 2016 (Fig. 6), although the difference between periods is only marginal (Fig. 6C, Table 3 and 4). The corresponding effect on Δ¹³C is positive in both periods. For Δ¹³C, sensitivities to both VPD and SPEI show no significant differences between periods; in effect-size terms, VPD responses shift from moderate pre-2016 to small post-2016, while SPEI effects remain small (Table 4).

Under the May-July climate window, climate sensitivity of Norway spruce BAI becomes more pronounced after 2016, indicating enhanced growth limitation. In contrast, Scots pine Δ13C shows more prominent and consistent climate sensitivity than Norway spruce, with large responses to VPD and small responses to SPEI, while Scots pine BAI remains comparatively buffered from climatic variability (Table 4).

Propagating uncertainty in breakpoint timing through Eq. (6) (Sect. 2.4.3) preserved the signs and significance status at the 0.05 threshold of the pre-period SPEI and VPD coefficients in both species and for both response variables (Table S21). The between-period slope differences also retained their signs but remained non-significant. Although a breakpoint was detected in approximately 95 % of the pooled Norway spruce GI realisations (Appendix A), this does not establish a change in climate sensitivity. For Norway spruce, the BAI intercept shift changed sign and remained non-significant, while the Δ13C intercept shift was no longer significant when realisations without a detected breakpoint were included (p=0.079). These results support climate associations, while the period differences in SPEI and VPD slopes and the Norway spruce intercept shifts in Eq. (6) remain statistically uncertain.

Table 3 Pre- and post-2016 changes in the sensitivity of the basal area increment (BAI) and carbon isotope discrimination (Δ¹³C) to climate variables for Scots pine and Norway spruce. Changes in model-estimated parameters for the Standardized Precipitation-Evapotranspiration Index (SPEI, β) and vapor pressure deficit (VPD, γ) for the period May–July (MJJ) based on model runs described in Eq. (6) (p<0.001 ”${}^{***}$”, 0.01 ”${}^{**}$”, 0.05,”^*”, 0.1 ”.”).

Table 4 Climate sensitivity of basal area increment (BAI) and carbon isotope discrimination (Δ¹³C) for Scots pine and Norway spruce before and after 2016. Slopes are derived from linear mixed-effects models for atmospheric vapor pressure deficit (VPD) and the Standardized Precipitation–Evapotranspiration Index (SPEI). Effect-size classes (small, moderate, large) are defined relative to observed variability (BAI: 0.14–1.28 m² ha⁻¹ yr⁻¹; Δ¹³C: 16.5 ‰–22.3 ‰) and represent the proportional contribution of each slope to the total observed range of variability per unit climate forcing. BAI slope = m² ha⁻¹ yr⁻¹ per unit and Δ¹³C slope = ‰ per unit.

4 Discussion

In this study, we combined regional GI estimates from the NFI with site-level tree-ring analyses of BAI and Δ¹³C to evaluate how Scots pine and Norway spruce respond to intensifying water stress across mineral and organic soils under contrasting water-availability conditions in southern Finland. Both species showed growth reductions during their breakpoint-defined post-change periods, but the magnitude and physiological response differed. Norway spruce showed particularly strong post-change declines in GI and BAI across both soil types, whereas Scots pine exhibited comparatively smaller growth reductions and maintained relatively stable growth under wet conditions. In contrast, Scots pine displayed a stronger increase in i_WUE under recent dry conditions, while Norway spruce showed weaker and less consistent i_WUE responses despite substantial growth decline. These contrasting growth-physiology relationships suggest species-specific differences in adjustment to ongoing shifts in increasing atmospheric moisture demand (Fig. 1), with important implications for future boreal forest dynamics, particularly in drought-prone environments such as southern Finland. Rising atmospheric dryness can induce considerable stress on boreal trees by increasing evaporative demand and reducing water availability, thereby constraining carbon uptake and growth (Adams et al., 2017; Frank et al., 2015). Consequently, species-specific differences in physiological regulation and growth sensitivity are likely to influence future forest dynamics, competitive interactions, and adaptive forest management strategies under ongoing climate change (Levesque et al., 2013).

After comparing the NFI data between soil types, the mineral and organic soils show similar temporal trends in GI (Fig. 2). Site-level analyses revealed additional contrasts associated with local moisture conditions, particularly in mean BAI and i_WUE (Fig. 3), although these differences were not always statistically robust. Scots pine tended to exhibit stronger growth reductions on dry upland sites, while also showing increased i_WUE, suggesting enhanced physiological adjustment under increasing water limitation. In contrast, Norway spruce showed consistent but non-significant tendencies toward stronger growth reductions on dry sites and generally higher mean BAI on dry rather than wet poorly drained sites, indicating reduced growth under wetter conditions. Its i_WUE responses were weaker and differences across moisture conditions not statistically significant. These patterns indicate that local hydrological conditions may exert a stronger control on growth and physiological responses than broad mineral-versus-organic soil classifications alone.

Together, the results highlight that future forest vulnerability under climate change will likely depend on the interaction between species-specific drought sensitivity and spatial heterogeneity in soil water availability across the boreal landscape in southern Finland.

4.1 Environmental drivers on recent growth declines

The NFI data show a clear decline in GIs for both Scots pine and Norway spruce over the most recent decade (Henttonen et al., 2024; Fig. 2, Appendix A) consistent with our independent plot-based growth records (Fig. S2). Declines across species and soil types point to environmental drivers operating beyond individual sites. The period contrasts indicate changes in slopes rather than baseline levels. After propagating uncertainty in the 2016 breakpoint through the mixed models, slope contrasts retained their direction and significance class, whereas changes in the Norway spruce intercepts were not significant. The results therefore support altered climate sensitivity after 2016, but not an abrupt shift in baseline growth or Δ¹³C. In Norway spruce, SPEI had contrasting relationships with the two responses. Δ¹³C increased with SPEI in both periods, consistent with reduced stomatal limitation in wetter summers, whereas BAI decreased with SPEI, particularly after 2016. The negative BAI–SPEI relationship therefore cannot be interpreted as evidence of soil water limitation. All Norway spruce plots were on drained organic soils, where growth is also influenced by drainage and water-table dynamics that were not measured here. Evidence for an atmospheric drought constraint is clearer for VPD, which was negatively associated with growth after 2016.

Biotic damage may also have contributed to the Norway spruce decline. Drought can increase Norway spruce susceptibility to the European spruce bark beetle (Ips typographus) and associated fungi (Netherer et al., 2024), while heat, drought, and windthrow can promote outbreaks (Netherer et al., 2021). Visibly damaged trees were avoided during sampling, but damage status was not recorded, and the tree-ring data include only trees surviving until sampling. Nevertheless, a similar post-2016 shift in the independent regional NFI growth index and the concurrent change in Δ¹³C indicate that undetected biotic damage alone is unlikely to explain the observed pattern. However, a compounding effect of bark beetles or pathogens cannot be excluded, and survivorship bias may cause stand-level growth losses to be underestimated.

According to our observations, increased VPD is a significant driver of the recent growth decline observed in Norway spruce in southern Finland. Based on site-level data, Scots pine showed a slight decrease in growth in dry conditions during the last decade. However, in both dry and wet conditions, growth variability was accompanied by a steep increase in i_WUE (Fig. 3). The large Δ¹³C–VPD slopes observed for Scots pine ($\approx -\mathrm{4.5}$ to −4.9 ‰ per unit VPD) indicate a strong stomatal response to atmospheric drought. However, this pronounced isotopic sensitivity, growth reductions in Scots pine remained moderate. This mismatch in sensitivities suggests that Scots pine was able to regulate water loss under high evaporative demand without a proportional reduction in radial carbon accumulation.

In contrast, Norway spruce exhibited a substantial decline in growth in both dry and wet conditions, especially after 2016, alongside relatively stable i_WUE. This different mismatch suggests that growth reductions were not tightly coupled to increased stomatal regulation. Instead, the decline may reflect non-stomatal limitations to carbon accumulation, potentially associated with hydraulic stress, reduced photosynthetic capacity, or altered carbon allocation. This interpretation is consistent with the lower stomatal sensitivity to atmospheric drought reported for Norway spruce compared to pine species (Ewers et al., 2001; Oren et al., 1999).

Correlation analyses indicate that Norway spruce growth is not only sensitive to current-year VPD but also shows stronger lagged and seasonal responses to VPD, particularly in previous-year autumn conditions. This response is consistent with documented carryover effects in spruce growth (Bouriaud and Popa, 2009), and whereby structural carbohydrate dynamics and storage of assimilates link current growth to climatic conditions in previous years. In particular, autumn carbon allocation and subsequent spring depletion of stored reserves can mediate delayed responses to climatic stress such as elevated VPD. Consistently, growth responses in this species became more strongly constrained after 2016, while Δ¹³C–climate relationships did not show comparably strong shifts between periods, indicating a divergence between physiological response and growth limitation. This pattern points to an increasing role of cumulative stress effects or non-stomatal limitations (e.g., carbon allocation or hydraulic constraints) in controlling growth.

Declining SPEI and increasing VPD intensify physiological stress by exerting additional pressure from drier soils and a heightened atmospheric moisture demand (Allen et al., 2010, 2015; Williams et al., 2013). This process accelerates soil moisture depletion and increases drought severity (Novick et al., 2016). When expressed per standard deviation of each predictor (Table S20), SPEI and VPD showed contrasting importance between species. VPD had the larger effect across both responses and periods in Scots pine, whereas their effects were more similar in Norway spruce, with SPEI larger in three of the four comparisons. This suggests that Scots pine responses were more strongly associated with atmospheric dryness, while Norway spruce responded similarly to atmospheric demand and summer water balance. The impact of these processes is evident in the responses documented in Scots pine and Norway spruce (Fig. 6). Norway spruce showed a notable decline in growth, particularly after 2016, and was more severely impacted by rising VPD, likely due to its physiological characteristics and adaptation to colder, more humid conditions (Lagergren and Lindroth, 2002; Levesque et al., 2013; Pichler and Oberhuber, 2007; Thum et al., 2008).

Overall, our results align with broader trends observed not only in boreal forests across Fennoscandia (Henttonen et al., 2024; Junttila et al., 2024; Laudon et al., 2024), but also throughout Europe (Knutzen et al., 2025; Römer et al., 2025), where soil water limitation and increasing VPD have been linked to increased tree mortality and growth declines (Adams et al., 2017; Grossiord et al., 2020). Recent studies further indicate that the current degree of atmospheric drying is unprecedented compared to the last 400 years (Treydte et al., 2023), making these observations a regional, emerging problem that is largely driven by ongoing climate change. This intensification of water stress not only heightens the vulnerability of these forest ecosystems but also poses a severe threat to tree species with limited adaptability to arid conditions (IPCC, 2023).

Gridded climate data represent regional forcing and do not resolve plot-scale variation in hydrology and microclimate. This scale mismatch may explain why local moisture classes captured variation in BAI and i_WUE not detected by the gridded predictors. SPEI should therefore be interpreted as a regional water-balance index rather than a measure of plot-level soil moisture (De Frenne et al., 2021; Lembrechts et al., 2019).

4.2 Δ¹³C reveals tree responses to environmental stress

Analysis of Δ¹³C in tree rings has shown that it offers valuable insights into the physiological adjustments of trees in response to environmental stress (Brooks et al., 1998; McDowell et al., 2003). Δ¹³C reflects the ratio of intercellular to atmospheric CO₂ concentrations, which in turn is influenced by stomatal conductance and photosynthetic rate (Farquhar et al., 1982). Under water-limited conditions, trees close their stomata to conserve water, but stomatal regulation can vary among species, resulting in interspecies discrepancies in reduced stomatal conductance and limited CO₂ uptake. This water-carbon balance, in turn, leads to changes in the Δ¹³C values fixed within the tree ring (Barbour and Farquhar, 2000; Francey and Farquhar, 1982).

Moreover, the observed changes in Δ¹³C values are indicative of adjustments in i_WUE and carbon assimilation under variable water availability (Ehleringer et al., 1993). Increases in i_WUE suggests that trees optimize water use by regulating stomatal aperture to balance water loss with carbon gain (Walker et al., 2020). The strong linear relationship found between Δ¹³C and moisture-related climate variables, even when using detrended time series, suggests that trees respond physiologically to the interannual variations in water availability and not only to the long-term trends (Figs. 5 and S17–S18).

Although detrended BAI and GI data may not show strong correlations with climate variables due to different limiting factors affecting growth (Cabon et al., 2020; Cuny et al., 2015), the physiological responses captured by Δ¹³C provide evidence of water stress impacting the long-term growth trends (Figs. 4, 5, and 6). This observation highlights an apparent contradiction, where growth data alone may under detect the impacts of water-related physiological stress, emphasizing the importance of integrating physiological indicators (such as Δ¹³C) with growth measurements to achieve a more comprehensive understanding of tree responses.

Moreover, our findings align with those of McDowell et al. (2010), who demonstrated that drought-induced mortality was related to a significant decline in growth alongside low i_WUE sensitivity to water limitation. Here, “i_WUE sensitivity” refers to the magnitude of the temporal trend or slope in i_WUE, where a significant sensitivity indicates a steep trend, and a low sensitivity indicates a nearly flat slope (Fig. 3D). McDowell et al. (2010) observed that surviving trees experiencing water stress exhibited increased stomatal closure, leading to reduced carbon assimilation and growth, and predisposing them to higher mortality risk. However, trees that died showed a larger growth decline but lower sensitivity in i_WUE. Consequently, a sharp decline in growth coupled with modest increase in i_WUE, as observed in Norway spruce relative to Scots pine, suggests that Norway spruce has lower stomatal regulation (i.e., a low water regulation buffer) during periods of increasing water demand. The corresponding abrupt decline in growth is consistent with the dual pattern of growth and Δ¹³C observed elsewhere (McDowell et al., 2010; Wang et al., 2021), suggesting that Norway spruce may be experiencing similar physiological stress that could increase its vulnerability to mortality under continued water-limiting conditions.

Scots pine showed no significant change points in Δ¹³C at any plot, whereas Norway spruce shifted to declining Δ¹³C after 2015–2016 at three of six plots. This difference reflects contrasting temporal response patterns rather than a stronger physiological sensitivity in spruce. Scots pine exhibited a relatively consistent increase in i_WUE throughout the study period, together with stronger Δ¹³C responses to interannual variation in VPD (−4.53 ‰ and −4.89 ‰ per unit VPD before and after 2016, against −1.92 ‰ and −0.57 ‰ in Norway spruce; Table 4), consistent with sustained stomatal regulation under increasing atmospheric water demand: the Scots pine response is gradual and persistent, so there is no structural break to detect. In contrast, Norway spruce maintained relatively stable Δ¹³C values for much of the record before shifting after 2015–2016, producing detectable breakpoints despite comparatively weaker i_WUE responses. The post-2015 decline in Δ¹³C coincided with the pronounced reduction in spruce growth, indicating that temporal changes in Δ¹³C trajectories and climatic sensitivity represent distinct aspects of species responses to increasing drought stress. Moreover, the observed stability in Scots pine growth suggests an effective water-use strategy, enabling it to sustain growth under a long trend of declining water availability, thus avoiding abrupt changes in secondary growth. The findings indicate that Scots pine shows resilience to water limitations, maintaining growth and carbon assimilation (as seen in BAI and GI) while regulating water loss through increased i_WUE.

The variations in both Δ¹³C and i_WUE of the two species underscore their divergent strategies for coping with increasing water limitations. The use of Δ¹³C as a physiological indicator helps our understanding of how trees respond to climatic variability from a mechanistic angle. The value of Δ¹³C in interpreting tree responses to environmental stress has been demonstrated in previous studies (Brooks et al., 1998; Saurer et al., 2014), reinforcing its applicability in the present research.

4.3 Species and soil influences on tree water stress and forest management implications

Norway spruce showed stronger growth sensitivity and weaker i_WUE response than Scots pine under increasing water limitation in southern Finland, highlighting important considerations for the management of these tree species in this boreal region. However, species responses also varied with site conditions, indicating that soil-water interactions modulate drought impacts at finer spatial scales.

At the regional scale, both species showed weak and non-significant tendencies toward stronger breakpoint-defined post-change growth declines on mineral soils compared with organic soils, likely reflecting lower soil water availability and reduced water retention capacity during dry periods (Launiainen et al., 2022; Zignol et al., 2025). Despite this, mean GI patterns differed between species, with Scots pine tending to maintain slightly higher values on mineral soils and Norway spruce on organic soils. Together, these results suggest that broad mineral-versus-organic soil classifications exert only a limited influence on recent growth dynamics relative to species identity and local hydrological conditions.

For Scots pine, the tendency toward slightly lower GI on organic soils may reflect hydrological constraints despite their generally high water-holding capacity. During dry periods, organic soils can experience reduced hydraulic connectivity due to surface drying and the formation of hydrophobic layers, which limit infiltration and root water uptake (Doerr et al., 2000; Hökkä et al., 2025). These processes may constrain access to near-surface water and potentially contribute to reduced growth under drought conditions. In contrast, Norway spruce exhibited slightly higher mean GI on organic soils, whereas mineral soils may expose the species to more frequent short-term water deficits due to lower water storage capacity. This vulnerability may be further exacerbated by the relatively shallow rooting system of spruce under stressful conditions, limiting access to deeper soil moisture (Kalliokoski et al., 2008; Puhe, 2003). Overall, these findings suggest that soil type alone does not exert a dominant control on post-change growth dynamics but instead may interact with species-specific traits and local hydrological conditions to modulate drought sensitivity. This highlights the importance of integrating both species-specific characteristics and fine-scale site conditions when assessing drought risk and developing adaptive forest management strategies.

Focusing on contrasting soil moisture conditions, namely upland versus lowland positions on mineral soils for Scots pine and poorly versus well-drained organic soils for Norway spruce, both species showed weak and non-significant tendencies toward stronger growth declines under drier site conditions, likely due to lower soil water availability. However, clearer differences emerged in mean responses (EMMs), particularly for Scots pine. Scots pine growing on dry upland mineral soils exhibited lower BAI but notably higher i_WUE during the breakpoint-defined post-change period, indicating a physiological shift toward tighter stomatal control under greater water stress. Conversely, Scots pine on wetter lowland mineral soils displayed higher BAI and lower i_WUE, consistent with more favourable growing conditions and reduced need for physiological water conservation. These divergent patterns underscore the critical role of soil moisture regimes in mediating species responses to water limitation (Bose et al., 2020; Green et al., 2019). From a forest management perspective, upland Scots pine stands may benefit from interventions such as reducing stand density, enhancing structural diversity through varied tree ages, sizes, and rooting depths, or promoting mixed-species composition to alleviate water stress and enhance resilience (Aldea et al., 2024; Hanby et al., 2025; Pardos et al., 2021). In contrast, lowland stands could be managed to optimize growth while closely monitoring soil moisture conditions for emerging drought risks. Integrating these spatially explicit strategies will be essential for sustaining Scots pine productivity and stability under intensifying climate change in southern Finland.

Norway spruce growing on poorly-drained organic soils showed lower BAI during the breakpoint-defined post-change period compared with well-drained conditions, indicating reduced growth under saturated soil conditions. Prolonged waterlogging in these soils, where roots tend to remain shallow to avoid anoxic layers, likely limited oxygen availability and impaired root functioning, nutrient uptake, and access to deeper water reserves, thereby increasing vulnerability to subsequent dry periods (Lampela et al., 2023). Consistently, trees on both soil-moisture types showed only modest increases in i_WUE, highlighting a limited physiological capacity to cope with intensified water stress. Consequently, current forest management approaches in drained organic soils of southern Finland, where spruce-dominated stands are widely favoured, may require reconsideration. Promoting mixed coniferous stands could enhance ecosystem resilience, while management strategies that reduce drainage intensity and promote greater water retention in peat soils may help buffer trees against drought impacts (Hökkä et al., 2025; Pardos et al., 2021). Further, practical adjustments to silvicultural regimes in drained peatlands, such as modifying thinning regimes and/or adopting alternative final-harvesting methods including continuous-cover forestry, could reduce intra-stand competition for water by maintaining heterogeneous stand structures and lower peak resource demand compared to even-aged systems (Lehtonen et al., 2023b; Leppä et al., 2020a, b; Tikkasalo et al., 2023). Modifying drainage systems through measures like block ditching and/or partial rewetting may also improve soil water storage (Laine et al., 2011), while continuous monitoring of soil conditions is critical to inform adaptive silvicultural decisions under variable hydrological regimes (Mäkipää et al., 2024). Collectively, these strategies offer promising pathways to improve tree performance and long-term ecosystem stability in drained organic soils as climate-induced water stress intensifies.

5 Conclusions

Our findings demonstrate that forests in southern Finland are increasingly affected by water limitation, particularly during the most recent decade. Norway spruce showed a more pronounced reduction in growth following the post-2016 period, while Scots pine exhibited smaller and more gradual growth declines, evident at the site level from around 2009 but not supported in the regional growth index. In parallel, i_WUE increased more steeply in Scots pine than in Norway spruce, suggesting a weaker physiological adjustment in spruce and reinforcing its greater vulnerability to rising water limitation and prolonged drought stress.

Interannual variability in both growth and i_WUE was generally consistent with changes in SPEI and VPD, highlighting the pivotal role of water limitation in shaping recent growth and physiological responses. Linear mixed-effects models further confirmed that species-specific sensitivity to these hydroclimatic drivers differed, with Norway spruce showing stronger growth reductions and a significant VPD constraint on growth after 2016, while Scots pine appeared more buffered in growth and exhibited stronger increases in i_WUE. Expressed per standard deviation of each predictor, VPD dominated the responses of Scots pine, whereas Norway spruce responded to atmospheric demand and to the summer water balance to a similar degree.

Tree responses to increasing water limitation were primarily structured by species, with weaker modulation by soil type and site-level moisture conditions. Across variables, soil-related differences were generally more apparent in mean values than in temporal trends. Breakpoint-defined post-change period declines in GI tended to be stronger on mineral soils for both species, although not significantly, while mean GI showed weak contrasting patterns, with Scots pine somewhat higher on mineral soils and Norway spruce on organic soils. At the site level, moisture effects were more pronounced: both species tended toward stronger growth declines under dry conditions, but Scots pine maintained higher mean BAI on wet plots, whereas Norway spruce showed higher mean BAI on dry plots, indicating reduced growth under poorly drained conditions. Patterns in i_WUE further emphasized species contrasts, with Scots pine showing higher i_WUE under dry conditions and Norway spruce displaying weak and non-significant responses across moisture regimes.

Taken together, our findings highlight that boreal forest responses to increasing water limitation are primarily structured by species-specific physiological strategies, while being modulated by local hydrological conditions. Accounting for both species differences and fine-scale variation in soil moisture conditions in adaptive forest management strategies will be essential for maintaining forest productivity and resilience under increasingly warm water-limited climatic conditions in northern Europe.

Appendix A: Uncertainty of regional growth indices and its implications on results

Growth indices (GI) were estimated based on tree-ring increment cores collected on Finnish NFI sample plots using linear mixed-effects models (Eq. A1).

$$\begin{array}{}\text{(A1)}& \log \left(g_{\mathrm{pit}}\right)=\mathit{\alpha }+{\mathit{\gamma }}_{\mathrm{t}}+\mathrm{\cdots }+a_{\mathrm{p}}+b_i+{\mathit{ϵ}}_{\mathrm{pit}}\end{array}$$

where g_pit is the basal area increment measured in the year t from the core of tree i on NFI sample plot p, α is the intercept term of the model, a_p and b_i are random plot and tree effects, respectively, and the residuals ϵ_pit are assume to follow an AR(1) temporal autocorrelation structure. The dots (⋯) represent model terms that are expected to explain within-year growth variability, including the impact of age, basal area and relative size of tree i, as well as site fertility class and average temperature sum at the location of plot p. Thus, terms γ_t are expected to capture the annual variation in log-transformed increments attributable to environmental factors such as climate variability.

The annual GI for year t was defined as:

$$\begin{array}{}\text{(A2)}& I_{\mathrm{t}}=\mathrm{100}\exp \left(\mathit{\gamma }{{}^{\prime }}_{\mathrm{t}}\right)\end{array}$$

where $\mathit{\gamma }{{}^{\prime }}_{\mathrm{t}}={\mathit{\gamma }}_{\mathrm{t}}-\stackrel{\mathrm{‾}}{\mathit{\gamma }}$ and $\stackrel{\mathrm{‾}}{\mathit{\gamma }}$ is the average of all γ_t's.

In linear models with categorical predictors, such as year t (treated as a categorical rather than continuous variable) in this case, all parameters are not significant. This is conventionally handled by fixing the effect associated to one of the classes to zero. In this case, γ_t was fixed to zero for the reference year 2022. Selection of reference year has no impact on the estimated GI values, but it must be taken into account when estimating its sampling variance, as explained below.

For simplicity, let us denote the years included in the model by $t=\mathrm{1},\mathrm{2},\mathrm{\dots },T$. As a result of fitting the linear mixed model, we obtained the estimates $\widehat{\mathit{\gamma }}={\left[{\widehat{\mathit{\gamma }}}_{\mathrm{1}}{\widehat{\mathit{\gamma }}}_{\mathrm{2}}\mathrm{\cdots }{\widehat{\mathit{\gamma }}}_{T-\mathrm{1}}\right]}^T$ (${\widehat{\mathit{\gamma }}}_T=\mathrm{0}$ not included) and their covariance matrix Σ.

Since

$$\begin{array}{}\text{(A3)}& {\widehat{\mathit{\gamma }}}^{\mathbf{\prime }}={\left[\widehat{\mathit{\gamma }}{{}^{\prime }}_{\mathrm{1}}\widehat{\mathit{\gamma }}{{}^{\prime }}_{\mathrm{2}}\mathrm{\cdots }\widehat{\mathit{\gamma }}{{}^{\prime }}_T\right]}^T=\mathit{A}\widehat{\mathit{\gamma }},\end{array}$$

where

$$\begin{array}{}\text{(A4)}& {\mathbf{A}}_{T\times \left(T-\mathrm{1}\right)}=\left[\begin{array}{cccc}\mathrm{1}-\frac{\mathrm{1}}{T}& -\frac{\mathrm{1}}{T}& \mathrm{\cdots }& -\frac{\mathrm{1}}{T}\\ -\frac{\mathrm{1}}{T}& & & \mathrm{1}-\frac{\mathrm{1}}{T}\\ \mathrm{⋮}& & & \mathrm{\ddots }\\ -\frac{\mathrm{1}}{T}& & & \mathrm{1}-\frac{\mathrm{1}}{T}\\ -\frac{\mathrm{1}}{T}& -\frac{\mathrm{1}}{T}& \mathrm{\cdots }& -\frac{\mathrm{1}}{T}\end{array}\right]\end{array}$$

the variance-covariance matrix of $\widehat{\mathit{\gamma }}{}^{\mathbf{\prime }}$ is

$$\begin{array}{}\text{(A5)}& {\mathbf{\Sigma }}^{\mathbf{\prime }}=\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^T.\end{array}$$

The non-linear transformation (Eq. A2) was handled with a Monte Carlo approach:

• 1000 realizations of parameter vector γ^′ were simulated from the multivariate normal distribution with mean vector $\widehat{\mathit{\gamma }}{}^{\mathbf{\prime }}$ and covariance matrix Σ^′,

• the transformation was applied to each simulated γ_t^′: $I_{\mathrm{t}}{}^{\prime }=\mathrm{100}\exp \left({\mathit{\gamma }}_{\mathrm{t}}{}^{\prime }\right)$, and

• 95 % confidence interval was constructed from the 0.025- and 0.975-quantiles of the 1000 simulated values of I_t^′.

The GIs and their 95 % confidence intervals were estimated by tree species (Scots pine and Norway spruce) and soil type (mineral and organic; Fig. A1). Although the standard errors of GIs are relatively large on organic soils, especially for Norway spruce, and increase towards the end of the series due to the decreasing sample size (Fig. 2), interannual variability in them is estimated with relatively high precision (Fig. A2 left). This is explained by the strong positive error correlations of the estimate for subsequent years (Table A1) resulting from the tree-level random effects (the core from one tree contributes to several years' GIs).

The effects of GI uncertainty on breakpoint detection (Fig. 4A–B) were assessed by fitting the segmented regression models to each of the 1000 simulated GI series in the same way as to the estimated GIs (Sect. 2.4.2). Detection rates differed strongly among spp. and soil types. For Scots pine, a significant breakpoint was detected in only 9 out of 1000 simulations (0.9 %) on mineral soils and in 320 of 1000 (32 %) on organic soils, where breakpoints were detected on organic soils they clustered near 2019 (median 2019.4, 95 % range 2019.0–2020.1) rather than near 2009. For Norway spruce, every simulation on mineral soils yielded exactly one significant breakpoint (100 %; median 2016.3, 95 % range 2016.0–2016.7). On organic soils the outcome was more variable: no breakpoint in 9.7 % of simulations, one breakpoint in 43.6 %, and two breakpoints in 46.7 %, with the breakpoint nearest 2016 having a median of 2016.8 and a 95 % range of 2009.3–2017.7 (Fig. A2 right). The probability that a breakpoint exists is therefore 1.00 for Norway spruce on mineral soils and 0.90 on organic soils, against 0.01 and 0.32 for Scots pine on mineral and organic soils, respectively. The species contrast reported in Sect. 3.1 and 3.3 is thus a property of the sampling distribution of the growth index and not of any single realisation.

Two features of the procedure should be stated. First, the breakpoint search implemented for the submitted analysis stopped after two breakpoints, so simulations recorded as having two breakpoints were right-censored rather than having stopped because no further break was significant. Re-running the same simulations without this limit shows that 33.0 % of the Norway spruce organic-soil simulations in which a breakpoint is detected support three or more, and that 98.7 % of those that stopped on the limit prove to have three or more once it is removed; these simulations describe a series with several changes of slope rather than a single shift around 2016. Detection on the observed NFI series is unaffected: without the limit the procedure returns exactly one breakpoint for both Norway spruce series and none for either Scots pine series, identical to Table S19. Second, the original simulation was run without a fixed random seed, so the 1000 realisations reported here are those preserved from the original run rather than a regenerated set; a seed has been added to the archived code.

Confidence intervals for segmented regression predictions (Fig. A3) were constructed from the 0.025- and 0.975-quantiles of the pointwise predictions over the 1000 models based on the 1000 simulated GI series. For spruce on organic soils, the slope was positive in 98 % of simulations in 2015, 66 % in 2016, and 18 % in 2017.

The same 1000 realisations were carried through the segmented regressions and through Eq. (6), so that the results reported in the main text are not conditional on a single detected breakpoint. Pooling by Rubin's rules over the 1000 realisations, the Norway spruce GI trend on mineral soils is $+\mathrm{3.11}\pm \mathrm{0.52}$ yr⁻¹ before the breakpoint (p<0.001) and $-\mathrm{8.69}\pm \mathrm{2.79}$ yr⁻¹ after it (p=0.006), with a breakpoint detected in every realisation (median 2016.3). On organic soils the pooled point estimates are of the same sign and comparable magnitude ($+\mathrm{4.62}\pm \mathrm{6.07}$ and $-\mathrm{8.55}\pm \mathrm{7.04}$ yr⁻¹) but neither differs significantly from zero (p=0.48 and p=0.26), a breakpoint is detected in 90 % of realisations, and its location is less well constrained (95 % range 2009.3–2017.7); the regional decline in Norway spruce growth after 2016 is therefore resolved by the NFI sample on mineral soils, whereas on organic soils the direction of change is consistent but its magnitude is not (Table A2). Carried through Eq. (6), the pooled fixed effects retain the sign and significance class of the estimates reported in Tables 3 and 4 for both species and both response variables: for Norway spruce the SPEI effect on BAI is −0.036 (95 % CI −0.065 to −0.008, p=0.013; significant in 100 % of realisations) and the VPD effect on Δ13C is −1.84 (95 % CI −3.08 to −0.60, p=0.004; significant in 98 %), and for Scots pine the VPD effects on both BAI and Δ13C remain significant in 100 % of realisations. The shift in intercept between periods is the one quantity that does not survive: for Norway spruce BAI it changes sign relative to the fixed-2016 fit (+0.033 fixed against −0.010 pooled, p=0.911), is significant in only 30 % of realisations, and has 66 % of its variance attributable to the breakpoint location, while for Δ13C it weakens from p=0.002 to p=0.045 conditional on detection and to p=0.079 under the full model average (Table S21). A fixed-split sensitivity analysis at 2015, 2016 and 2017 gives the same picture (Table S22).

Figure A1 NFI-based standardized growth index (GI) for Scots pine and Norway spruce from 1999 to 2022 on mineral and organic soils (cf. Fig. 2). Error bars represent the 95 % confidence intervals depicting the uncertainty due to NFI sampling.

Figure A2 Left: Ten realisations of GI series for Norway spruce on organic soils, simulated from the distribution determined by the NFI sampling variance and correlations of the GI estimates. Right: Relative frequencies of all break points of segmented regression models fitted to the 1000 simulations from the uncertainty distribution of GIs for spruce on organic soils.

Figure A3 NFI-based standardized growth index (GI) and confidence intervals (95 %) for segmented regression predictions, accounting for NFI sampling uncertainty. Points show GI; lines denote confidence intervals.

Table A1 Correlation matrix of sampling errors of the GI estimates for Norway spruce on organic soils for the period 2015–2022.

Table A2 Segmented-regression estimates for the NFI growth index with NFI sampling uncertainty propagated. Each of the 1000 simulated GI series was fitted in the same way as the estimated GI (Sect. 2.4.2); slopes are pooled by Rubin's rules over the realisations in which a breakpoint (BP) was detected and located inside the 1999–2022 window. P(BP) is the proportion of the 1000 realisations yielding at least one significant breakpoint. The last column gives the number of realisations contributing to the pre- and post-breakpoint slopes; it is smaller than 1000×P(BP) where a detected breakpoint fell outside the analysis window. Post-breakpoint slopes for Scots pine rest on very short segments and few realisations and should not be interpreted as evidence of a regional trend.