Increasing atmospheric dryness reduces boreal forest tree growth

Domain: earth_science
Published: 2023-10-30
Editorial outcome: minor_revision

Abstract
Rising atmospheric vapour pressure deficit (VPD) associated with climate change affects boreal forest growth via stomatal closure and soil dryness. However, the relationship between VPD and forest growth depends on the climatic context. Here we assess Canadian boreal forest responses to VPD changes from 1951-2018 using a well-replicated tree-growth increment network with approximately 5,000 species-site combinations. Of the 3,559 successful growth models, we observed a relationship between growth and concurrent summer VPD in one-third of the species-site combinations, and between growth and prior summer VPD in almost half of those combinations. The relationship between previous year VPD and current year growth was almost exclusively negative, while current year VPD also tended to reduce growth. Tree species, age, annual temperature, and soil moisture primarily determined tree VPD responses. Younger trees and species like white spruce and Douglas fir exhibited higher VPD sensitivity, as did areas with high annual temperature and low soil moisture. Since 1951, summer VPD increases in Canada have paralleled tree growth decreases, particularly in spruce species. Accelerating atmospheric dryness in the decades ahead will impair carbon storage and societal-economic services.

Introduction
Climate change will alter boreal forests capacity to store carbon, all the more so because they are dominated by cold-tolerant species and are warming faster than most other land areas (0.2 °C to 0.5 °C per decade over 1961–2015)1,2. Warming is expected to alter the net annual carbon uptake of boreal forests through changes to temperature-related variables (growing season length, drought severity, vapour pressure deficit, freeze frequency) that affect tree mortality, recruitment, physiology and growth3–5. Across Canada’s boreal forest, land carbon storage capacity varies due to local and regional differences in climate, vegetation, soils, surficial geology, and disturbance regime6. The complexity of the boreal environment and the diversity of boreal species produce notable contrasts in productivity responses to warming6–8. Across boreal Canada, forest productivity is anticipated to decrease in the west because of water stress, and increase northeastward where temperature is currently the main factor limiting tree growth9–11. Such mechanisms have been raised as potential drivers of forest greening or browning in Canada12,13. Warming can thus enhance plant growth in cool and wet environments, but impair it above a threshold of water stress7,14. Additionally, differential trajectories among tree genera, species and ecoclines are possible, depending on their stress tolerance and plasticity15–18. Changes in growth rates of boreal forests are of great concern due to the critical roles these forests play in global carbon storage, human population support and economic-industrial services19.

Atmospheric vapour pressure deficit (VPD), the difference between the amount of water vapour in the atmosphere and the potential amount held at saturation20, is a major temperature-related determinant of plant physiology. Without a corresponding increase in the actual amount of atmospheric water vapour, VPD will increase with warming because the saturation vapour pressure is a curvilinear function of air temperature21. High VPD can induce high leaf and soil water loss, increasing the risk of plant water deficits and stress22,23. Water deficits increase xylem tension and tissue cavitation: these negative effects are minimized by stomatal closure, which decreases water loss at the cost of carbon uptake, leading to growth reductions and eventually tree mortality20,21,24,25. From a hydrometeorological perspective, higher atmospheric dryness also increases atmospheric demand for water from the land surface, increasing evaporation and reducing available soil moisture26. Hence, plant growth responses to atmospheric dryness can result from direct effects on stomatal conductance and indirect effects on soil moisture through increasing evapotranspiration27.

Tree species have evolved different strategies to cope with atmospheric dryness and therefore have differential responses to VPD changes, modulated by environmental conditions20. There is thus a need for in-situ, multi-species comparisons to quantify these differences across a broad range of environmental conditions25,28. Recent remote sensing studies have produced mixed results about the relative role of VPD on productivity and carbon uptake26,29. Tree-ring studies offer an important opportunity to understand VPD impacts, especially pertaining to aboveground carbon uptake. To date, tree-ring analyses directly linking changes in dryness with boreal forest growth have only focused on a few species and areas30. Integrating tree-ring sampling into national forest inventories has enabled large-scale ecological assessments related to forest health and carbon cycling, including in Canada31–33. Quantifying large-scale tree growth responses to VPD would improve our understanding of forest growth responses under climate change and reduce uncertainties in the model predictions used to evaluate strategies for mitigating and adapting to climate change.

