Comparing the Cell Dynamics of Tree-Ring Formation Observed in Microcores and as Predicted by the Vaganov–Shashkin Model

New insights into the intra-annual dynamics of tree-ring formation can improve our understanding of tree-growth response to environmental conditions at high-resolution time scales. Obtaining this information requires, however, a weekly monitoring of wood formation, sampling that is extremely time-intensive and scarcely feasible over vast areas. Estimating the timing of cambial and xylem differentiation by modeling thus represents an interesting alternative for obtaining this important information by other means. Temporal dynamics of cambial divisions can be extracted from the daily tree-ring growth rate computed by the Vaganov–Shashkin (VS) simulation model, assuming that cell production is tightly linked to tree-ring growth. Nonetheless, these predictions have yet to be compared with direct observations of wood development, i.e., via microcoring, over a long time span. We tested the performance of the VS model by comparing the observed and predicted timing of wood formation in black spruce [Picea mariana (Mill.)]. We obtained microcores over 15 years at 5 sites along a latitudinal gradient in Quebec (Canada). The measured variables included cell size and the timing of cell production and differentiation. We calibrated the VS model using daily temperature and precipitation recorded by weather stations located on each site. The predicted and observed timing of cambial and enlarging cells were highly correlated (R 2 = 0.8); nonetheless, we detected a systematic overestimation in the predicted timing of cambial cells, with predictions delayed by 1–20 days compared with observations. The growth rate of cell diameter was correlated with the predicted growth rate assigned to each cambial cell, confirming that cell diameter developmental dynamics have the potential to be inferred by the tree-ring growth curve of the VS model. Model performances decrease substantially in estimating the end of wood formation. The systematic errors suggest that the actual relationships implemented in the model are unable to explain the phenological events in autumn. The mismatch between the observed and predicted timing of wood formation in black spruce within our study area can be reduced by better adapting the VS model to wet sites, a context for which this model has been rarely used.

New insights into the intra-annual dynamics of tree-ring formation can improve our understanding of tree-growth response to environmental conditions at high-resolution time scales. Obtaining this information requires, however, a weekly monitoring of wood formation, sampling that is extremely time-intensive and scarcely feasible over vast areas. Estimating the timing of cambial and xylem differentiation by modeling thus represents an interesting alternative for obtaining this important information by other means. Temporal dynamics of cambial divisions can be extracted from the daily tree-ring growth rate computed by the Vaganov-Shashkin (VS) simulation model, assuming that cell production is tightly linked to tree-ring growth. Nonetheless, these predictions have yet to be compared with direct observations of wood development, i.e., via microcoring, over a long time span. We tested the performance of the VS model by comparing the observed and predicted timing of wood formation in black spruce [Picea mariana (Mill.)]. We obtained microcores over 15 years at 5 sites along a latitudinal gradient in Quebec (Canada). The measured variables included cell size and the timing of cell production and differentiation. We calibrated the VS model using daily temperature and precipitation recorded by weather stations located on each site. The predicted and observed timing of cambial and enlarging cells were highly correlated (R 2 = 0.8); nonetheless, we detected a systematic overestimation in the predicted timing of cambial cells, with predictions delayed by 1-20 days compared with observations. The growth rate of cell diameter was correlated with the predicted growth rate assigned to each cambial cell, confirming that cell diameter developmental dynamics have the potential to be inferred by the treering growth curve of the VS model. Model performances decrease substantially in estimating the end of wood formation. The systematic errors suggest that the actual relationships implemented in the model are unable to explain the phenological events in autumn. The mismatch between the observed and predicted timing of wood formation in

