Decadal Variations in Eastern Canada’s Taiga Wood Biomass Production Forced by Ocean-Atmosphere Interactions

Across Eastern Canada (EC), taiga forests represent an important carbon reservoir, but the extent to which climate variability affects this ecosystem over decades remains uncertain. Here, we analyze an extensive network of black spruce (Picea mariana Mill.) ring width and wood density measurements and provide new evidence that wood biomass production is influenced by large-scale, internal ocean-atmosphere processes. We show that while black spruce wood biomass production is primarily governed by growing season temperatures, the Atlantic ocean conveys heat from the subtropics and influences the decadal persistence in taiga forests productivity. Indeed, we argue that 20–30 years periodicities in Sea Surface Temperatures (SSTs) as part of the the Atlantic Multi-decadal Oscillation (AMO) directly influence heat transfers to adjacent lands. Winter atmospheric conditions associated with the North Atlantic Oscillation (NAO) might also impact EC’s taiga forests, albeit indirectly, through its effect on SSTs and sea ice conditions in surrounding seas. Our work emphasizes that taiga forests would benefit from the combined effects of a warmer atmosphere and stronger ocean-to-land heat transfers, whereas a weakening of these transfers could cancel out, for decades or longer, the positive effects of climate change on Eastern Canada’s largest ecosystem.