Here, we assess boreal forest responses to changes in atmospheric VPD using a well-replicated tree-ring network covering Canada’s forests over the period 1951–2018. The growth-VPD relationships enabled mapping of spatially-explicit VPD responses across Canada’s boreal zone. We then explored the main drivers of differential growth responses to VPD, including species, local precipitation and temperature, elevation, and tree age and size. Finally, we determined how VPD and growth are changing over time for the most responsive species. Our results indicate that increased VPD has already had a broad-scale negative effect on tree growth in Canada’s forests over recent decades, specifically for younger trees and widespread spruce species.

Results
Tree response to atmospheric vapour pressure deficit
Utilizing tree-ring data obtained from 32,189 trees spanning Canada (Fig. 1), we developed 4931 Generalized Additive Mixed Models (GAMMs) based on species and site factors. These models were employed to evaluate the impact of VPD on the annual basal area increment (BAI, the increment in cross-sectional area of trees). We achieved convergence in 3559 (72%) of the species-site GAMM models, indicating successful outcomes. The average goodness-of-fit between year-to-year observed growth fluctuations and growth predicted from site-specific growth-VPD models was r²=0.51 (σ±0.21, n=3559 convergent models). Among the convergent species-site GAMMs, 58% (2057 models) showed a significant relationship between BAI and VPD_t, VPD_t-1, or both across the nine species analyzed (Picea mariana, Picea glauca, Pinus banksiana, Populus tremuloides, Pinus contorta, Pseudotsuga menziesii, Picea engelmannii, Abies lasiocarpa, and Pinus resinosa).

Among the convergent species-site GAMM models, 31% (1,096 models) had a t-value for the relationship between annual BAI and VPD of the year of ring formation (VPD_t) achieving the p<0.05 threshold. Of these 1,096 significant relationships, 752 (about two-thirds) were negative, indicating that increasing VPD is detrimental to tree growth for these sites-species combination (Fig. 2, Table S1). The majority of negative t-values were found near warm, dry margins of the boreal forest in south central Canada (i.e, southern part of the Boreal Shield ecoregion), specifically for Picea glauca and Picea mariana (Figs. 1–3, and Supplementary Materials Fig. S1). A third of the significant BAI-VPD t-value were positive (344 models), indicating that increasing VPD was positively related with current year tree growth (Fig. 2, Table S1). The positive t-values were found in western Canada, i.e Boreal Cordillera, and in cooler moister areas of central and eastern Canada, i.e Hudson Plain and Taiga Shield (Figs. 1–3, Fig. S3). Both negative and positive t-values were found for all species except Pseudotsuga menziesii, the distribution area of this species does not overlap the previously cited ecoregions. Furthermore, a significant relationship with VPD of the year prior to ring formation (VPD_t-1) was observed in 47% of the convergent species-site GAMM models (1,660 models). In these models, t-values were almost exclusively negative (96%, Fig. 2, Table S1). Collectively, our findings suggest that VPD, whether from the previous year or the current year, had a significant and detrimental impact on BAI in approximately 51% of the 3,559 species-site combinations where models converged (roughly 37% of the initial 4,931 candidate species-site models). Based on these observations, VPD responsiveness of tree growth was found in one-third to half of the species-site combinations. The consistency of the results was maintained in the partial BAI-VPD GAMM models, which aimed at controlling for the impact of summer soil moisture index (SMI) on growth, as demonstrated in the Supplementary Materials.

Identify the determinants of VPD growth response
We assessed the determinants of growth response to VPD through the relationship between the significant t-values of growth-VPD relationship and seven environmental, biological and forest structure variables, namely site elevation, mean annual temperature (MAT), mean annual precipitation (MAP), summer SMI, tree species, mean tree age and BA at the site level.

Among the significant variables, average depth in the random forest decision trees was lowest for species identity. Species identity was selected as root node in 158 of 500 trees, meaning it is the first variable determining the direction and strength of growth-VPD t-values (Fig. 4, Table 1). The species showing more negative VPD growth responses were, in order of most to least sensitive, Pinus resinosa, Pseudotsuga menziesii, Pinus contorta, Picea glauca, Picea mariana, and Pinus banksiana, while Populus tremuloides, Picea engelmanii and Abies lasiocarpa showed a less negative response to VPD. Site-level mean tree age displayed the second lowest average depth and was selected as a root node in 68 of 500 decision trees (Fig. 4, Table 1). Mean tree age was generally positively correlated with t-values (Fig. 5), meaning a more negative response to VPD occurred in young to mature trees (0–100 years) compared to mature and old trees (100–250 years); a low sample size hinders the capacity to infer the typical growth response to VPD beyond 250 years of mean tree age. Site-level mean BA displayed high average depth and very low MSE, meaning it makes a minor contribution to t-value variations.