INTRODUCTION
Modeling permits the description of complex biogeochemical processes that occur in nature (Danis et al., 2012), particularly as a suite of factors drive tree-growth response to climate. Tools such as MAIDENiso, TreeRing2000, and the Vaganov-Shashkin (VS) model are mechanistic models for predicting tree growth that account for the endogenous and exogenous factors shaping tree growth and productivity Danis et al., 2012). Among these models, the VS model requires the smallest number of inputs. Furthermore, these inputs include data that are widely used and easily available, such as tree-ring width chronologies and daily mean precipitation and temperature. The availability of these data has led to an increased use of the VS model, which can now be parameterized using a user-friendly interface, the VS-oscilloscope (Shishov et al., 2016) or new MATLAB version of the model (Anchukaitis et al., 2020). Nonetheless, the intra-annual predictions of the VS model continue to lack validation with long-term field observations that, unlike tree-ring chronologies, are scarce.
The main prediction of the VS model is the daily tree-ring growth rate, which is obtained by integrating three partial growth rates based on day-length (photoperiod), temperature, and soil water content. The variation of these environmental factors over the entire year, including the growing season, affects xylem cell production and development to result finally in different wood increments Balducci et al., 2016). The VS model has been applied to simulate long-term climate responses based on the variation of treering width chronologies at a regional scale; study locations include the southeastern United States Evans et al., 2006) and the Tibetan Plateau Yang et al., 2017). Parameterization of the model is necessary to account for the endogenous components of tree-growth response and represents a critical but necessary step of the modeling process, a step that can provide important information about tree growth at the local scale (Evans et al., 2018;Tychkov et al., 2019).
The core of tree-ring growth is cambial activity, and the VS model has been developed to establish not only the start and end of the growing season but also the timing of xylem cell production and the final size of these cells (Popkova et al., 2018). This particular aspect of the VS model has heightened interest in its application with the abundance of literature discussing intra-annual tree-ring dynamics within the cambial zone and xylem; these dynamics are related to both endogenous (developmental patterns and hormones) and exogenous (weather and seasonality) factors (Buttò et al., 2020). The variation of xylem cell traits, i.e., cell diameter and cell wall thickness, provides important information about the trade-off between hydraulic safety and efficiency, a factor that allows plants to adapt and survive in a changing environment (Hacke et al., 2017). Cell traits depend on the temporal dynamics (duration, rate) of their differentiation phases, and environmental factors, i.e., temperature and precipitation, all of which have an influence that varies over the growing season (Fonti et al., 2010;Cuny et al., 2019). By disentangling the effects of environmental factors on tree-ring development at a daily scale, the VS model may have the potential to simulate xylem cell differentiation, even though the existing version of the model focuses mainly on cambial cell production.
The dynamics of cell development-computed via the tree-ring growth rate curves of the VS model-are emergent properties that must be validated via observations of xylem formation. The interpretation of these emergent properties of the models should benefit from the observations (Cook and Pederson, 2011); however, unlike the traditional inter-annual tree-ring widths, long intraannual chronologies of secondary growth are rare and data sets that are limited to a couple of years of observation do not ensure a complete picture of tree-growth response to environmental factors as these responses are often nonlinear (Rossi, 2015). Apart from the study of Popkova et al. (2018), which was based on three years of observation of wood formation, to our knowledge no study has compared the timing of cell development simulated by a VS model with data obtained directly through repeated observations of xylogenesis using microcores.
To validate the VS model predictions with observations, we used our existing 15-year, weekly scaled chronologies of wood formation from across the boreal forest of Quebec, Canada. We aimed to compare the predicted and observed timing of wood growth at both a tree-ring and xylem-cell resolution by using microcores collected from black spruce (Picea mariana Mill.). Black spruce is the dominant species in the Canadian boreal forest and it grows within a great diversity of stand structures, from Alaska to Newfoundland. In Quebec, the distribution of this species extends to 58°N and forms extensive, closed forests in northeastern North America, including some of the wettest and coldest boreal forest stands. The ubiquity of black spruce leads to a very diversified tree-growth response to climate, reflecting the role of various local environmental drivers (Walker and Johnstone, 2014;Nicault et al., 2015).
For model calibration at the tree-ring scale, we relied on the existing literature and field information to assess the performance of the VS model for predicting the variability of factors and parameters that affect black spruce growth. Then, we used observations of xylogenesis to validate the timing of wood formation at the tracheid scale, i.e., cell scale, for which variability is well represented by our long time series of observations.

Study Sites and Tree Selection
Samples were collected from five sites in the coniferous boreal forest of Quebec (Canada) along a latitudinal gradient stretching between 48°N and 53°N ( Table 1). The sites Simoncouche (SIM) and Bernatchez (BER) are located in the balsam fir [Abies balsamea (L.) Mill.)]-white birch (Betula papyrifera Marsh.) bioclimatic domain, while Mistassibi (MIS) and Camp Daniel (DAN) lie in the black spruce-moss bioclimatic domain. Mirage (MIR), the northernmost site, lies in the black spruce-lichen domain and is characterized by a lower tree density and growth than the more southern sites. Mean annual temperature ranges between −3.4 and 1.9°C, with the southernmost and northernmost sites being the warmest and the coldest, respectively ( Table 1). Precipitation ranges from 626 to 906 mm along the latitudinal gradient with drier conditions toward the north ( Table 1). At each site, we selected ten dominant or co-dominant trees, avoiding individuals having polycormic stems, evident parasite damage, reaction wood, or partially dead crowns ( Table 1).

Climate Measurements
Precipitation, temperature, and soil water content at 30 cm soil depth were collected by automatic weather stations equipped with CR10X data loggers (Campbell Scientific Corporation, Canada); these stations were installed in a forest gap within each site. We averaged hourly measurements to obtain daily time series. We filled any minor data gaps caused by short-term technical problems using the ANUSPLIN model (Hutchinson et al., 2009;Hopkinson et al., 2011;McKenney et al., 2011).