Results
Wood biomass production (MW) in EC's taiga is primarily controlled by growing season temperatures (Fig. 2). The correlation coefficients between regionally averaged MW and May to August temperatures are all positive (all r > 0.3, p < 0.05), indicating that warm springs/summers positively influence wood biomass production. The spatial coherency of the response to temperatures across the investigated area is another indication of the influence of growing season temperatures on black spruce forest productivity. Indeed, the correlations with May to August temperatures are mostly positive (Fig. 3), irrespective of the site's position in EC (no clear pattern emerges, however). Consequently, the regionally averaged MW is imprinted by a strong and significant regional surface temperature signal. During the growing season, MW correlates well with surface temperatures over much of EC, from the western margin of the Hudson Bay to the Labrador Strait (Fig. 4a). Significant correlations with MW also extend eastward and reach the western portion of the Atlantic ocean (about 1000 km East of Newfoundland, Canada). Our analysis only reveals a weak influence of precipitation on MW both at the site scale and the regional scale.
Not all components used to calculate MW respond equally to temperature variability during the growing season ( Fig. 4b-e). Earlywood densities (EW D ) and latewood widths (LW W ) are both very poorly correlated with growing season temperatures or precipitation Fig. 3. Conversely, earlywood widths (EW W ) and latewood densities (LW D ) both respond similarly to temperatures and precipitation over the 1901-2000 period, although with different strengths. Compared to MW, EW W presents slightly weaker relationships with temperatures and tends to be better associated with spring temperatures (Figs 2, 3 and 4b). However, LW D appears to be very responsive to growing season temperatures (r = 0.69, p < 0.01), even more than MW or EW W . In fact, LW D presents the strongest and most spatially coherent response of all proxies forming MW. The strength of this response is mostly attributable to a clear and highly significant positive effect of August temperatures on LW D (Fig. 2), visible across Quebec-Labrador (Fig. 3). Indeed, the site-specific correlation coefficients are all systematically higher Scientific RepoRts | 7: 2457 | DOI: 10.1038/s41598-017-02580-9 than 0.55 (p < 0.01) for August temperatures. The average LW D also carries a strong regional surface temperature signal, even stronger than the one imprinted into the MW series (Fig. 4e). The regional surface temperature signal resonates over much of the North Atlantic ocean, where temperatures offshore of Newfoundland appear to be significantly correlated with LW D .
Growing season temperatures that control MW and its components across EC are positively associated with variations in SST across the Atlantic. The summer AMO May-Aug index that summarizes variations in Atlantic SSTs during the warm season correlates well with the growing season temperatures over EC (r = 0.5, p < 0.05) during the twentieth century (Fig. 5a). In fact, in the circum-Atlantic area, EC represents the sector where growing season surface temperatures are most strongly related to SST variability, as described by the summer AMO index. Regions of EC correlated with summer AMO variability include both sides of the Hudson Bay and extend eastward to the coast of the Labrador Sea, where maximal associations between surface temperatures and summer AMO are found. No correlation could be found between AMO May-Aug and growing season precipitation variability in the study area. Conversely, the NAO pattern has a strong influence on winter (Dec-Feb) temperatures in the circum-Atlantic area, generating colder conditions in EC during periods of positive winter NAO (Fig. 5b). However, the winter NAO does not appear to be significantly related to growing season temperatures in EC (Fig. 5c).   Because EC's surface temperatures are highly responsive to the AMO, biomass production in this region is also affected by SST variability. Indeed the correlation between summer AMO and MW is 0.34 (p < 0.05) over the last century. Interestingly, a cross-spectral coherency analysis also highlights a fairly high (>0.7) and significant (<0.1) coherency for periodicities between 20 and 30 years during the same period (Fig. 6a). Thus, our analysis suggests that MW and AMO May-Aug might share a comparable decadal variability component, despite the fact that the coherency at the inter-annual scale remains modest. The coherency between spectra is even more visible when lower frequencies (f > 30 years) are filtered out from both MW and the AMO series (Fig. 6c). Ultimately, isolating frequency bands comprised between 20 and 30 years (Fig. 6d) allows magnification of the spectral coherency. Moreover, a lagged correlation analysis reveals that filtered (20-30 years) variations in summer AMO systematically precede those of MW by approximately 4-5 years within the same frequency band (Fig. 6e). Interestingly, winter NAO also shares coherent 20-30 year periodicities with both MW (Fig. 6f-j) and summer AMO series ( Fig. 6k-o) and could be identified as a potential source of decadal variability and persistence in black spruce growth at the centennial scale. Ocean-atmosphere interactions that occur in the North Atlantic region significantly affect wood biomass production over decadal to inter-decadal time scales, resulting in spatially coherent responses of black spruce forests to temperatures across EC. The 20-30 years periodicities common to MW, summer AMO and winter NAO time series suggest that significant amounts of heat that contribute to the growth of EC's taiga forests could originate from the North Atlantic Ocean. However, numerous processes evolving at different spatial and temporal scales may control heat transfers towards EC, all of which may be imprinted in the AMO series. In order to highlight these processes, we determined the three main modes of SST variability in the North Atlantic through a principal component (PC) analysis (Fig. 7) and related them to AMO/NAO variability. Indeed, the first PC (PC1, correlation with AMO = 0.65, p < 0.01) represents 31% of the overall North-Atlantic SST variability and clearly summarizes the dominant, basin-wide, low-frequency summer AMO variability (Fig. 7a). Fluctuations in PC1 result in generalized warmer surface air temperatures across the subpolar North-Atlantic ( Fig. 7a) with warm anomalies forming a horseshoe pattern, reaching northern Europe and then extending westward to the coast of Labrador. PC1 presents a dominant low frequency component with most spectral power around 50 to 70 years (Fig. 7b). The pattern of correlation of SSTs with PC1 is quite similar to the one found between surface air temperature and AMO (Fig. 5a). The second component (PC2, 11% of SST variability, correlation with AMO = 0.42, p < 0.01), represents decadal to inter-decadal fluctuations of SSTs over the subpolar region of the North Atlantic ocean. The second PC highlights a particularly strong response of surface temperatures over the western portion of the North-Atlantic basin, south of Greenland. The imprint of PC2 on surface air temperatures across the North-Atlantic basin extends westwards to reach the Quebec-Labrador peninsula in EC (Fig. 7c). Interestingly, PC2 has most of its spectral power between 20 and 30 years (Fig. 7d). Lastly, PC3 (9% of SST variability, correlation with AMO = 0.02, p > 0.05) also shows considerable spectral power between 15 and 25 years and reveals a pattern of correlation ( Fig. 7e) that closely resembles the one found for NAO and surface air temperatures (Fig. 5b,c).