Local mean annual temperature (MAT) was the primary environmental and climatic determinant of t-values, with the third lowest average depth among random forest trees and second highest selection as root node (for 117 of 500 trees). MAT was negatively correlated with t-values, meaning VPD had a more negative effect on growth in higher MAT environments or during warmer years (Fig. 5). Long-term mean summer SMI showed the fourth lowest average depth and was chosen as the root node in 96 decision trees (Fig. 4, Table 1). Mean summer SMI was positively correlated with t-values, meaning the most detrimental effects of VPD were seen in sites with low available soil water (Fig. 5).

Spatiotemporal changes in VPD and growth of widespread tree species
Central and northern Canadian boreal regions have experienced significant increases in summer VPD since 1951 (Fig. 6a). The affected regions extend across the ranges of Picea mariana and Picea glauca (Fig. 1, Fig. S1). At these locations, mean VPD values increased from slightly above 0.5 kPa in the 1950s–1960s to >0.6 kPa in the 2000s. We found that percent growth changes (GC, see Methods) in terms of BAI, averaged by species, were linearly and inversely related to prior-summer VPD (Fig. 6b, c). From the slope of this relationship, we estimated that for every ~0.1 kPa increase in VPD, there was approximately a 10 to 11% decrease in mean BAI throughout the geographic ranges of these spruce species (Fig. 6b, c). Decreases in tree growth from 1951 to the present were apparent in both species, paralleling the increases of summer VPD (also see Fig. S5). Although the northern part of the region was sparsely sampled among our ~3500 sites, because the ring-width of both species is negatively affected by high prior-summer VPD values (Figs. 3, 6), we might also expect a negative growth response for these species in this region in the decades ahead. We also examined the same relationship for two other species that showed significant growth sensitivity to VPD, Pinus contorta and Picea engelmannii (Fig. 6d, e). Again, the percent change in growth averaged across Pinus contorta was linearly correlated with prior-summer VPD, but to a much lesser extent than in the two spruce species. For Picea engelmannii, there was no evidence for a negative relation of range-wide averaged growth change to prior-summer VPD (r=‒0.18 with 95% CI [−0.44; 0.12]), but rather a moderate positive relationship to current-summer VPD (Fig. 6d).

Methods
Study area
The study area covers eleven forested ecozones (Ecological Stratification Working Group (ESWG), 1996) in Canada (Fig. 1). These ecozones are the Montane Cordillera in the west coast, Boreal Cordillera and Taiga Cordillera that contain the Canadian Rockies, the Boreal Plains and Taiga Plains located east of the Rockies, and Prairies, Taiga Shield, Boreal Shield and Hudson Plains in the continental interior to the east. Bounded by three Great Lakes in the east are the Mixedwood Plains, and adjacent to the east coast is the Atlantic Maritime (Fig. 1). The climate in the boreal/hemi-boreal zone is predominantly high-latitude continental, with long cold winters, short cool summers, and relatively low annual precipitation, but with significant regional variation9. Summer mean air temperature averages from 10 °C in the Taiga Cordillera to 14 °C in the Atlantic Maritime, whereas winter mean air temperature ranges from −22.0 °C in the Taiga Cordillera to −1.5 °C in western Montane Cordillera. Annual precipitation totals average from 200-500 mm in the Taiga Plains to over 4,000 mm in western Montane Cordillera.