Xylem Formation Dynamics
Microcores were collected weekly or fortnightly between April andOctober (2002-2016) from 10 individuals per site. Sampling was performed using a surgical bone needle (2002)(2003)(2004)(2005)(2006)(2007) or Trephor (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016). Sampling at MIR took place from 2012 to 2016. Microcores were dehydrated through successive series of immersions in ethanol and D-limonene. The microcores were embedded in paraffin, cut into 8 µm transversal sections, and stained with cresyl violet acetate (0.16% in water). Xylem cell development was detected by counting the number of cells undergoing each stage of cell differentiation. We counted cambial and xylem cells along three radial lines and identified the differentiation stages of (I) enlargement, (II) cell wall thickening and lignification and (III) mature tracheids. During wood formation, cambial derivatives start dividing and differentiating in xylem and phloem cells by increasing in size. The maturation of the xylem cells is achieved once secondary wall deposition and lignification are completed (Plomion et al., 2001). Cells undergoing different differentiation stages are identified by their different shape, size, color and glistening under polarized light. Cambial cells are irregularly shaped and smaller than xylem cells (Skene, 1969). Due to the different reactions of the components of cell walls to cresyl violet acetate, enlarging cells show pinkish coloration while cells undergoing secondary wall deposition are stained in violet, turning blue when mature (Deslauriers et al., 2003). Mature cells glisten under polarized light.
We estimated the daily sequence of dividing, differentiating, and maturating cells by fitting generalized additive models (GAM) with splines to the number of cells counted for each sampling day, thereby assessing cell production and the timing of division and differentiation of each tree-ring cell at a daily scale (Cuny et al., 2013). Timing represented the day of the year (DOY) when each cell entered into a new developmental stage (Buttò et al., 2019). Accordingly, we estimated the DOY in which each cambial cell stopped dividing and differentiated into an enlarging xylem cell i.e. timing of enlargement, and the DOY in which each cell stopped enlarging and started cell wall deposition i.e. timing of secondary wall deposition and lignification. The duration of enlargement was computed as the difference between the timing of secondary wall deposition and the timing of enlargement for each cell. The period considered as wood formation spanned from the DOY when the number of cambial cells increased in spring to the DOY when the last formed tracheids in the tree-ring was fully mature.

Wood Anatomy
In summer 2017, we collected additional microcores from 10 individuals per site (Rossi et al., 2006). We prepared these microcores using the abovementioned procedure, stained them in safranin (1% in water), and permanently fixed the samples on slides using Permount ™ . We obtained pictures of the transversal sections using a camera fixed on an optical microscope at a magnification of 20×. Radial lumen diameter and cell wall thickness (single wall) were measured for all the study years (15 for all sites except for MIR, for which we had 5 years of observation) using WinCELL (Regent Instruments, Canada). We calculated the tracheidograms of cell diameter, relying on their relative position across the tree ring, and fit our results with GAMs to obtain values representative of each site and year, accounting for cell production as assessed by xylogenesis monitoring (Buttò et al., 2019). We assessed the growth rates of cells as a ratio between cell diameter and the duration of enlargement; the value was scaled between 0 and 1 to be compared with the nondimensional tree-ring growth rate computed by the VS model ( Table 2).

Tree-Ring Time-Series Analysis
We measured tree-ring width using WinCELL (Regent Instruments Inc., Canada) on histological samples obtained from ten microcores collected at all five sample sites. To compare trees with different growth rates, we detrended the chronologies, removing the effects of tree age, genetic growth potential, microsite characteristics, and stand history (Cook et al., 1990). We applied 67% cubic splines with a 1/2 cut-off time-series length to detrend and produce standardized chronologies for the 2002-2016 period using the detrend and chron functions of the dplR package in R (Bunn, 2008).

VS Model Calibration and Validation
We used standardized tree-ring width series to calibrate the VS model for 2002-2016 using the VS-oscilloscope (Shishov et al., 2016). We performed parameterization to avoid any contradictions with field observations and available information for the sites (Tychkov et al., 2019). Accordingly, the calibration of the environmental parameters in the VS oscilloscope, was performed considering the measures of our weather stations during the different moments of the growing season, that we identified by means of the xylogenesisis monitoring. Site-specific parameters were calibrated considering site features, while for parameters linked to black spruce ecology, like root deepness, we relied on literature references to check if the values we established were realistic. Simulations of daily soil water content were improved by activating the soil melting block in the VSoscilloscope, a parameter that was designed originally to include permafrost melting (Shishov et al., 2016). As we were dealing with wet sites, we used this block to consider the marked amount of water released by snowmelt at the start of the growing season. Nonetheless, the current version of VS model does not estimate soil thawing and soil moisture for the dormancy, because these values are assumed to be constant. For each site, we evaluated the effect of the environmental factors on tree-ring growth via graphical interpretation of the partial growth rates patterns provided by the VS model simulations . The graphical representation of the daily average partial growth rates simulated by the VS model allows to determinate the most limiting factor to growth at daily scale, which corresponds to factors linked to the lowest growth rate (Shishov et al., 2016). Indeed, tree-ring growth rate is a function of three partial growth rates that depend on daily temperature, soil water content, and photoperiod; these factors are integrated into the treering growth rate and represent an important result of the model simulations ( Table 2). The start of the growing season as predicted by the model depends on the crossing of a critical threshold for all three partial growth rates; this occurs when each environmental factor activates or permits tree-ring growth . In addition to the minimum temperature for growth to start (T min ), the sum of temperature for growth initiation (T beg ) can be parameterized, in order to consider the forcing needed for cambial resumption. In the current version of VS-oscilloscope, the end of growth occurs when the integral growth rate of the tree ring falls below a critical threshold (critical growth rate Vcr) (Tychkov et al., 2019). We obtained the general pattern of the partial growth rates using the default arguments of the function geom_spline from the R's package "ggformula" (Kaplan and Pruim, 2020).
We computed the cambial cell growth rate, i.e., the average growth rate corresponding to cambial cells produced each year, and then determined their timing of division and enlargement following Popkova et al. (2018) ( Table 2). The protocol proposed by Popkova et al. (2018) extrapolates the timing of cell division and enlargement by the tree-ring growth rate, assuming that the production of a new cambial cell occurs only once the previous cells have left the cambial zone. Cell production would indeed occur without overlapping. The observed number of cells produced each year is necessary to calculate temporal dynamics of cell production, which we obtained through our observations of the microcores.
We assessed the robustness of the VS model simulations using the actual and predicted tree-ring width indices; we relied on the Pearson correlation (R) and root mean square error (RMSE). We retained only models showing a significant R (p < 0.05), producing the smallest RMSE to minimize the variance of simulated indices and to select the simulation with the smallest average prediction error (James et al., 2013). The model validation involved comparing the simulations (named "predicted data") with our results from the microcores (named "observed data"). Pearson correlations and linear regressions served to compare the predicted and observed data. We transformed the data when necessary to meet the assumption of normality, and we fit LOESS functions (span = 0.7) to compare the general patterns visually.