Discussion
This study highlights the critical importance of growing season temperatures in determining the pace at which black spruce boreal forests produce and store wood biomass in tree trunks. The evidence for a widespread and spatially coherent dependence of wood biomass production on summer temperatures is based on the investigation of EC's most extensive tree ring multi-proxy network ever assembled in the unmanaged boreal forest of Quebec-Labrador. The sensitivity of high-latitude species to growing season temperatures is a well-known feature of dendroclimatological studies. For this reason, boreal tree ring series have commonly been used as proxies for past growing season temperatures both in the Canadian [19][20][21] and Eurasian context [22][23][24][25] . In the study area, growing seasons are most often short, starting in May and ending in late August-early September 26 . Warmer conditions, particularly during growth initiation and growth completion, increase the photosynthetic efficiency and thus synthesis/allocation of carbonaceous materials to tree trunks 27,28 .
Previous studies that investigated the response of black spruce forests relied strictly on ring width chronologies [29][30][31] . Based on those measurements, important heterogeneities were found in the relationship between tree ring series and climate across the full boreal biome 32 . In the taiga forest of EC, significant and positive dependence on summer temperatures was restricted to the northernmost populations (>54 °N) 29 . South of 52 °N, temperature have had a negative impact on growth and induced widespread growth declines driven by an increase in both evaporative demand and autotrophic respiration 31 . Further south, in the closed-crown boreal forest, black spruce forest stands shift from precipitation-dependant to thermo-dependant populations around the 49 th parallel 30 . Our analysis reveals a considerably more coherent response to temperatures across the entire taiga because MW is capable of tracking temperature controls throughout the growing season. Indeed, the proxy used to calculate wood biomass production (MW) integrates variations in ring width and wood density, during both the earlywood and latewood growth phases. For example, our results revealed that two proxies involved in the calculation of MW are strongly thermo-dependent: earlywood widths (EW W ) and latewood densities (LW D ). The thermo-dependence of EW W implies that early spring temperatures play a crucial role in the production of wood biomass. Indeed, a detailed investigation of black spruce wood micro-cores revealed that warmer (earlier) springs favor an enhanced growth 33 , mainly because growth initiation occurs sooner and earlywood xylem cells grow faster under such conditions. The response of MW to temperatures also appears to be very coherent because LW D reflects the strong influence of late summer conditions on latewood formation. Relationships between LW D and late summer temperatures are well established in the literature: thick-walled latewood xylem cells form during warm (and perhaps also dry) summers 34 , therefore augmenting wood biomass produced per unit surface (volume). LW D has a notable influence on the strength and coherency of the response of MW to summer temperatures, but the overall effect on wood biomass is dampened by LW W (eqs 5 and 6), a very complacent proxy that represents less than a third of the full ring width (31%).
Our study revealed that EC's growing season temperatures are among the most responsive to AMO variability across the circum-Atlantic area (Fig. 5a), especially at decadal to inter-decadal time scales. The spatial correlation with both growing season temperatures and MW are unaltered when periodicities longer than 30 years are filtered out from the summer AMO series (not shown), confirming the relevance of decadal to inter-decadal timescales for the study of processes controlling heat fluxes in the western sub-polar North Atlantic Ocean. Based on previous works 35, 36 , we hypothesize that climatic and ecological impacts triggered by summer AMO variability reflect the interactions between two dominant modes of summer SST variability and their respective impacts on surface temperature patterns across the Atlantic basin. For example, from an analysis of observational records (both marine and continental) and GCM modeling results, Frankcombe et al. 36 have highlighted decadal to inter-decadal (20-30 years long) fluctuations in SSTs and subsurface temperatures at subtropical/subpolar latitudes in the North-West Atlantic Ocean. These fluctuations could correspond to those exemplified by our PC2. In fact, 20-30 years fluctuations in summer SST patterns have been attributed to the westward (zonal) migration of surface and subsurface temperature anomalies. The mechanisms invoked to explain the westward migration of temperature anomalies across the Atlantic basin involve the development of thermal Rossby waves forced by antecedent temperature and salinity gradients. Such processes have been theorized in ocean-only models 37 and represent an important component of SST variability 38 . Such 20-30 year periodicities are clearly discernible in SST and subsurface temperature series 35 , GCM simulations [39][40][41] and proxy-based reconstructions of the strength of the Atlantic Meriodional Overturing Circulation (AMOC) 42,43 which has commonly been invoked as an important driving force affecting the AMO 42,44,45 . Overimposed on decadal to inter-decadal SST variability, Frankcombe et al. 36 have highlighted an additional source of multi-decadal variability in the North Atlantic visible on frequency bands lower than 50 years and most typically associated with the dominant pattern of the AMO. It also presents a characteristic horseshoe pattern much similar to the pattern found by PC1 in the present study. The physical mechanisms causing these basin-wide, low-frequency SST variations in the North Atlantic however remain debated 13 . Frankcombe et al. 36 argued that they most likely result from the modulation of the decadal component of SST variability by meridional temperature and salinity exchanges between North-Atlantic and Arctic waters, further affecting ocean circulation at the scale of the basin. However, from an analysis of slab-ocean models, Clement et al. 13 found that those variations can be triggered by atmospheric circulation only, without any intervention from ocean circulation.
Here, we argue that the North-Atlantic ocean can store the heat affecting decadal to inter-decadal changes in surface air temperatures of surrounding continental masses. However, other processes need to be invoked to explain the redistribution of heat over EC. Interactions with sea ice cover are particularly relevant in that context. This is particularly true for EC, which is surrounded by the Hudson Bay and Labrador Sea, two oceanic systems that undergo significant seasonal freezing during winter. At the scale of the circum-Atlantic region, warmer SSTs contribute to reducing sea ice concentrations in both the Labrador Sea and Hudson Bay during winter. The link between North Atlantic SSTs and sea ice concentrations is supported by state-of-the-art 1000 year modeling experiments, which confirm the negative impact of summer SSTs on the ice cover of these subpolar seas 41 . Sea ice cover has a well-known impact on temperatures at subpolar latitudes, but the effect appears to be strongest in late fall and winter 46,47 . However, feedbacks from ice cover on atmospheric conditions may also extend up to the warm season. Indeed, recent research has demonstrated that decadal to inter-decadal variations in sea ice concentrations covary with both temperatures and growth of black spruce forests in EC, south of 52 °N 31 .
Whereas the North Atlantic stores the heat required to stimulate EC's taiga productivity, NAO-driven winter climate variability could contribute to initiate, amplify or increase the persistency of heat transfers from oceans to adjacent lands. Our study has shown that, just like North Atlantic SSTs, the NAO shares coherent decadal to inter-decadal periodicities with both MW (Fig. 6f-j) and summer AMO variability (Fig. 6k-o). More precisely, surface air temperature response to PC3 across the Atlantic basin suggests that the winter NAO imprints surface air temperature patterns across the North Atlantic. Surprisingly, common periodicities between MW, summer AMO and winter NAO exist despite the fact that NAO variability remains mostly a winter phenomena, with virtually no impact on summer temperatures (Fig. 5c). Hence, NAO-driven climate variability cannot be invoked as a potential factor that directly stimulates wood biomass production in the study area. Nevertheless, recent research Scientific RepoRts | 7: 2457 | DOI:10.1038/s41598-017-02580-9 has shown that NAO and AMO-type variability interact on decadal to inter-decadal time scales, although some debate remains as to whether the former is a cause or a consequence of the latter 48 . NAO variability can "excite" the ocean oscillator, drive and amplify changes in the vigor of the AMOC 39 which in turn influences AMO-type patterns of SST variability in the North Atlantic. For example, positive NAO conditions resulting from large surface pressure differences between Iceland and the Azores lead to the strengthening of westerly winds. Stronger westerlies increase surface turbulent heat fluxes over the Labrador Sea and the subpolar gyre. This results in an intensification of deep water formation in the Labrador Sea 40 , in the stimulation of the oceanic circulation and in significant sea ice loss in subpolar Atlantic. Contrastingly, it is well known that the noisier atmosphere has a limited intrinsic capacity to generate, by itself, decadal climatic persistence. Thus, the NAO might, in turn, be imprinted by changes in SST patterns resulting from positive ocean-atmosphere feedbacks 48 . Indeed, warming SSTs during positive summer AMO tend to favor the development of NAO-like pressure gradients over the North Atlantic. The development of such pressure gradients contributes to maintain and intensify the NAO over decadal to inter-decadal time scales. Here, we put forward the idea that winter-based NAO atmospheric variability might modulate SST patterns in the North Atlantic sector, but ultimately, the timescale at which heat is transferred to EC is more likely to be set by the speed at which SST anomalies travel westward, as part of the AMO variability. Further investigation of this ocean-atmosphere feedback mechanism needs to be undertaken by means of coupled GCMs modeling experiments to fully explore its causes and possible consequences on EC's climate and ecosystems. A promising framework (the delayed oscillator model) has recently been proposed by Sun et al. 49 and could be used to investigate these complex ocean-atmosphere interactions and their impacts on climate variability at the circum-Atlantic scale.
One question remains to be answered: why is the response of black spruce forests to North Atlantic to SST variability lagged by 4-5 years? The delay could have a biological origin and result from a preferential allocation and re-utilization of photosynthetic products stored by trees during years of favorable growth. Even if we took good care of removing the autocorrelation component in our series, it is difficult, for now, to completely rule out this effect without a deeper investigation of high-resolution dendroecological data and associated ecophysiological modeling. Another, perhaps more realistic explanation for this delay would be that ocean-atmosphere interactions that take place in the North-Atlantic sector may take up a few years to amplify to a level that significantly impacts wood biomass production across EC's taiga. From an ocean-only perspective, several works have confirmed that AMO variability leads changes in air temperature by a few years (3 to 5 years in average) and have associated this characteristic feature to the propagation of a thermal Rossby mode across the sub-polar North-Atlantic 37,50,51 . Delays could also be expected with respect to interactions between NAO, AMO and ice cover dynamics at decadal to inter-decadal timescales. For example sea ice anomalies can imprint summer North Atlantic SSTs in a way that increases the persistence of severe sea ice conditions during following years 52 , therefore slowing down the pace at which North-Atlantic SST anomalies may propagate to adjacent lands.
Our findings clarify the response of EC's largest ecosystem to long-term temperature variability. Based on our analysis, warming temperatures anticipated for the 21st century 53 are expected benefit growth of EC's taiga forests. However, increasing temperatures will have a positive effect only up to a certain limit. In the southernmost part of the boreal forest of eastern Canada, excess heat was associated with greater temperature-induced water stress, resulting in large-scale growth declines 31 . In the present study, no such growth decline was observed in the MW series, probably because average temperatures are considerably colder in the taiga ecosystem than in the southern boreal region, making temperature-induced growth declines less likely in the unlogged part of the boreal forest. Our study also implied that thermo-dependent biological processes controlling black spruce growth are linked to the state of the ocean and heat transfers from lower latitudes to subarctic environments. Thus, a possible lack of phase between anthropogenic warming trends and AMO-like patterns of SST variability in the North Atlantic could attenuate the positive effects of climate change on EC's taiga forests. Indeed, proxy-based reconstructions and state-of-the-art modeling studies have pointed to a general slowdown of the AMOC since the beginning of the 20 th century, -a slowdown that should continue at least during the first part of the 21 st century 43 . Hypothesizing that AMOC controls part of the AMO (this common belief has however been questioned recently 13 ), thus, a slowdown of the overturning circulation associated with colder conditions in the North-West Atlantic could cancel out, for decades or longer, the beneficial effects of human-induced warming on EC's largest ecosystem.