Climatic data
Daily maximum and minimum temperature (°C), precipitation (mm) and relative humidity (%) were obtained for the period from 1950 to 2018 using BioSIM v10.3 software, which interpolates site-specific estimates from all of Environment and Climate Change Canada’s historical daily weather observations58,59 (compiled weather data available at ftp://ftp.cfl.forestry.ca/regniere/Data/Weather/Daily/; last accessed 2023-06-20). Using daily minimum and maximum temperature as well as precipitation data, we employed the BioSIM software to calculate the daily vapour pressure deficit (VPD, in kPa) following the approach described in Allen et al.60. For dewpoint temperature calculation, BioSIM utilizes the methodology outlined in Kimball et al.61 and leverages daily temperature and precipitation data. The soil moisture index (SMI, in % of water holding capacity) was estimated in BioSIM using the quadratic-plus-linear (QL) formulation procedure described in Régnières et al.58, which accounts for water loss through evapotranspiration (simplified Penman–Monteith potential evapotranspiration) and gain from precipitation. Following Girardin et al.10,31, we established the parameters for critical and maximum available soil water as 300 mm and 400 mm, respectively. A low SMI is an indication of low available soil water at the site. We acknowledge that the SMI metric used is a simplified representation of reality and entails inherent uncertainty due to the lack of detailed soil attribute representation for parameterizing site-specific water holding capacity. Daily temperature, precipitation, VPD, and SMI data were all interpolated to a 1° x 1° grid (n=1705 grid points). In contrast to many other regions worldwide, Canada has a limited density of weather stations (Fig. S6), which hampers the ability to derive precise climatological information at specific sites. Although broader-resolution gridded data may not adequately capture elevational climate gradients, they do offer the advantage of homogeneity and long-term time series, which are essential for establishing climate–growth relationships62. VPD and SMI data were averaged for the summer season (June to August) by year, and attributed to the closest tree-ring sampling sites’ grid point. In addition, long-term means of annual climatologies (1951–2018) were also computed for each site (mean annual January to December temperature, MAT; mean annual January to December sums of precipitation, MAP; and mean summer soil moisture index, SMI).

Tree-ring data
Annual tree-ring width data were retrieved from the Canadian Forest Service Tree-Ring repository (CFS-TRenD 1.0)31, developed with the goal of combining data from different sources and making them available in a consistent format for large-scale analyses. CFS-TRenD contains measurements from 40,206 samples from 4,594 sites and 62 tree species. The primary national-scaled dataset in CFS-TRenD 1.0 is increment cores sampled since 2001 during the establishment of Canada’s National Forest Inventory (NFI)63. The NFI network, designed to represent species distributions and their range of growing conditions in Canada, comprises 6,010 core samples from 870 sites. Other important data are more regional and include data from a network of permanent sample plots established by the Alberta Biodiversity Monitoring Institute (ABMI), a network of temporary sample plots established by the Ministère des Ressources Naturelles et de la Faune du Québec (MFFPQ; Programme d’inventaire écoforestier nordique)31,64, and the Climate Impacts on Productivity and Health of Aspen (CIPHA) network65. Additional contributions to CFS-TRenD are from smaller-scale projects carried out by individual researchers, and data extracted from the International Tree-Ring Data Bank (ITRDB).

Analyses were limited to the nine dominant species (seven genera) in the dataset (Fig. 1) to maintain a sampling density that balances the local variations in growth across sampling sites within regions, and detects growth sensitivity to climate at the regional scale. These nine species represent 80% of the samples (32,189 trees) in the CFS-TRenD repository: Picea mariana (black spruce; 25% of the dataset with 1.0e4 trees sampled), Picea glauca (white spruce; 12%, 4.7 e3 trees), Pinus banksiana (jack pine; 11%, 4.4 e3 trees), Populus tremuloides (trembling aspen; 10%, 3.9 e3 trees), Pinus contorta (lodgepole pine; 7%, 2.9 e3 trees), Pseudotsuga menziesii (Douglas fir; 7%, 2.8 e3 trees), Picea engelmannii (Engelmann spruce; 4%, 1.7 e3 trees), Abies lasiocarpa (subalpine fir; 2%, 9.3 e2 trees), and Pinus resinosa (red pine; 2%, 8.7 e2 trees). We excluded Abies balsamea (balsam fir, mostly dominant in the eastern ecozones) because its growth dynamics are strongly influenced by spruce budworm (Choristoneura fumiferana) defoliation18,66. All of the nine species occur in Canadian boreal/hemi-boreal forests6. Half of the sites included two or three sampled trees, and 83% of the sites included a single species. Fifty percent of the trees displayed between 40 and 90 measured tree rings (i.e 25th and 75th percentiles of the age distribution). To ensure the accuracy of the measured tree rings and the assignment of the correct calendar year to each tree ring, both visual and statistical quality control procedures were implemented10,31. Individual contributors to CFS-TRenD typically conducted quality control using the COFECHA program67, or analogous approaches that assess tree sample correlations against species-specific reference chronologies. Subsequently, the harmonization of the different contributions was performed within the R v4.2.2 environment31. In an analysis of sample coherency, Girardin et al.31 found high spatial synchronicity in the interannual growth fluctuations among most samples.

Responsiveness of tree growth to vapour pressure deficit
Calculation of basal area increments. Tree-level average ring-width series were converted to annual basal area increments (BAI; cm2 yr-1), which is recognized to be closely related to tree productivity68. BAI was estimated from the relation:

BAI_t = \pi R_t^2 - \pi R_{t-1}^2 (1)

where R_t and R_{t-1} are the stem radii (cm) at the end and beginning of a given annual ring increment. Minimum tree age was determined from ring counts, starting from the outermost ring. Considering the atypical response to environmental drivers, particularly the influence of overstory tree competition on juvenile stems, tree rings formed during the initial ten years were excluded from the analysis. This decision was made based on the understanding that these early rings often exhibit distinct growth patterns that are less representative of the major environmental drivers69.

Growth and climate relationships. Species-specific growth-VPD relationships were determined using generalized additive mixed models (GAMM) fitted at each site between the log-transformed BAI series, tree basal area of the year prior to ring formation and atmospheric VPD values during summer (June, July and August) of the prior and current year of ring formation, with tree considered as a random effect:

log(BAI_{tjk}) = \alpha_j + s_1(BA_{t-1jk}) + s_2(age_{tjk}) + \beta_{1j}VPD_t + \beta_{2jk}VPD_{t-1} + Tree ID_{jk} + \epsilon_{tjk} (2)

where j stands for the species, k for the site, and t for the year. BA is the basal area, BAI is the basal area increment, age is the age in years, and s is a cubic regression spline smoothing parameter whose degree of smoothness was determined through an iterative fitting process. Temporal autocorrelation was considered with AR1, an autoregressive term of order 1 accounting for year t and Tree ID (each tree’s unique identifier). Approximately 20% of the sampled cores either did not contain tree pith or were estimated to be located more than a centimeter away from it. Since we do not always have access to the original samples or scanned images, we were unable to apply a correction method for the pith offset. Nevertheless, we anticipate that the tree-level random effect (Tree ID) would account for this missing information. If the error magnitude is important or if there is substantial noise resulting from forest management activities or abiotic/biotic disturbances, the convergence of the species-site model is anticipated to be unsuccessful. As a result, the affected site will be excluded from subsequent analyses. The significance of variables at the 5% level was determined from t-tests in GAMM models (t-value’s p<0.05). The growth model was fitted using the ‘mgcv’ v1.8.41 R package70.

Drivers of growth sensitivity to VPD. Environmental, tree species and forest structure variables responsible for the observed distribution in VPD-growth relationships (i.e t-values) were assessed as follows. These analyses were guided by the hypothesis that the strongest negative responses to VPD would be observed at sites exposed to high temperatures or with low soil water availability. We also postulated that there is differentiation among species’ responses to VPD resulting from their evolutionary strategy for stomatal conductance. We considered four environmental variables (elevation, mean annual temperature (MAT), mean annual precipitation (MAP), and summer SMI), tree species, and two forest structure variables (mean site age and basal area (BA)). First, we used the Random Forests (RF) algorithm, a machine learning method adapted to complex and potentially non-linear relationships71. Our RF was based on the bootstrap of 500 training decision trees aggregated to predict the t-values of the significant (p<0.05) growth-VPD_t or growth-VPD_t-1 relationships, using these seven predictor variables as the potential explanatory variables. The running and interpretation of the RF were performed using ‘randomForest’ v4.7.1.1 and ‘randomForestExplainer’ v0.10.1 R packages72,73. Variable importance was measured as the decrease of unscaled mean square error (MSE decrease) observed when the variable is randomized, and averaged among training trees74. We retrieved variables’ average depth among training trees and their occurrence as the first node. Finally, to get the direction of the relationship between variables and t-values, we drew the partial dependence plot representing the marginal effect of each variable on the t-values. To untangle the direct VPD effect of atmospheric dryness on stomatal conductance from its indirect effect on soil moisture, we examined partial BAI-VPD GAMM models obtained after controlling for the effect of soil moisture index (SMI) on growth (see Supplementary Materials and Fig. S4).

Trends in VPD and growth. Annual fluctuations in growth over time were computed for Picea mariana, Picea glauca, Pinus contorta and Picea engelmanii that are among the primary harvested species in Canada, and which geographic distribution is particularly threatened by climatic changes19. Annual growth fluctuations were obtained from the detrending of BAI series that removed growth trends due to tree age and size75. We detrended time series of BAI using species-specific GAMMs fitted to the log-transformed BAI:

log(BAI_{tj}) = \beta_{jk} : log(BA_{t-1j}) + s(age_t) + corAR1_{j} (3)

Annual growth changes (GC), expressed as the percent deviation from predicted values generated by the GAMMs, were then computed following Girardin et al.76. GC was average by species and year to show yearly temporal variability in species-specific GC since 1951, see Fig. S5 for models fitting statistics. The average GC 95% confidence interval was computed at the species level, for each calendar year, by correcting the effective degrees of freedom (n′) based on first-order (i.e, lag=1) autocorrelation estimates of Moran’s I77 (moran.test function in R). Correlations (r) between growth changes and VPD were computed from a bootstrapping technique that accounted for autocorrelation and trends in data78, and regression slope (\beta) from ordinary-least-square regressions.

Linear trends of summer VPD were examined for 1951–2018 using least squares linear regressions. A derivation of a t-test with an estimate of the effective sample size accounting for serial persistence in data returned the slope’s significance against the null hypothesis that the trend is zero79.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_3": "Second Review of Mirabel et al. \"Increasing atmospheric dryness reduces boreal forest tree growth\"\nI appreciate that the authors have made an effort to alleviate my concerns. The required additional information regarding the climate data has now been added to the Methods section and (some of) the figures have been improved.\nBased on the information provided, I think it is very likely that the results are robust w.r.t. the climate data, but some of the methodological choice are still not ideal in my opinion (primarily w.r.t. gridding).\nMajor Comments:\nMy concerns about the source of the climate data have been largely addressed with additional explanation and references in the Methods section, although some mention of the fact that most of the humidity data are estimated from minimum temperature data in an earlier section would benefit the reader, I think.\nMy concern about the resolution and gridding of the climate data was mostly addressed and I appreciate that the authors did provide results based on direct interpolation of climate data to sites, but I do not agree with the author's rational for gridding.\nBased on the data provided, the choice of interpolation does not have a large impact, but I would have liked to see, for example, how many of the GAMM models converged with the site-based interpolation (I suspect maybe more). It would have been good to see t-value maps like those in Fig. 3 with site-based estimates as well.\nThere is some change apparent in Fig. R3 and I would consider these results closer to the truth than what is shown in Fig. 4 in the manuscript, but I would agree that for the most part, the differences are probably not relevant/meaningful. It is interesting, though, that Mean Age notably decreased in importance (as measured by MSE).\nMore broadly, however, I strongly disagree with the authors rational for choosing gridding in the first place, and I believe that using station-based interpolation would be the methodologically more appropriate approach (and that is what I meant to suggest, not interpolation to a very dense, hi- res grid, although I guess I could have been clearer). It is true that weather stations in Canada are sparse, but if one is interested in specific sites, there is no point in interpolating to a regular grid first and then choosing the closest point, when one could just use, e.g. IDW to average available stations in the vicinity. This is even more so the case when lapse-rate adjustment is applied. Every additional step of interpolation introduces more uncertainty.\nInterpolation to a regular and appropriately-spaced grid is only necessary, when area averages are being considered, where representativeness of a station for an area is of concern.\nConcerning the correlation between climate variables on the regular grid and at site locations: I am almost surprised that it is not higher. Correlation only measures variability and does not capture constant biases, but the main difference I would expect between site values and regular grid points would be in the lapse-rate adjustment, which primarily manifest in relatively constant deviations. RMSE or bias would have been a better measure here.\n(I am assuming the authors are referring to temporal correlation.)\nSimilarly\nMinor Comments:\nFig. R2: I am not sure what the authors are trying to show... there are clearly areas in the southern part of the domain where station density is much higher than 60 - 100 km and a finer grid would be justified (Fig. R1). At the same time Fig. R2 also clearly shows that there are regions in the Rocky Mountains where the coarser grid misses regionally differing trends that are present at higher resolution, and which would reduce noise in the analysis, if properly captured.\nIf the authors are concerned about some regions having lower levels of precision (more noise) than others due to station density: that is already the case due to strong topographic gradients in the Rocky Mountains, and gridding to a 1 degree rid exacerbates this issue.\nFigure 1: The figure looks significantly better, but the resolution could still be improved.\nFigures 2 & 3: These figures are much improved and I think the decision to move most of the species-dependent plots into SI was good.\nWhat was the selection criterion for the species shown in the main manuscript?\nFigure 6: I am not sure if the change in color map does anything to improve the readability of the map. In particular, the large water bodies on the map are shown in a color that is also part of the color map, and at first glance it was not apparent to me that this was simply a geographic overlay - I think a color that is not part of the color map should be used here, or simply a thin black outline. Furthermore, the authors state that the resolution had been increased, however, I struggle to see that.\nI think it would also be appropriate to mention in the context of Fig. 6 that the underlying humidity data were largely estimated from Tmin (using Kimball's method).", "reviewer_4": "As requested by the editor I do not give an additional full review of the paper - I feel the three reviewers before did already an excellent job and identified a whole set of unclarities. A fully agree with my co-referees that the paper provides novel and important information showing that not only in rather warm and already (from the view of the atmosphere and the soil) dry regions (as e.g. by Williams et al. 203) increases in VPD are already today (and even more in future) a threat for tree growth and functioning. The amount of tree ring data is impressive and it is central that not only species difference but also effects of age and region were accounted for.\nThe response of the authors to the points raised by the three referees:\nAs mentioned below the referees not really identified errors in interpretation or flaws in methodology bur mainly ask for clarification of several points were either some information was missing in the original manuscript or was described in a rather ambiguous way. In the revised version the author did and excellent job and clarified - in my opinion - all points that were put forward by the reviewers. This holds true for the technical questions of Paolo Cherubini, related to BAI calculations, cross dating and additional inclusion of ITRDB data as well as for the clarification of the VPD and humidity data as well as of the BioSIM estimates as requested by referee 3. From my point of view the results are novel and show the VPD dependency has been shown based on a really large data set and the authors reacted very adequately to the points raised by the reviewer."}, "author_response": "We appreciate your time and advice regarding the information added after the first round of reviews. Below, we will address your concerns and provide all the necessary materials to properly support our methodological choices.\nMajor Comments:\nMy concerns about the source of the climate data have been largely addressed with additional explanation and references in the Methods section, although some mention of the fact that most of the humidity data are estimated from minimum temperature data in an earlier section would benefit the reader, I think.\nWe are pleased to hear that the changes we made to the methods have enhanced readability. In accordance with your advice to clarify our VPD calculation method as early as possible in the manuscript, we have applied the suggested changes to the caption of Figure 6 (please see the later comment for details).\nMy concern about the resolution and gridding of the climate data was mostly addressed and I appreciate that the authors did provide results based on direct interpolation of climate data to sites, but I do not agree with the author's rational for gridding.\nBased on the data provided, the choice of interpolation does not have a large impact, but I would have liked to see, for example, how many of the GAMM models converged with the site-based interpolation (I suspect maybe more). It would have been good to see t-value maps like those in Fig. 3 with site-based estimates as well.\nTo further justify our choices, we tested the convergence rate of models using the site-based interpolation and corresponding maps. With this site-based dataset, 72% (3,572 models) of the species-by-site models converged, which is 13 more sites compared to our previous results where we had 3,559 converging models (same percentage).\nWe also reproduced the maps as was done in Figure 3, following the same procedure. The two maps are nearly indistinguishable from those included in the manuscript.\nFig. R1. Pointwise t-values of the regression between annual growth fluctuations estimated from tree rings (BAI) and vapour pressure deficit (VPD) using a site-based interpolation. Maps are displaying site-species t-values; a bidimensional interpolation was performed on a spatial resolution of 1 x 1 degree, using the inverse distance weighting method based on the 12 closest neighbours. Interpolations were bounded using boreal mask (all species map)\nAs suggested, we proceeded with a double-checking process and obtained a similar number of convergent sites and similar t-value maps. It can be concluded that the use of climate data extracted at the site scale has no different effect compared to the use of climate data extracted at the nearest grid point for the specific statistical models used in our study.\nThere is some change apparent in Fig. R3 and I would consider these results closer to the truth than what is shown in Fig. 4 in the manuscript, but I would agree that for the most part, the differences are probably not relevant/meaningful. It is interesting, though, that Mean Age notably decreased in importance (as measured by MSE).\nIt is true that the validation dataset yielded slightly modified results. However, these modifications did not lead us to alter our conclusions, as the relative importance of each variable remained fundamentally unchanged. Since site-specific data were utilized in the validation process, slight variations in MAT and MAP values could affect the random forest modeling. Pending the editor's approval, we would prefer to maintain the gridded approach, given the robustness of the results.\nMore broadly, however, I strongly disagree with the authors rational for choosing gridding in the first place, and I believe that using station-based interpolation would be the methodologically more appropriate approach (and that is what I meant to suggest, not interpolation to a very dense, hi-res grid, although I guess I could have been clearer). It is true that weather stations in Canada are sparse, but if one is interested in specific sites, there is no point in interpolating to a regular grid first and then choosing the closest point, when one could just use, e.g. IDW to average available stations in the vicinity. This is even more so the case when lapse-rate adjustment is applied. Every additional step of interpolation introduces more uncertainty.\nInterpolation to a regular and appropriately-spaced grid is only necessary, when area averages are being considered, where representativeness of a station for an area is of concern.\nThe use of gridded data in dendroclimatological studies is a standard practice. This format makes it easier to distribute and share the data with other researchers, as it is structured and standardized. Gridded data often results in a smaller dataset compared to using raw, point-level data. This reduction in data size can lead to more efficient computing, As we have illustrated here, this choice has no impact on the results. In the context of the site-level modeling approach described in this work, the choice of gridded data does not introduce any significant bias or affect the outcomes. We acknowledge that the impact of using gridded data may vary depending on the modeling approach. For instance, if a different modeling approach, such as using Generalized Additive Mixed Models (GAMM), were applied to range-wide species data, the results might be slightly noisier (less precise or more variable) when using gridded data.\nConcerning the correlation between climate variables on the regular grid and at site locations: I am almost surprised that it is not higher. Correlation only measures variability and does not capture constant biases, but the main difference I would expect between site values and regular grid points would be in the lapse-rate adjustment, which primarily manifest in relatively constant deviations. RMSE or bias would have been a better measure here.\n(I am assuming the authors are referring to temporal correlation.)\nSee responses to earlier comments in regard to the use of gridded data. Yes, we were referring to temporal correlation.\nMinor Comments:\nFig. R2: I am not sure what the authors are trying to show... there are clearly areas in the southern part of the domain where station density is much higher than 60 - 100 km and a finer grid would be justified (Fig. R1). At the same time Fig. R2 also clearly shows that there are regions in the Rocky Mountains where the coarser grid misses regionally differing trends that are present at higher resolution, and which would reduce noise in the analysis, if properly captured.\nThere is seemingly greater density, but the number of stations that cover the full period is low. Note too that we don’t have extensive tree-ring data in those areas of higher station density. Hence, there represent little weight in our analyses.\nIf the authors are concerned about some regions having lower levels of precision (more noise) than others due to station density: that is already the case due to strong topographic gradients in the Rocky Mountains, and gridding to a 1 degree rid exacerbates this issue.\nSee responses to earlier comments in regard to the use of gridded data.\nFigure 1: The figure looks significantly better, but the resolution could still be improved.\nWe appreciate your comment and have attempted another method to export the maps. We hope the current version is improved. If it is the PDF saving process that is causing file compression, we can provide separate files as per the editor's suggestion.\nFigures 2 & 3: These figures are much improved and I think the decision to move most of the species-dependent plots into SI was good.\nWhat was the selection criterion for the species shown in the main manuscript?\nThe species featured in the main manuscript hold significant economic importance as they are among the primary species harvested. Additionally, their geographic distribution is particularly threatened by climatic changes. For instance, black spruce, which grows in permafrost soils, is expected to be affected by warming and may experience more frequent insect outbreaks. (Gauthier, et al. 2014)\nWe mentioned these criterion in the methods section:\nl. 424: “Annual fluctuations in growth over time were computed for Picea mariana, Picea glauca, Pinus contorta and Picea engelmanii that are among the primary harvested species in Canada, and which geographic distribution is particularly threatened by climatic changes (Gauthier et al., 2014). Annual growth fluctuations were obtained from the detrending of BAI series that removed growth trends due to tree age and size (Zhang et al. 2018).”\nFigure 6: I am not sure if the change in color map does anything to improve the readability of the map. In particular, the large water bodies on the map are shown in a color that is also part of the color map, and at first glance it was not apparent to me that this was simply a geographic overlay - I think a color that is not part of the color map should be used here, or simply a thin black outline. Furthermore, the authors state that the resolution had been increased, however, I struggle to see that.\nThe figure was revised accordingly with the suggestions.\nI think it would also be appropriate to mention in the context of Fig. 6 that the underlying humidity data were largely estimated from Tmin (using Kimball's method).\nFollowing your advice, we added the method used for VPD calculation in the figure legend :\nl. 704: « Temporal trends in summer VPD (kPa / year) since 1951, estimated from daily temperature and precipitation using Kimball’s method (Kimball et al., 1997).”\nWe appreciate your thorough review of our manuscript and your response to the previous round of reviews. We would also like to express our gratitude for placing our work in perspective and highlighting its relevance."}]

Decision letter