Climate Along the Latitudinal Gradient
For the 2002-2016, our weather stations recorded warmer monthly mean annual temperatures than the long-term average (Tables 1 and 3). For this same period, mean annual precipitation was higher in the north (715 mm, MIR) and lower in the south (622 mm, SIM), with a greater northward interannual and daily intra-annual variability during wood formation (Table 3). However, mean annual and mean daily soil water Growth rate linked to environmental factors; their integration produces the tree-ring growth rate Nondimensional (min 0, max 1)

Tree-ring growth rate
The main output of the VS model and represents the daily tree-ring growth rate Nondimensional (min 0, max 1)

Cambial cell growth rate
Growth rate corresponding to the production of a new cambial cell Nondimensional (min 0,max 1) content decreased from south to north, spanning from 0.17 V/Vs at MIR to 0.36 V/Vs at SIM over the year (Table 3). Daily soil water content during wood formation ranged from 0.17 to 0.42 V/Vs at DAN and BER, respectively.

Duration and Rate of Xylem Growth
For all sites, correlations between the observed and predicted treering width indices were positive and highly significant (SIM and MIS: P < 0.05, N = 15; BER, DAN: P < 0.01, N = 15; MIR: P < 0.01, N = 5) and R-values ranged from 0.54 to 0.90 (Supplementary Figure S1). RMSE ranged between 0.06 (DAN) and 0.1 (MIR), staying below the RMSE threshold of 0.3 and attesting to the good fit between the indices and the simulated results (Supplementary Figure S1). Fixed parameters of temperature for tree growth ranged from 4 to 29°C, showing a 1-5°C difference depending on the site (Supplementary Table S1). Minimum and maximum soil moisture levels were similar along the entire latitudinal gradient, although MIR had the lowest maximum soil moisture (Supplementary Table S1). Partial growth rates showed marked inter-annual variability, in particular for temperature growth rate ( Figure 1). In general, the temperature growth rate peaked at the end of August (DOY 240), while the photoperiod growth rate peaked around the summer solstice (DOY 170) (Figure 1). The maximum temperature growth rate was highest at SIM, where it reached 0.8 relative units, and decreased to 0.62 at the northernmost site, MIR. The water growth rate reached 1, the maximum value, at SIM, MIS, and DAN, around the middle of July (DOY 200) for SIM and MIS and 20 days later at DAN. At MIR, the water growth rate peaked around DOY 240, although it never attained the maximal value ( Figure 1).
From the partial growth rate patterns (Figure 1), wood formation started when temperatures were suitable for the resumption of cambial activity, which according to predictions, occurred between the end of May and the onset of June. Cambial activation was earliest at the southernmost site (SIM, DOY 143) and latest in northernmost site (MIR, DOY 157) ( Figure 2, Table  4). The predicted start of wood formation occurred 4-13 days after the observed start ( Figure 2, Table 4). The linear relationships between the predicted and observed start of the wood formation produced R 2 values that ranged from 0.4 (MIR) to 0.5 (DAN) and showed a systematic overestimation for the start of wood formation (Figure 2). The end of wood formation-occurring when the photoperiod was the most limiting factor-was predicted latest for the southernmost site (SIM, on average DOY 267) and earliest for the northernmost site (MIR, on average DOY 255, Table 4). The predicted end of wood formation was generally overestimated with wood formation lasting 9-19 days longer than observed ( Table 4). The observed end of wood formation showed a higher interannual variation of up to 2 weeks, whereas the predicted end of the growing season varied on average 6 days over the years ( Table 4). The simulated and observed ends of wood formation showed a weak relationship, having R 2 values of only 0.02 (DAN) to 0.3 (MIS) (Figure 2).