Methods
Tree Ring Series. The tree ring network used in this study is composed of 48 sites spread out across the taiga forest of EC, north of 52 °N (Fig. 1). At each site, stems were selected using strict criteria of tree age, topographic homogeneity and structure. Only dominant trees with a symmetrical shape and free from major anomalies were selected. Collected cross-sections were dried and finely sanded. Tree rings were cross-dated and then measured along two or three radii using a micrometer with an accuracy of 0.001 mm (Velmex Inc., Bloomfield, NY) under a binocular magnifying glass. The accurate dating of each ring was then verified with the program COFECHA 54 . Density measurements were performed at each site following established procedures 55 . Only discs without any anomalies (reaction wood, branches, rotten wood, and so forth) were selected. Three wood samples per tree were precisely cut in 1 mm probes and then placed in a Soxhlet apparatus with ethanol for resin extraction. After a certain acclimation time to ensure constant size and hygroscopic conditions, the probes were X-rayed. To separate latewood from earlywood measurements, the X-ray micrograph was analyzed on a DENDRO 2003 microdensitometer (Walesch, Switzerland). A cellulose acetate calibration wedge 56 was used to convert the lightness of measurement into density (g/cm3) values. This densitometry analysis provided density and width measurements for both earlywood (EW D and EW W ) and latewood (LW D and LW W ). Cross-dating statistics are presented in Table 1 for each measured proxy. Wood Biomass Index (MW). In the present study, tree ring width and density measurements (for both latewood and earlywood) were used to calculate the annual biomass of wood produced by black spruce forests at each site investigated. Basal area increments (BAIs) were calculated from tree ring widths as follows (assuming a circular cross-section): where R t and R t−1 are the stem radial increments at the end and beginning of a given annual ring, respectively. For thin slices of wood, it is straightforward to show that BAI (cm 2 ) is directly proportional to basal volume increment (BVI) (cm 3 ). Therefore: t t BVI measurements can therefore be decomposed in earlywood and latewood fractions: One can finally convert the respective ring volumes (cm 3 ) into mass units (g) using density measurements (D, in g cm −3 ). Therefore, the mass of wood (MW, in g yr −1 ) can be obtained from the dimensionally correct equation: MW series were standardized using 60 years spline functions, and the first-order autocorrelation was removed from the series. We acknowledge that tree ring standardisation and removal of autocorrelation in our tree ring series inevitably censors its lower frequency component (>50-60 years) and complicates the interpretations of teleconnection phenomena that dominate on those time scales. Ideally, to investigate those teleconnections, standardisation should be conducted based on the Regional Curve Standardisation method that was conceived precisely to overcome this problem 57 . However, it cannot be applied in our case because we do not meet the specific criteria for that method 58 (eg. use of trees of comparable cambial ages but that lived during different periods).
Response of MW to climate. Annual biomass increments (MW) and other density and ring width parameters at all sites were correlated to gridded land temperature and precipitation datasets over a common period: 1901-2000 (N = 100). For this purpose, we used gridded surface air temperatures (at 2 m), precipitation and sea surface temperatures from the ERA reanalysis project (1°) 59 . MW values averaged across EC were also compared to regional temperature and precipitation averages based on the same method. Analyses of teleconnections between MW and NAO 60 and AMO 61 over the last century  were based on standard correlation analysis and Morlet cross-spectral coherency analysis. In order to guard against the risk of spurious spectral coherencies, 90% significance levels were calculated based on the generation of 300 Monte Carlo AR(1) randomizations. Correlations were calculated to examine the mechanisms responsible for the teleconnection. Effective degrees of freedom were calculated in order to adjust p-values for the presence of temporal autocorrelation in the series 62 .