Model Parameters and Environmental Properties
Pearson correlations between the predicted and observed daily soil water content ranged between 1 and −1 (Supplementary Table S2). During dormancy, an average 16% of the simulated years produced a strong correlation (R > 0.4) between the predicted and observed daily soil water content (Supplementary Table S2). MIS showed the best performances, where 33% of the years showed a strong and positive correlation between the predictions and observations. We observed the lowest correlations at SIM, where the predicted and observed soil water content correlated strongly in only 7% of the years (Supplementary Table S2). The model simulations for winter resulted in a constant soil water content that was either an underestimate or overestimate depending on the year and site ( Figure 3).
The observed soil water content showed a large intra-annual variation, in particular at the more southern sites (Figure 3). At SIM and BER, soil water content varied in winter between 0.2 and 0.4 V/Vs, culminating at the end of April (DOY 121, Figure  3). From the beginning of May, increased variability in soil water content mirrored a decrease in the observed soil water content, which dropped from 0.5 V/Vs at SIM and 0.6 V/Vs at BER to 0.3 V/Vs at both sites during the first week of October (Figure 3). During the second half of October, soil water content for these sites began to increase slightly and produced a smaller peak in mid-November. The soil water content at MIS and DAN varied little, peaking on DOY 121, although at a much lower intensity than at other sites; MIS and DAN soil water content varied between 0.2 and 0.3 V/Vs and 0.1 V and 0.2 V/Vs, respectively ( Figure 3). In contrast, the predicted soil water content at all sites remained constant until DOY 162, generally underestimated at Wood formation occurred from May to September, although specific dates varied between sites and years (see Table 4 for details). SWC, soil water content.
the southern sites (SIM, BER, and MIS), and overestimated at DAN (Figure 3).
During wood formation, we noted a strong correlation between the observed and predicted soil water content, with 32% of years having a correlation R > 0.4, whereas MIS only showed 7% of the years above this strong correlation threshold (Supplementary Table S2). We observed the best performances at MIR and BER where correlations between the predicted and observed summer soil water content were 60, and 80%, respectively (Supplementary Table S2). In most years, the predicted soil water content was overestimated at DAN and underestimated at SIM and BER (Figure 3). At MIS and MIR, the predicted vs. observed data points fell on the 1:1 intercept, although we also observed a marked variability (Figure 3). FIGURE 1 | Intra-annual partial tree-ring growth rates linked to photoperiod (yellow), temperature (red), and water (blue). General trends for the three partial growth rates are represented by splines. The lowest partial growth rate represents the most limiting factor for each DOY.

Variation in the Timing of Wood Formation Along the Latitudinal Gradient
For all sampled years, the predicted and observed timing of cell division and cell enlargement were highly correlated (R > 0.9, Supplementary Table S2). On average, the general pattern of the predicted timing of cell division showed a constant trend for cell positions 19 to 2 for DAN and MIR, respectively ( Figure 4). Nevertheless, the predicted and observed timing for cambial cells demonstrated a strong linear relationship and R 2 values between 0.8 and 0.9 (Figure 4).
We obtained similar strong correlations and relationships for the predicted and observed timings of cell enlargement (Supplementary Table S2, Figure 5); at the beginning of the growing season, however, the differences between the observed and predicted timings of cell enlargement were initially quite small-from 1 to 2 days depending on the site-and increased consistently with cell position, attaining a difference of 40 DOY ( Figure 5).
The relationship between predicted cambial cell growth rate and observed cell growth rate for all five sites showed R 2 values of 0.4 at MIR, DAN, and SIM to 0.3 at MIS ( Figure 6). The predicted cambial cell and observed cell growth rates were highly correlated for 89% of the years (Supplementary Table  S2). Contrary to what we observed at the other sites, the correlation values at MIR varied considerably and were occasionally strongly negative (2015) (Supplementary Table  S2). On average, the cambial cell growth rate was highest for the first cells of the tree ring, being 0.5 at all sites-except at SIM where it was 0.6-and decreased to 0.1 for the final cells ( Figure  6). The difference between the predicted cambial cell growth rate and the observed cell growth rates ranged between 0.3 (BER) and 0.4 (DAN), having a systematic overestimation of the predicted values. The overestimation of the cell growth rates occurred for the first cells at all the sites in particular and became smaller when cell growth slowed ( Figure 6). Inter-annual variations of the model predictions were similar to those of our observations ( Figure 6). Accordingly, a greater variation of the intra-annual growth rate predicted in MIR matched with a greater variation in observed values ( Figure 6). Regardless, the greater variability in growth rates in MIR could also be linked to the reduced number of available observations for this last site.

DISCUSSION
Understanding the effect of the environmental drivers on black spruce growth is crucial for predicting how environmental change will affect wood formation for this economically and ecologically important species. In Canada, correlations between tree-ring growth and the temperature and precipitation follow two different gradients, unraveling complex patterns in black spruce growth responses to these environmental factors (Nicault et al., 2015). Correlation with temperature is strongest in the North, whereas tree growth is most correlated with precipitation in the Western  Canada (Huang et al., 2010;Walker and Johnstone, 2014).
Correlations of tree growth with temperature and precipitation have also been confirmed by dendro-anatomical analyses at the cellular scale and have demonstrated that environmental factors affect xylem conductivity and long-term patterns of intra-annual cell traits (Puchi et al., 2019). According to the existing literature, the correlation between temperature and tree-ring growth is positive when temperature represents the most limiting factor and is negative in response to warming during the previous growing season and the current spring-probably reflecting drought stress (Huang et al., 2010;Walker and Johnstone, 2014;Nicault et al., 2015;Puchi et al., 2019). Correlation with precipitation, when significant, is always positive (Walker and Johnstone, 2014;Girardin et al., 2016;Puchi et al., 2019). The variation in the local responses to environmental factors might indicate a long-term effect of the environmental factors, which is consistent with the conservative black spruce growth strategy (Chen et al., 2019). In this sense, VS model simulations improvements targeted on black spruce and VS model application on a larger territory could disentangle the common drivers of black spruce adaptation to environmental factors and their effect on treering growth and wood formation at local scale. According to Vaganov et al. (2006), the algorithms underlying the VS model are based on the assumptions that intra-seasonal dynamics are the FIGURE 3 | From the left to the right, average pattern of predicted (red) and observed (blue) soil water content with 95% confidence intervals (left), observed and predicted variations of soil water content during dormancy (middle) and during growing period (right). The white background on the graphs represents the dormancy period, whereas the shaded background represents the period of wood formation; Predicted versus observed soil water content for the dormant and wood formation periods. Sites are organized according to latitude (the southernmost site (SIM) as the top row of graphs, the northernmost (MIR) is on the bottom row).
result of current environments, while the variability in cell production is determined by long-term conditions experienced by the trees, which are integrated to model's simulations by parameterization. The effect of previous years conditions on growth is thus implicit, but the model has never been formulated to split the carryover effect of the environmental factors from the effect of the current environmental conditions. The partial growth rates for the three main environmental factors predicted by the VS model revealed that different limiting factors shape tree-ring production in black spruce. At the beginning of the growing season, which is predicted by VS model for the end of May, temperature is the most limiting factor at all the sites. In cold environments, temperature is the major limiting factor of tree growth. Mean air temperature FIGURE 4 | From the left to the right, average patterns for observed (blue) and predicted (red) timing of cell division in relation to cell position, gray shading represents the 95% confidence interval (left); regression between the predicted and observed timing of cell division (right). Sites are organized according to latitude (the southernmost site (SIM) as the top row of graphs, the northernmost (MIR) is on the bottom row). thresholds during xylogenesis range between 6 and 8°C (Rossi et al., 2007). Accordingly, the average temperature of the lower range of optimal temperatures for our sites (T opt1 , T min ), as estimated by the parameters, was between 6 and 7.6°C (Supplementary Table S1). The predicted end of the growing season occurred around the middle of September when photoperiod became the most limiting factor. In boreal and temperate forests, short days induce cambium dormancy by triggering the physiological mechanisms that lead to growth cessation (Buttò et al., 2020). At high latitudes, the shortening of day length also allows trees to anticipate lower temperatures and induces cold acclimation responses (Wingler, 2015).
The growth rate determined by water, i.e., water growth rate, was never a limiting factor along the gradient-except at BER where the predicted soil water content was heavily underestimated. This result confirms both field and greenhouse observations that the cell production and cell trait sizes of black spruce in eastern Canada do not change significantly with changes in soil moisture (Belien et al., 2012;Balducci et al., 2016). In contrast, Girardin et al. (2016) observed that stands growing FIGURE 5 | From the left to the right, average patterns for observed (blue) and predicted (red) timing of cell enlargement in relation to cell position, gray shading represents the 95% confidence interval (left); regression between the predicted and observed timing of cell enlargement (right). Sites are organized according to latitude (the southernmost site (SIM) as the top row of graphs, the northernmost (MIR) is on the bottom row).
in western Canadian ecoregions are very responsive to variations in soil water availability; this affects tree productivity by reducing the capacity of black spruce to fix carbon (Girardin et al., 2016). Thus, the lack of response of black spruce growth to variations in soil moisture must be ascribed to site conditions, which were well represented by the water growth rate of the VS model.

Site Conditions
The best model performances in terms of prediction of tree-ring width, soil water content, as well as the onset and end of the growing season were obtained at the northernmost site (MIR), one of the driest sites along the gradient. These results confirm that the VS model parameter specification is based mainly on sites FIGURE 6 | From the left to the right, average patterns for observed (blue) and predicted (red) cell growth rates in relation to cell position, gray shading represents the 95% confidence interval (left); regression between the predicted and observed cell growth rate (right). Sites are organized according to latitude (the southernmost site (SIM) as the top row of graphs, the northernmost (MIR) is on the bottom row). characterized by very cold conditions, such as in Siberia or on the Tibetan Plateau (Vaganov et al., 2011;He et al., 2017) and semiarid or arid conditions, such as those in the southwestern US or northern Africa Touchan et al., 2012). Therefore, the VS model ensures better simulations in dry environments, where it has been successfully applied to the study of the long-term response of tree growth to climate Touchan et al., 2012;Zhang et al., 2016;Yang et al., 2017). However, our application of the VS model also provides useful guidance for improving model parameterization and its performance for wetter sites by involving the boreal forest characteristics of soil water balance and soil properties. According to the assumptions of the VS model, soil water content-deemed as the water available for tree growth-depends on precipitation, snowmelt, evaporation, and runoff (Vaganov et al., 2011). In the Canadian boreal forest, however, soil water balance and, in general, water and nutrient distribution within the soil horizons, are correlated strongly with organic layer composition and thickness-boreal soils can reach a depth of 150 cm (Simard et al., 2007;Laamrani et al., 2014). The forest floor of black spruce stands is composed generally of mosses and lichens, which modify the soil thermal conductivity and water content (O'Donnell et al., 2009). Due to its thickness and speciesspecific composition, the forest floor promotes water retention in the surface soil, thereby acting as a thermal buffer between the atmosphere and soil to decrease thermal conductivity (Sharratt, 1997;Turetsky et al., 2010). Accumulation on the forest floor entails a drop in soil temperature, and as such the decomposition rate of the organic matter slows; these conditions result in decreased tree growth Lavoie et al., 2005).
All soils at our study sites are podzols, although they differ in terms of depth to bedrock, soil thickness, and species composition on the forest floor. These peculiarities likely explain the mismatch between precipitation and soil water content that we observed along the gradient between 2002 and 2016. Despite receiving less precipitation, the southernmost sites of BER and SIM remained the wettest of the gradient. SIM is characterized by a thin organic layer (10-20 cm thick) and a very shallow bedrock . As mentioned previously, the soil organic layer can negatively affect tree growth in boreal forest stands. However, the thickness of the organic layer at SIM was too thin to affect tree productivity (Lavoie et al., 2007); thus, the site at SIM is particularly favorable to black spruce growth. The other sites had deeper soils, 20-40 cm in depth, and differed in their forest floor characteristics . The presence of mosses and a deeper soil at MIS and DAN foster a greater soil retention in the upper layer of the soil. In contrast, MIR has a forest floor dominated by lichens and contains a very thin soil depth; these properties favor the establishment of a thinner and drier forest floor (Waldron et al., 2013).
In this study, predictions of soil water content have improved via the activation of the "soil melting block" in the VS model, a specific module of the model that considers the contribution of snowmelt to soil moisture at the beginning of the growing season (Shishov et al., 2016). Rossi et al. (2011) had previously highlighted the importance of snowmelt at our sites by observing the effects of snowmelt on the duration of xylogenesis and cell production. The presence of accumulated snow prevents soil warming in the spring, and these colder temperatures inhibit root reactivation, delaying spring rehydration and, in turn, cambial reactivation (Turcotte et al., 2009). A later snowmelt corresponds to delays in soil warming, cambium activation, and cell production (Vaganov et al., 1999). Thus, the date of snowmelt influences the timing of wood formation and may explain some of the difference between the observed and the predicted onset of the growing season.

Timing of Cell Division and Differentiation
Compared to the observed data, the predicted timings of cell division and enlargement were delayed from a few days to two weeks. This overestimation was more constant for the timings of cell division, and gradually higher across the tree ring for the timings of cell enlargement. For this, we advance the hypothesis that the upward shift of the predicted timings of cell division relative to observations can be improved by a more precise prediction of the start of the growing season. The increased mismatch between the predicted and observed timing of enlargement can be reduced further by including some information. Popkova et al. (2018) illustrate that the time that dividing cells spend in the cambial zone can be determined by splitting the predicted growing season into a number of periods equal to cell production. To leave the cambium, each cell must attain an average growth rate computed by the ratio between the cumulative sum of the daily growth rate obtained by VS model and cell production. However, to obtain a more reliable reconstruction of the temporal dynamics of cell division, a new assumption can be added to the model, one that considers the natural variation of cell production rates over the growing season. Furthermore, an accurate estimate of the number and order of cells entering the xylem should account for the number of cells that will never differentiate, i.e., the number of dormant cambial cells that constitute the mother cells and the cambial cells that undergo phloem differentiation (Plomion et al., 2001).
Cell production rates, depending by cambial activity, interact with the residence time of each cell within the differentiation zones. Variations in these cell production rates stem from the complex seasonal patterns that occur during tree-ring formation (Cuny et al., 2013;Balducci et al., 2016). An improved performance of the VS model timing procedure can be achieved by considering the specific pattern of the cell production seasonal rate. Cell production rate is characterized by a bell-shaped curve, with the peak and following decrease being in response to intraseasonal dynamics and internal factors (Deslauriers and Morin, 2005). The peak in the cambial growth rate i.e. the number of cells produced per time units, matches a peak of auxin, the main hormone promoting cell division, and occurs during earlywood formation (Buttò et al., 2020). At the same time, soil water content and temperature affect cambial growth rate and cell production, the first prevailing in dry environments and the second in boreal climates (Deslauriers et al., 2008;Ziaco et al., 2018;Ren et al., 2019).

VS Tree-Ring Growth Rate and Cell Differentiation
We compared cell growth rates with the cambial cell growth rates predicted by the protocol of Popkova et al. (2018). In particular, we tested the possibility of inferring the cell growth rate from the tree-ring growth curve obtained by the VS model by adjusting the procedure of Popkova et al. (2018) by focusing on the differentiation zones rather than the cambial zone. Once cell production is known, Popkova et al. (2018) infer the cambial cell growth rate using the curve of tree-ring growth rate and perform the most suitable regression between cambial cell growth rates and the measured cell diameters from actual tracheidograms. With the coefficients from the abovementioned regression, Popkova et al. (2018) then obtain those parameters useful for producing synthetic tracheidograms. However, cell diameter can be computed without using observations by integrating the relationship between cell traits and their temporal dynamics of differentiation (Cuny et al., 2014;Buttò et al., 2019). The timings of enlargement computed by the protocol proposed by Popkova et al. (2018), identify the start of cell differentiation i.e. the timings of cell enlargement, which is the day in which a cell starts increasing in size ( Table 2). Xylem cell differentiation temporal dynamics assessment requires the computation of the timings of secondary wall deposition and cell maturation to estimate the duration of cell enlargement and the duration of secondary wall deposition. The duration of these differentiation phases depends on exogenous factors i.e. photoperiod, soil water content and temperature, which are already included in the model, and on endogenous factors i.e. hormonal signaling, sugars concentration (Rathgeber, 2017;Cartenì et al., 2018;Buttò et al., 2020). Then, a new framework that considers the nonlinear relationship between the duration of cell enlargement and cell diameter can be implemented, which also includes the effect of the secondary wall on the process of cell expansion (Buttò et al., 2019). During cell differentiation, secondary wall deposition influences cell diameter by altering the plastic properties of cell walls. Cell wall thickening, which occurs during cell trait differentiation, constrains cell expansion, especially during latewood formation when cell deposition lasts longer than during earlywood formation (Cartenì et al., 2018). These new findings that assess the quantitative relationships between cell traits, developmental temporal dynamics, and environmental factors should be implemented into the next versions of the VS model to improve simulations (Ziaco, 2020).

The End of the Growing Season
The end of growth was one of the most critical portions for our comparisons, and in a wider sense, for studies of both xylem and bud phenology (Gallinat et al., 2015). At our sites, the end of wood formation occurred latest at the southernmost sites and showed a clear latitudinal pattern. Among the climatic factors, photoperiod and temperature are the main known drivers of the end of the wood formation, and the end of growth occurred earliest at the northernmost sites (Rossi et al., 2014;Cuny et al., 2019). However, other factors may also affect the end of the growing season, including cell production and sugar availability, the latter being crucial for cell wall deposition. Cell wall thickening at the end of the growing season lasts 30-40 days after the last cell is differentiated and drives the timing of latewood formation Cartenì et al., 2018). Cell production, which is the result of cambial activity, also influences the kinetics of wood formation, since competition for resources increases with the number of cells, and sugar availability shapes cell traits by modifying their differentiation phases (Cartenì et al., 2018). The end of the growing season can be parameterized during VS model calibration thought the critical growth rate, a coefficient representing the threshold under which the integral growth rate of the VS model does not allow wood formation anymore. Being based on the integration of the three partial growth rates linked to the environmental factors, the critical growth rate may not be sufficient for establishing a realistic end of the growing season, a moment that is not highly dependent on weather conditions, unlike at the start of the growing season .

CONCLUSIONS
In this study, we compared the predictions of intra-annual treering formation dynamics estimated by the Vaganov-Shashkin (VS) model with field observations based on a 15 years-long monitoring of xylogenesis across a latitudinal gradient. Our results show that the model successfully describes the influence of climate variables on black spruce tree-ring growth. However, algorithms can still be improved to ensure a more reliable estimate of the timing of wood formation. Considering that the model has been tested for the first time on the Canadian boreal forest, it is difficult to relate the mismatch between predictions and observations to issues linked to the use of the model on wet environments or to the underlying algorithms computing the timings of cell development. VS-model parameterization is indeed based on semi-dry, dry or permafrost environments, which better support the assumptions of a constant soil water content during dormancy and of a closer relationship between tree growth and water availability at growth resumption. In our case, we tested the model on wet sites, and observed a low representativeness of the soil water content in spring, leading to cumulative errors in soil water estimation during the growing season. A recalibration of the model on other wet sites could provide more robust elements to improve the performance in terms of predicted water content and partial growth rate. Moreover, the algorithms computing the timings of cell production should consider that the cambial growth rate varies during the growing season with consequent effects on the residence time of each cell within the cambial zone. The development of a new procedure must involve phenology and dynamics of each cell during development, based on the underlying biological processes. As a consequence of cell differentiation dynamics, the computation of the cell traits, i.e. cell diameter and cell wall thickness, should be included in form of new module of the VS model based on larger datasets, encompassing species from different environments to obtain more extensive and exhaustive growth simulations.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
VB, SR, and VS conceived and planned the study. VB carried out the analysis and wrote the manuscript. MP provided the initial version of the script to compute the timing of wood formation with the VS model. MP, IT, and VS provided technical support for the analyses. VB, VS, MP, IT, SR, AD, MH, and HM contributed to the interpretation of the results. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
Special thanks are extended to D. McKenney and P. Papadopol for sharing their data set with the temperature and precipitation chronologies and to Murray Hay for verifying the English of the text.