# High-Resolution X-Ray Computed Tomography: A New Workflow for the Analysis of Xylogenesis and Intra-Seasonal Wood Biomass Production

Romain Lehnebach, Matteo Campioli, Jozica Gričar, Peter Prislan, Bertold Mariën, Hans Beeckman, Jan van den Bulcke

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Romain Lehnebach, Matteo Campioli, Jozica Gričar, Peter Prislan, Bertold Mariën, et al., High-Resolution X-Ray Computed Tomography: A New Workflow for the Analysis of Xylogenesis and Intra-Seasonal Wood Biomass Production. Frontiers in Plant Science, 2021, 12, 10.3389/fpls.2021.698640. [https://doi.org/10.33537/pps.2021.698640](https://doi.org/10.33537/pps.2021.698640).

HAL Id: hal-03353742

https://hal.inrae.fr/hal-03353742

Submitted on 24 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Check for
updates

# High-Resolution X-Ray Computed Tomography: A New Workflow for the Analysis of Xylogenesis and Intra-Seasonal Wood Biomass Production

Romain Lehnebach¹,²*, Matteo Campioli³, Jozica Gričar⁴, Peter Prislan⁴, Bertold Mariën³, Hans Beeckman⁵ and Jan Van den Bulcke¹

¹ UGCT-UGent-Woodlab, Laboratory of Wood Technology, Department of Environment, Faculty of Bioscience Engineering, Ghent University, Gent, Belgium, ² AMAP Laboratory (botany and bio-informatics of plant architecture and vegetation), Université Montpellier, CIRAD, CNRS, INRAE, IRD, Montpellier, France, ³ Research Group PLECO (Plants and Ecosystems), Department of Biology, University of Antwerp, Antwerp, Belgium, ⁴ Department of Yield and Silviculture, Slovenian Forestry Institute, Ljubljana, Slovenia, ⁵ Royal Museum for Central Africa, Service of Wood Biology, Tervuren, Belgium

## OPEN ACCESS

**Edited by:**

Fabio Gennaretti,
Université du Québec en Abitibi-
Témiscamingue, Canada

**Reviewed by:**

Emanuele Ziaco,
Uppsala University, Sweden
Valentina Buttò,
Université du Québec à
Chicoutimi, Canada

**\*Correspondence:**

Romain Lehnebach
lehnebach.romain@hotmail.fr

**Specialty section:**

This article was submitted to Technical Advances in Plant Science, a section of the journal Frontiers in Plant Science

**Received:** 21 April 2021

**Accepted:** 28 June 2021

**Published:** 06 August 2021

**Citation:**

Lehnebach R, Campioli M, Gričar J, Prislan P, Mariën B, Beeckman H and Van den Bulcke J (2021) *High-Resolution X-Ray Computed Tomography: A New Workflow for the Analysis of Xylogenesis and Intra-Seasonal Wood Biomass Production*. *Front. Plant Sci.* 12:698640. doi: 10.3389/fpls.2021.698640

Understanding tree growth and carbon sequestration are of crucial interest to forecast the feedback of forests to climate change. To have a global understanding of the wood formation, it is necessary to develop new methodologies for xylogenesis measurements, valid across diverse wood structures and applicable to both angiosperms and gymnosperms. In this study, the authors present a new workflow to study xylogenesis using high-resolution X-ray computed tomography (HRXCT), which is generic and offers high potential for automatization. The HXRCT-based approach was benchmarked with the current classical approach (microtomy) on three tree species with contrasted wood anatomy (*Pinus nigra*, *Fagus sylvatica*, and *Quercus robur*). HRXCT proved to estimate the relevant xylogenesis parameters (timing, duration, and growth rates) across species with high accuracy. HRXCT showed to be an efficient avenue to investigate tree xylogenesis for a wide range of wood anatomies, structures, and species. HRXCT also showed its potential to provide quantification of intra-annual dynamics of biomass production through high-resolution 3D mapping of wood biomass within the forming growth ring.

**Keywords:** high-resolution X-ray computed tomography, microtomy, tree growth, xylogenesis, secondary growth phenology

## INTRODUCTION

Forests represent the main CO₂ sink of terrestrial ecosystems. Understanding tree growth and carbon sequestration is of crucial interest to forecast the responses of forests to climate change and their climate mitigation potential. The forest carbon sink can be assessed with direct measurement of tree growth (i.e., Hubau et al., 2020) or indirectly through Dynamic Global Vegetation Models (DGVM). In DGVM, tree growth is often assumed to be limited primarily by photosynthesis (i.e., gross primary production) (Fatichi et al., 2019).
However, experimental evidence shows that tree growth does not directly depend on the photosynthesis and that tree growth is more sensitive to environmental conditions than the photosynthesis (Fatichi et al., 2014). Therefore, more attention should be devoted to determine the direct constraints on growth exerted by the environmental conditions at different temporal scales. One way to do this is by studying wood formation (i.e., xylogenesis) that defines seasonal tree growth and carbon sequestration (Friend et al., 2019).

Xylogenesis is driven by the cambium, that is, a layer of meristematic cells, forming a cylinder along the different axes of the plant body. For temperate species, when growing conditions become favorable in spring, the cambium cells start to divide and the newly formed cells undergo different structural changes (i.e., differentiation). Xylogenesis can be divided into five sequential processes: (1) cell division and (2) cell enlargement, which drive the growth in size; (3) cell-wall thickening and (4) cell-wall lignification, which are mainly responsible for the increase in mass and C sequestration (Cuny et al., 2015); and (5) programmed cell death, which allows the woody cell skeleton to be used as mechanical support and conduit for water transport. The wood formation and its dynamics are strongly influenced by environmental conditions, such as the temperature, day length, and water availability (Delpierre et al., 2016; Dox et al., 2021). Therefore, climate change is expected to affect the tree growth and carbon accumulation through the modulation of the timing and duration of xylogenesis.

An appealing method to study the phenology and the seasonal dynamics of wood formation as well as the climate-wood growth relationships consists in monitoring wood development in small cores (i.e., microcores) collected at regular time intervals during the growing season (Rossi et al., 2006). Seasonal observations of anatomical features characterizing the wood formation (e.g., cell morphology, cell developmental stage, or number of cells in a given stage) are made possible by brightfield microscopy after histological slicing and staining. Most of the xylogenesis studies and the development of the histological approach (referred hereafter as “microtomy”) have been performed on conifers with a homogeneous wood structure (Rossi et al., 2009, Deslauriers et al., 2008, Camarero et al., 2010). This intensive work on conifers also permitted to refine the understanding of the asynchrony between xylem size increase and biomass production through a detailed quantification of cell development kinetic during xylogenesis (Cuny et al., 2015). This approach allowed to develop a clear view on xylogenesis in conifer species (Rathgeber, 2017), while its knowledge on angiosperm species remains limited (ref. Delpierre et al., 2016, but see Cufar et al., 2008; Gričar et al., 2018; Prislan et al., 2018; Dox et al., 2020). In fact, the structural heterogeneity of angiosperm wood challenges the acquisition of xylogenesis data from a methodological point of view. Indeed, the diversity of cell types encountered in angiosperm wood (vessels, fibers, and parenchyma) makes cell-counting and morphometrics difficult, while limiting the acquisition of information on the xylogenesis dynamics. Therefore, even if some authors have been able to quantify the proportion of growth increment attributed to the different stages of

development (Prislan et al., 2018), intra-seasonal studies of secondary growth in angiosperms have been often limited to the estimation of the timing of qualitative events (e.g., earlywood vessel formation, latewood formation). From the biomass accounting point of view, bypassing cell-level information by characterizing the biomass accumulated during the season at the tissue level is a promising perspective. In this respect, the methodology of Andrianantenaina et al. (2019), based on an estimation of the apparent density of the forming wood, should allow to overcome the issue of structurally heterogeneous wood and thus to quantify the dynamics of intra-annual biomass production in angiosperms. The global predominance of angiosperms stresses the importance of further methodological development allowing straightforward xylogenesis measurements across the diversity of wood structures and species found within the angiosperm group, and across the different climates and biomes.

Investigating how tree growth interacts with environmental factors at a global scale requires the development of a high-throughput methodology for xylogenesis studies applicable also to angiosperms. Such methodology should avoid the histological preparation and analysis that are labor intensive, requiring specific expertise, and several manual steps, such as paraffin embedding, slicing, staining, mounting, and observation using bright field microscopy. So far, the labor intensity of the manual steps during sample preparation (requiring skilled technicians) and data collection (observations, requiring skilled wood anatomists) has limited the automation of the microtomic approach, pointing out the importance to develop an automated methodology to study seasonal wood formation.

High-resolution X-ray computed tomography (HXRCT) is becoming popular in the analysis and imaging of plant anatomy (Brodersen and Roddy, 2016). Different research fields related to tree growth and functioning, such as dendrochronology (De Mil et al., 2016; Van den Bulcke et al., 2019) and ecophysiology (Knipfer et al., 2015; De Baerdemaeker et al., 2019), seized the opportunity of HXRCT for substantial methodological improvements. However, HXRCT has never been used in the demanding field of xylogenesis. Producing three-dimensional gray scale images with pixel values (attenuation coefficients) directly related to the density of the material, HXRCT represents a powerful mean to study xylogenesis. More precisely, HXRCT can be used to track cell wall density change occurring during the formation and thickening of the secondary cell wall in order to differentiate the mature wood from the forming wood and to assess secondary growth phenology. In addition, HRXCT can quantify the biomass accumulated by trees throughout the course of the season.

The use of HXRCT in the framework of xylogenesis requires both high contrast and resolution to discern the forming cells and the cell wall. Obtaining a sufficient contrast between cell lumens and cell walls with HRXCT requires water removal from the sample. For mature xylem, water is generally removed by oven-drying the sample prior to scanning (De Mil et al., 2016). However, drying without collapsing and destroying the soft tissue, such as cambial and newly forming xylem cells, cannot be
achieved through oven- or air-drying (Stuppy et al., 2003; Leroux
et al., 2009).

Therefore, monitoring of xylogenesis with HXRCT requires
the development of an entire workflow including sample
preparation, image acquisition, and an analysis protocol,
alongside a methodology to follow the cell wall density change
within the forming growth ring.

Here, a new methodology is presented and validated to
study xylogenesis using HXRCT based on a new workflow to
monitor cell wall density changes occurring during xylogenesis
(Figure 1), applied to three tree species (black pine, common
beech, and pedunculate oak) with contrasted wood anatomy
growing in a temperate region. Several relevant parameters,
such as the timing and duration of size growth, mature xylem
production, and overall xylogenesis, as well as rates of size
growth and mature xylem production, were estimated with
HXRCT and validated with the classical microtomy approach.
In addition, the gray scale images produced by HRXCT were
used to account for seasonal biomass production dynamics.
Concerning the latter aspect, the authors ask whether the intra-
seasonal stem biomass production and size growth differ between
conifers and angiosperms, specifically among pine, beech, and
oak. While information on intra-seasonal biomass production
exist for conifers (Cuny et al., 2015; Andrianantenaina et al.,
2019) to the best of our knowledge, no study has investigated the
intra-seasonal biomass production in angiosperms.

MATERIALS AND METHODS

Selected Species, Study Sites, and
Sampling

To assess the applicability of the HXRCT method on tree
species with different functional strategies, three species growing
in Belgium were selected: black pine (Pinus nigra), common
beech (Fagus sylvatica), and pedunculate oak (Quercus robur).
With 778 mm.year⁻¹ of precipitation and an average monthly
temperature ranging from +3°C to 17°C (average annual
temperature of 10.1°C), Belgium is characterized by a maritime
temperate climate (Campioli et al., 2012). Beech and oak trees
were sampled in the “park of Brasschaat” (51°12′ N–4°26′ E)
while pine trees were sampled in the military domain “Klein
Schietveld” (51°21′ N–4°37′ E) (Dox et al., 2020; Marchand
et al., 2020). The sites are at a distance of 9.5 km in a lowland
(18–22 m a.s.l.) region with a similar climate. The three species
present contrasted xylem anatomy: beech and oak are diffuse
and ring porous species respectively, whereas black pine has the
typical conifer wood structure made mostly of tracheids. Four
healthy and (co)dominant trees per species were selected. From
the beginning of April to mid-November 2019, two microcores
were collected on each tree every two weeks using a Trophor
tool resulting in 18 sampling dates. The sampling was performed
between 1.3 m and 2 m above the ground following an upward
spiraling pattern with a distance of 4 cm between consecutive
sampling points to avoid wound reactions. The microcores were
immediately transferred to an Eppendorf tube filled with 70%
ethanol solution.

X-Ray CT Workflow

Dehydration, Critical Point Drying, and Sample
Mounting

After dehydration by increasing ethanol concentrations, the
samples for HXRCT were dried to obtain a sufficient contrast
between cell lumens and cell walls. Therefore, a gentle drying
method was required in order to preserve the structure of the soft,
newly forming wood cells. Therefore, critical point drying (CPD),
a technique commonly used for the preparation of scanning
electron microscopy samples (Bray, 2000), was used, but has
proven its efficiency in the drying of fresh plant structures also
for CT imaging (See **Supplementary Figure 1** for a comparison
of an air dried and CPD processed pine microcore) (Leroux et al.,
2009; Gutiérrez et al., 2018).

The samples were critically point dried with liquid CO₂ using
a CPD device (Baltec CPD 030). Each drying cycle consisted
of a mixing-draining procedure (repeated 10 times) allowing
the gradual replacement of the dehydration media (i.e., 100%
ethanol) by liquid CO₂. Then, the samples were dried by
removing the CO₂ in a supercritical state (> 31°C, > 73 bar).

Six samples were dried at the same time. After drying, the
samples were mounted on a cylindrical carbon stick (5 cm
x2 mm) for scanning. The mounted samples were stabilized in the
scanner room for 12 h, at 25°C and a relative humidity of 40%, to
avoid sample movement during scanning.

X-Ray CT Scanning and Reconstruction

Samples were scanned with the Nanowood X-ray CT scanner,
developed at the Ghent University Centre for X-ray tomography
(UGCT; http://www.ugct.ugent.be) (Dierick et al., 2014). The
directional X-ray tube (Microfocus, Hamamatsu) was operated at
a tube voltage of 70 kV and a target current of 100 µA. A large (25
x 20 cm, 1920 x 1536 px) Varian Si flat detector was used allowing
to capture the region of interest (i.e., innermost phloem layers,
cambium, and forming xylem growth ring) in a single field of
view. A total of 2001 projections were taken with an exposure
time of 700 ms per projection, resulting in 25 min per scan.

Scanned volumes were reconstructed with the Octopus
reconstruction software (Vlassenbroeck et al., 2007), including
Paganin phase contrast filtering (Paganin et al., 2002). The
reconstructed volume, with an approximate voxel pitch of
2.49 µm, consisted of 16-bit tangential images running from the
bark to the inner side of the xylem.

Processing and Analysis of the
Reconstructed Volume

Selection of the Regions of Interest and Gray-Value
Extraction

When sampling, the axis of the Trophor tool can substantially
deviate from both the radial direction and the transverse plane of
the stem. Therefore, the reconstructed volume was reoriented to
realign the cambium surface parallel to the volume edge in both
transverse and longitudinal plane (Supplementary Figure 2)
using ImageJ (Rasband, 2012). A volume of interest (VOI)
containing the inner layers of the phloem, the cambium, and the
<figure><img src="image_4.png" /><figcaption>**FIGURE 1** | Comparison between the microtomy and HXRCT workflow developed in the present study. The flowchart covers the different steps from field sampling to data analysis through sample preparation and data extraction.</figcaption></figure>

forming xylem growth ring was selected while avoiding selection of wide rays, resin ducts, or cracks.

This VOI was then processed with Octopus analysis, formerly distributed by the UGCT spin-off company XRE (now TESCAN-XRE, part of the TESCAN ORSAY HOLDING a.s.) (Brabant et al., 2011). For each tangential slice of the VOI, a first thresholding was applied to create a mask covering the area containing the sample. The masks were used to create a first region of interest (ROI) image stack containing xylem pixels (xylem ROI stack). For each tangential slice of the xylem ROI stack, a second thresholding was applied to create a mask excluding the pixels associated to cell lumens. The masks were used to create a second region of interest image stack containing cell wall pixels exclusively (cell wall ROI stack). Both xylem and cell wall ROI stacks of each microcore were then processed with R software (R Development Core Team, 2019), using the packages "magick" (Ooms, 2020) and "imager" (Barthelme, 2019), to extract the mean gray value per ROI resulting in two radial gray-value profiles per microcore, herein after referred to as "xylem gray profile" and "cell wall gray profile," respectively.

## Cell Wall Gray-Value Profiling: Rationale

The gray-value profiling methodology developed here relies on the positive and linear relationship between the gray value of HXRCT-derived images and the density of the scanned material. Therefore, cell wall density changes occurring within the forming growth ring can be assessed through cell wall gray value changes. From a structural point of view, the forming growth ring consists of a maturing xylem zone (where xylogenesis occurs), and a mature xylem zone, where xylogenesis already ended (Figure 2).

Within the maturing xylem, the secondary cell wall is formed and lignin is deposited resulting in an increase of the cell wall density from the cambium until the boundary with the mature xylem zone.

Apart from secondary changes occurring many years later, such as heartwood formation, the structure of the fully mature cell wall does not change in time (Lehnebach et al., 2019), its density is assumed to vary both between and within species in a very narrow range (1.4–1.5 g/cm3). Cell wall density, and consequently cell wall gray value are therefore assumed to remain constant within the mature xylem zone. Consequently, the cell wall density profile along the growing ring follows a two-phased pattern with (1) an increase from cambium to the end of the maturing xylem zone followed by (2) a constant density toward the innermost layers of the mature xylem zone.

Identifying the distance from the cambium at which the phase transition occurs, allows estimating the width of both maturing and mature xylem. An analytical approach is developed in this study, aiming to quantify the width of both maturing and mature xylem.

## Cell Wall Gray-Value Profiling: Implementation

For each sample, the cell wall gray profile was plotted over a transversal reference image, with each data point corresponding to the average gray value of the corresponding tangential slice (Figure 3). We used the smoothing property of a generalized additive model (GAM) to analyze the cell wall gray profile of each microcore. GAMs are extremely versatile in fitting non-linear and non-monotonic data series, and offer the possibility to
<figure><img src="image_5.png" /><figcaption>**FIGURE 2** | Variation of cell wall gray value or density within the forming growth ring. Following cell wall maturation during xylogenesis, the cell wall density (red line) is expected to increase from the cambium (green dotted line) to the most recent mature xylem layer (blue dotted line). Within the mature xylem, the cell wall is fully developed and its basic density is assumed to be constant until the growth ring boundary (yellow dotted line) and beyond. The distance from the cambium at which the phase transition occurs, allows estimating the width of both mature and maturing xylem.</figcaption></figure>

estimate and analyze their derivative. GAMs are therefore well-suited to analyze cell wall gray profile of diverse tree species. All GAMs were fitted with a Gaussian distribution and identity link function. The goodness of each fit was assessed visually.

Consisting of primary-walled cells, the cambium and the enlarging phloem and xylem cells have the lowest cell wall density (Figures 2, 3, 4A,B). The cell wall density increases at both sides of the cambium because of phloem and xylem cell wall differentiation (Figures 4A,B). We took advantage of this feature to define the cambium position semi-automatically by fitting a GAM and analyzing its derivative within a manually defined subset of the sample containing the cambium and the newly differentiating phloem and xylem cells (Figure 4B). Within this subset, the pointwise 95% confidence intervals of the GAM derivatives were used to detect the part where the GAM predicted values were constant (i.e., where the fitted GAM derivative values were significantly different from 0, i.e., 0 lying outside the

confidence interval). The median radial position of this part was fixed as the cambium position (Figure 4B). The radial position of the previous growth ring boundary was indicated manually on the reference image (Figure 4A). It allowed to estimate the width of the increment, $w_{incr-CT}$.

The distance from the cambium for each data point was calculated and fitted another GAM over the data spanning from the cambium (i.e., distance from the cambium = 0) to the previous growth-ring boundary. As for cambium detection, the 95% confidence intervals of the GAM derivatives allowed the identification of growth ring parts where the gray values are constant. The constant part closest to the cambium was selected as the region where the transition between maturing and mature zone occurs (Figure 4C). The radial position of the left edge of the transition region was considered as the width of the maturing xylem, $w_{maturing-CT}$. Finally, the width of the mature zone was computed as, $w_{mat-CT} = w_{incr-CT} - w_{maturing-CT}$.
<figure><img src="image_6.png" /><figcaption>**FIGURE 3** | Plot of the cell wall gray profile. The volume of interest (VOI) containing the inner layers of the phloem, the cambium, and the forming xylem growth ring (top) is sliced tangentially (middle). Averaging the cell wall gray values of each tangential slices results in a cell wall gray value profile plotted over the transversal reference image (bottom). The reference image is in the red frame (top and bottom image), the cambium and previous ring boundary are highlighted by the green and yellow lines respectively.</figcaption></figure>

GAM fitting and derivative analysis were performed using the R packages "mgcv" (Wood, 2017) and "gratia" (Simpson, 2020). A step by step description of the gray-value profiling methodology is available in **Supplementary File 1**.

## Microtomy Workflow

The microcores for microtomy were dehydrated in graded ethanol series, infiltrated, and embedded in paraffin blocks and finally sliced using a rotary microtome. The thin sections were stained in a water solution of safranin and astra blue, staining lignin in red and cellulose/hemicellulose in blue, respectively. Histological sections were prepared at the Slovenian Forestry Institute (Ljubljana), following the detailed protocol of Prislan et al. (2014). Observation under light microscope allowed to distinguish between maturing and mature xylem. For each microcore, the width of the cambium, and the width of the maturing and mature xylem were measured at three positions to consider potential within-sample heterogeneity. The maturing zone was measured as the ring part between the cambium and the first mature cells. Xylem cells were considered mature when the entire cell wall was stained red and the lumen was empty of any cellular content. In order to be comparable to the HXRCT methodology, the increment width was computed as the sum of both mature and maturing xylem width and the width of the cambium divided by two. When the maturing zone was absent (i.e., for samples collected at the very end of the growing season), the mature zone was set equal to the increment width.

## Estimation of Xylogenesis Parameters and Statistical Analysis

### Computation of Xylogenesis Critical Dates and Duration, and Growth Rate

Computation of xylogenesis phenological parameters was first performed through a "model-driven" approach, such as the fitting and derivation of the Gompertz function (e.g., Rossi et al., 2003). However, a recent methodological breakthrough suggests that a "data-driven" approach is able to capture properly the variation of xylogenesis dynamics (i.e., durations and rates) while the model-driven approach fails (Cuny et al., 2013). Preliminary analysis of the data (not shown) using a Gompertz function produced improper fitting and led us to opt for the GAM data-driven approach.

As the increase of the increment width during the growing season is monotone, a shape constrained additive model (SCAM) was used (Pya and Wood, 2015) to analyze the phenology of wood formation. SCAM was fitted using monotone increasing P-splines, Gaussian distribution, and identity link function. Circumferential growth heterogeneity was accounted for using the correction formula for cell counting data adapted to widths at the individual tree scale:

$$w_{cor,j} = w_{obs,j} * \frac{w_{prev,j}}{\frac{1}{n} \sum w_{prev,j}} \quad (1)$$
<figure><img src="image_7.png" /><figcaption>**FIGURE 4** | Implementation of the gray-value profiling. (A) The average gray value of each tangential slice is plotted on the reference cross-section image. At this stage, the previous growth ring boundary is identified (yellow vertical line). (B) Within the subset containing the cambium, a GAM fit (red) and its derivative allows automatic identification of the cambium position (green vertical line). (C) Within the growth ring, a GAM fit (red) and its derivative allows the detection of a constant gray value (blue filled circles). The constant gray value closest to the cambium is set as the boundary between maturing and mature xylem (blue vertical line). This procedure allows to estimate maturing xylem width (*w*<sub>maturing</sub>), mature xylem width (*w*<sub>mature</sub>), and increment width (*w*<sub>incr</sub>). The gray scale is unitless. It varies from 0 to 65,535 for 16 bits images used here.</figcaption></figure>

where, *n* is the number of sampling dates, and for each *j*th sample, *w*<sub>cor,j</sub> is the corrected width; *w*<sub>obs,j</sub> is the observed width; *w*<sub>prev,j</sub> is the observed width of the previous growth year (Rossi et al., 2003; Prislan et al., 2013).

First, the constrained GAM was fitted over the width of the increment to analyze the size growth dynamics at the tree level (Figure 5). The final increment width was computed as the median of the increment width predicted for the last three sampling dates (i.e., the dates for which the increment width had already stabilized), in order to account for circumferential growth heterogeneity. We defined the date of beginning (*tb*<sub>incr</sub>) and ending (*te*<sub>incr</sub>) of size growth as the day of year for which, respectively, 5% and 95% (i.e., *wb*<sub>incr</sub> and *we*<sub>incr</sub>) of the increment is formed. The duration of size growth was computed as the time span between both dates:

$$d_{incr} = te_{incr} - tb_{incr} \quad (2)$$

The average daily rate of size growth was calculated as:

$$r_{incr} = (we_{incr} - wb_{incr})/(te_{incr} - tb_{incr}) \quad (3)$$

The differences between the increment width predicted by GAM on the consecutive days, spanning from *tb*<sub>incr</sub> to *te*<sub>incr</sub>, were calculated to assess the distribution of size growth daily rates.

The same procedure was applied for the width of the mature xylem in order to define the date of beginning (*tb*<sub>mat</sub>) and ending (*te*<sub>mat</sub>), the duration (*d*<sub>mat</sub>) and the average rate (*r*<sub>mat</sub>) of mature xylem production. Finally, the duration of the whole xylogenesis was defined as:

$$d_{xylo} = te_{incr} - tb_{incr} \quad (4)$$

These parameters were estimated for both microtomy and HXRCT datasets.

## Consistency Between HXRCT and Microtomy Measurements

The consistency between increment, mature, and maturing xylem width measured with HXRCT and microtomy was assessed using major axis regression. The value of both intercept and slope and their respective 95% confidence intervals (CI) were used to detect significant differences compared to a 1:1 regression. An intercept differing from 0 (i.e., 0 lies outside the CI of the intercept) indicates a fixed bias across the entire range, that is, the HXRCT provides greater (or lower) values than microtomy. A slope differing from 1 (i.e., 1 lies outside the CI of the slope) indicates a proportional bias, that is, the HXRCT provided progressively greater (or lower) values with increasing values of the considered measurement than microtomy. Major axes were fitted with the "lmodel2" package in R (Legendre, 2018).
<figure><img src="image_8.png" /><figcaption>**FIGURE 5** | Description of the xylogenesis parameters estimated for both HXRCT and microtomy datasets. Parameters related to size growth (green) and mature xylem production (brown) are reported. The date of beginning (*tb*<sub>incr</sub>) and ending (*te*<sub>incr</sub>) of size growth are defined as the day of year for which, respectively, 5% and 95% (i.e., *wb*<sub>incr</sub> and *wc*<sub>incr</sub>) of the increment is formed. The date of beginning (*tb*<sub>mat</sub>) and ending (*te*<sub>mat</sub>) of mature xylem production are defined as the day of year for which, respectively, 5% and 95% (i.e., *wb*<sub>mat</sub> and *wc*<sub>mat</sub>) of the mature xylem is formed. The duration of size growth (*d*<sub>incr</sub>), mature xylem production (*d*<sub>mat</sub>) and xylogenesis (*d*<sub>xylo</sub>) are computed as the time span between *te*<sub>incr</sub> and *tb*<sub>incr</sub>, *te*<sub>mat</sub> and *tb*<sub>mat</sub>, and *te*<sub>mat</sub> and *tb*<sub>incr</sub> respectively. The average daily rate of size growth and mature xylem production are calculated as, $r_{incr} = (w_{incr} - w_{incr}) / (te_{incr} - tb_{incr})$ and $r_{mat} = (w_{mat} - w_{mat}) / (te_{mat} - tb_{mat})$.</figcaption></figure>

## Comparison of Parameters Estimated With HXRCT and Microtomy

To assess the agreement between critical dates, durations, and growth rates estimated using HXRCT and microtomy, the significant differences between the two datasets by case-based bootstrap resampling were tested. For a given species and for each parameter of interest (Figure 5), the difference between the median of the four trees based on HXRCT data on one hand, and microtomy data on the other hand was computed. The observed median difference was compared to the median reference distribution simulated under the null hypothesis (H0: the parameters extracted from HXRCT do not differ from the parameters extracted from microtomy). The reference distribution resulted in 10,000 resampled median differences. A two-tailed test was performed by computing the achieved significance level (ASL) as the proportion of absolute values in the reference distribution that are equal, or greater, than the observed absolute median difference. ASL is interpreted in the same way as the p value (P).

## Accounting for Biomass Production Dynamic Through Xylem Gray Value

The ability of HRXCT to assess biomass production dynamics throughout the course of the growing season was assessed by using the gray-value information. Here, the xylem gray profile (i.e., including lumen and cell wall pixels) was preferred to the cell wall gray profile as it accounts for both the proportions of the xylem layer allocated to cell wall (biomass) and lumens (void). The sum of the gray values of each xylem gray profile

was calculated. Being linearly related to the biomass accumulated in a given volume of the forming xylem ring, this variable was considered as a proxy of the biomass accumulated at the trunk base (*b*<sub>trunk</sub>) of the tree from the beginning of the growing season up to the considered sampling date. As for xylem widths, we accounted for circumferential growth heterogeneity by correcting *b*<sub>trunk</sub> by the observed width of the previous growth year at the individual tree scale. We fitted constrained GAM over *b*<sub>trunk</sub> to analyze the biomass production dynamics at the species level. The biomass production and the size increase were compared by analyzing the seasonal variation of both *b*<sub>trunk</sub> and the increment width (*w*<sub>incr</sub>). For convenience, both variables were normalized (i.e., ranging between 0 and 1) using the median of *b*<sub>trunk</sub> and *w*<sub>incr</sub>, respectively, predicted for the last three sampling dates. This provided the relative completion of each process (i.e., size growth and biomass production) and allowed to compute the time lag between biomass production and size growth on the course of the growing season. For this purpose, we computed the time lag between the relative completion of a given proportion of the final *b*<sub>trunk</sub> and the final *w*<sub>incr</sub> (e.g., Day of year (DOY) when a given relative completion of the final *b*<sub>trunk</sub> is produced - DOY when a given relative completion of the final *w*<sub>incr</sub> is produced). The time lags were computed for relative completion values ranging between 0 and 1 and every 0.01. The differences between both *b*<sub>trunk</sub> and *w*<sub>incr</sub> predicted by GAM on the consecutive days were also calculated. This allowed to assess the seasonal variation of size growth rate and biomass production rate and estimate their seasonal maximum.

# RESULTS

## Seasonal Growth Pattern

The patterns of seasonal xylem growth, at the species level, derived from HXRCT and microtomy strongly overlap (Figure 6). For both pine and beech, increment and mature xylem width increased along the growing season following a typical sigmoidal curve (Figure 6). For oak, the width of the increment follows a curvilinear increase while the mature xylem width follows a bisigmoidal pattern. Major axis fits did not detect any significant bias in the estimation of increment width by HXRCT compared to microtomy (Supplementary Figure 3; Supplementary Table 1). Regarding the mature xylem width, no biases were detected for beech and oak, however, a slight slope departure from 1 for pine (Supplementary Table 1) is found. Microtomy and HRXCT described similarly the seasonal variation of maturing xylem width. For beech, the maturing xylem width followed a bell-shaped pattern while for oak and pine the pattern was bimodal, with a first peak in spring and a second peak in July for oak and in September for pine.

Overall, the estimated maturing xylem width by HXRCT is significantly lower than for microtomy for pine and beech (see Supplementary Table 1), especially during summer (Figure 6). These represent the only cases in our study when HXRCT and microtomy produce different estimates (see the third paragraph of the discussion section).
<figure><img src="image_9.png" /><figcaption>**FIGURE 6** | Seasonal growth pattern assessed with HXRCT (red) and microtomy (blue). Variations of increment and mature xylem were fitted with constrained GAMs (top and middle row), while maturing xylem width was fitted with unconstrained GAMs (bottom row). The shaded band represents the standard error of the fit.</figcaption></figure>

## Wood Formation Calendar

Size growth, and therefore xylogenesis, began ($tb_{incr}$) at the end of March in oak (one week before the first sampling date, $tb_{incr-micro}$ = DOY 89, while $tb_{incr-CT}$ = DOY 88, see **Supplementary Table 2**), while it started between mid- and end of April for pine ($tb_{incr-micro}$ = DOY 100, while $tb_{incr-CT}$ = DOY 106) and beech ($tb_{incr-micro}$ = DOY 122,

while $tb_{incr-CT}$ = DOY 108, the difference is not significant). Microtomy and HXRCT gave similar estimates of $tb_{incr}$ (**Figure 7; Supplementary Table 2**). Size growth was completed in mid-September for beech ($te_{incr-micro}$ = DOY 256, while $te_{incr-CT}$ = DOY 251), and oak ($te_{incr-micro}$ = DOY 260, while $te_{incr-CT}$ = DOY 264) while it lasted until the beginning of October for pine ($te_{incr-micro}$ = DOY 277, while $te_{incr-CT}$ = DOY 270). Again,
microtomy and HXRCT provided similar estimate of $te_{incr}$ and therefore the duration of size growth ($d_{incr}$) was highly similar for both methods ($d_{incr-micro}$ = DOY 177, $d_{incr-CT}$ = DOY 162 for pine, $d_{incr-micro}$ = DOY 134, $d_{incr-CT}$ = DOY 140 for beech and $d_{incr-micro}$ = DOY 163, $d_{incr-CT}$ = DOY 171 for oak, the differences are not significant, see **Supplementary Table 2**).

The mature xylem production began in April for oak ($tb_{mat-micro}$ = DOY 105, while $tb_{mat-CT}$ = DOY 104), and from the end of May to beginning of June for beech ($tb_{mat-micro}$ = DOY 145, while $tb_{mat-CT}$ = DOY 141) and pine ($tb_{mat-micro}$ = DOY 152, while $tb_{mat-CT}$ = DOY 154). However, the end of mature xylem production, and therefore the end of xylogenesis, occurred by the end of September for beech ($te_{mat-micro}$ = DOY 265, while $te_{mat-CT}$ = DOY 264) and oak ($te_{mat-micro}$ = DOY 262, while $te_{mat-CT}$ = DOY 262), and by the end of October for pine ($te_{mat-micro}$ = DOY 304, while $te_{mat-CT}$ = DOY 294, the difference is not significant). Overall the duration of mature xylem production ($d_{mat}$) is shorter for beech ($d_{xylo-micro}$ = 142 days, while $d_{xylo-CT}$ = 154 days, the difference is not significant) than for pine ($d_{xylo-micro}$ = 206 days, while $d_{xylo-CT}$ = 188 days, the difference is not significant) and oak ($d_{xylo-micro}$ = 166 days, while $d_{xylo-CT}$ = 169 days, the difference is not significant). Dates ($tb_{mat}$, $te_{mat}$) and durations ($d_{mat}$) related to mature xylem formation as well as the duration of xylogenesis ($d_{xylo}$) were similar across methods (**Figure 7; Supplementary Table 2**).

## Growth Rates

The estimated size growth rate ($r_{incr}$) ranged from 5 to $6 \mu m.day^{-1}$ in pine to 21 to $22 \mu m.day^{-1}$ in beech. The mature xylem production rate ($r_{mat}$) presented similar values. Overall, estimates of size growth and mature xylem production rate by microtomy and HXRCT are highly similar (**Figure 8; Supplementary Table 3**). Moreover, the distributions of $r_{incr}$ and $r_{mat}$ estimated with HXRCT and microtomy strongly overlap (**Supplementary Figure 3**).

## Seasonal Pattern of Biomass Production and Size Growth and Their Time Lag

The patterns of biomass production and size growth are very similar at the species level (**Figure 9**). For both pine and beech, the biomass increases during the growing season following a sigmoidal curve, while biomass in oak increases curvilinearly. However, biomass production lags behind size growth. The time lag observed in respect of both size and biomass completion of the xylem ring is more important for pine (mean = 9 days, maximum=15 days) than for beech (mean = 4 days, maximum = 9 days) and oak (mean = 5 days, maximum = 7 days). For pine and beech, the time lag reaches its maximum during the latewood formation (at the end of the growing season, when 80–90% of the ring width is formed). In oak, the time lag is maximum from the end of the earlywood formation up to half of the latewood formation (between 40% and 70% of the ring formation).

The time lags observed at the peak of both size growth and biomass production rate contrast with the time lag observed in respect of both size growth and biomass completion. For both pine and beech, the relative size growth rate, as well as the biomass production rate peaked at the end of June. The

difference between maximum size growth (DOY = 173) and biomass production rates (DOY = 177) was four days in pine. This difference was even smaller in beech, that is, two days (DOY = 171 vs. DOY = 173). For oak, the relative size growth rate was maximal at the end of March (DOY = 89), while the biomass production rate peaked at the beginning of June (DOY = 152), resulting in a time lag of 63 days.

## DISCUSSION

This is the first time that xylogenesis dynamics and intra-seasonal growth rate parameters are assessed using HRXCT. The development of an entire HRXCT workflow from sample preparation, to reconstructed volume analysis, allowed the estimation of relevant xylogenesis parameters (i.e., dates, durations, and rates) and seasonal time series of increment width. Moreover, the benchmarking of HRXCT with the classical microtomy proved the HRXCT approach to produce accurate estimates of xylogenesis dynamics over species with contrasted xylem anatomies.

The ability of HRXCT to encompass diverse xylem structure relies on the gray-value profiling methodology. While 2D cell morphometric analysis is not applicable on various anatomical structure, the gray-value profiling overcomes this issue by assessing the maturation state from monitoring of density at the xylem layer scale. The HRXCT methodology also proved to produce accurate estimates of actual ring width along the whole growing season ranging from ca. 30 to 5000 µm (see **Supplementary Table 4**). Given this result and the high resolution of HRXCT (i.e., 2.5 µm/pixel), the authors assume that HRXCT will succeed in estimating xylogenesis dynamics of slow-growing trees. Therefore, the generic and data-driven HRXCT method showed to offer a very interesting perspective to fill the knowledge gap about secondary growth phenology and intra-annual wood production over a structural and taxonomical diversity.

It is important to pinpoint that the xylem dimensions estimated with microtomy and HRXCT slightly differ in their nature. HRXCT differentiates the mature xylem from the maturing xylem on the basis of cell wall density differences. Therefore, cells in the very late stages of wood formation, for which the biomass accumulation process is almost accomplished (i.e., cells at the end of secondary wall formation and lignifying cells) or terminated (i.e., fully lignified cells but still alive) might be considered as mature cells with HRXCT. On the other hand, with microtomy such cells are considered as still maturing since the mature zone consists exclusively of fully lignified dead cells. Therefore, the slight overestimation and underestimation of the mature and the maturing zone, respectively, with HRXCT may rely on the different nature of xylem dimensions estimated with HRXCT and microtomy. However, these differences are not sufficiently important to bias the estimation of xylogenesis dynamics.

As gray scale information can be directly converted into estimates of biomass, X-ray computed tomography has been successfully used in the field of dendrochronology and biomass
<figure><img src="image_10.png" /><figcaption>**FIGURE 7** | Xylogenesis calendar and duration. Left panels: Estimates of size growth beginning (*tb*<sub>incr</sub>) (i.e., beginning of xylogenesis) and size growth ending (*te*<sub>incr</sub>) (green) as well as mature xylem production beginning (*tb*<sub>mat</sub>) and ending (*te*<sub>mat</sub>) (i.e., ending of xylogenesis) (brown) are presented for HXRCT and microtomy. Green and brown lines range from minimum to maximum estimated dates, while second and third quartile are represented by short vertical lines. Right panels: Estimates of xylogenesis duration (*d*<sub>xylo</sub>) (median) for HXRCT (red) and microtomy (blue) are presented. Black lines range from minimum to maximum estimated duration, while second and third quartile are represented by short vertical lines. The significance of the bootstrap resampling test (P) assessing the difference between HXRCT and microtomy is reported.</figcaption></figure>

estimation (Bastin et al., 2015; De Mil et al., 2016). A significant step forward in the quantification of the intra-annual dynamics of carbon sequestration, for both conifers and angiosperms, relies on the calculation of the apparent density of the forming wood (i.e., the product of cell wall density and the ratio of cell wall area to xylem area) (Andrianantenaina et al., 2019). Although the apparent density approach proved to be more direct and

accurate than the cell morphometric approach (Cuny et al., 2015), it still relies on physical 2D slicing of the sample. Here, it is showed that HRXCT can also be used in the field of intra-seasonal biomass production and carbon sequestration research. In fact, the sum of the xylem gray value, which is actually directly related to density (attenuation coefficients are recorded), and therefore to the biomass accumulated at the trunk base provided a very
<figure><img src="image_11.png" /><figcaption>**FIGURE 8** | Size growth ($r_{incr}$) and mature xylem ($r_{mat}$) production rates (median) estimated with CT (red) and microtomy (blue). Black lines range from minimum to maximum estimated duration, while second and third quartile are represented by short vertical lines. The significance of the bootstrap resampling test (P) assessing the difference between HXRCT and microtomy is reported.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**FIGURE 9** | Relative completion of size growth and biomass production. The top line presents the relative completion of size growth and biomass production during the growing season as well as the time lag between size growth and biomass production (red line). For each species, a reference image and its associated xylem gray-value profile (blue line) is shown. The bottom line presents the variation of the relative completion rate, with the vertical lines highlighting the DOY at which the maximum rate occurs.</figcaption></figure>
interesting insight into the capability of HRXCT to assess biomass production dynamic throughout the course of the growing season. The observed seasonal variation of the time lag between biomass production and size growth agrees with the anatomical signature of the studied species. For both pine and beech, the time lag reaches its maximum during the latewood formation, when 80-90% of the growth ring width is already achieved and thick-walled tracheids or fibers are produced. In oak, the time lag shows the largest values from the end of earlywood formation to the half of the completion of the latewood, thus between 40% and 70% of the growing ring part consists mainly of thick-walled fibers and narrow vessels. Therefore, for all species, the time lag is largest when the denser part of the growth ring is forming and suggests that wall formation and maturation of thick-walled cells require a longer period than the formation and maturation of thin-walled cells. This pattern already observed in conifers (Cuny et al., 2014) seems therefore to apply to angiosperms.

A very short time lag was found (four days) between the maximal size growth rate and the maximal biomass production rate in pine. This result strongly contrasts with the one-month time lag generally acknowledged in conifers (Cuny et al., 2015). However, Andrianantenaina et al. (2019) reported a two weeks lag between biomass production and size growth in spruce growing in north-east France. They speculated that this reduced time lag could be accounted for by the occurrence of a long summer drought during the measurement period in 2015, resulting in the shortening of the growing season and less investment in cell production. The growing season 2019 (when the samples were collected) experienced a very intense drought in late spring-early summer. Also because of the water deficit accumulated during the very dry 2018, the soil in late spring- early summer of 2019 was dry as never observed before in Belgium (Mariën et al., 2021). Therefore, wood formation in black pine might be sensitive to drought as in other conifer species (Andrianantenaina et al., 2019).

At the same site of our sampling, Dox et al. (2021) found that the end of wood formation in oak occurs earlier under drought conditions while this parameter is rather insensitive to drought in beech. In 2019, the authors found that wood formation in oak stopped earlier (DOY = 262) than in 2018 (DOY = 282 in Dox et al., 2021). As the drought in 2019 was even more pronounced than in 2018 (see above), the results corroborate the drought sensitivity of wood formation in oak.

The time lag in beech (two days) was even lower than for pine, while it was much longer in oak (63 days). However, we note that for oak, the peaks of size growth rate and of biomass production rate are much less pronounced than for the other species. During earlywood formation, from late March till early June, both size growth rate and biomass production rate were maximal and rather constant in oak. We believe that the large time lag in oak results from the ring-porous structure of its wood, which implies a fast size growth, due to the enlargement of the vessels, at the beginning of the growing season. Given the uniqueness of such parameters computed for angiosperms, any comparison with previous studies is not possible.

We believe that further development of HRXCT will open the way for more insight into biomass production in angiosperms.

Therefore, we expect that the new HRXCT workflow presented here will be the foundation for major breakthroughs also in the field of intra-seasonal biomass production and carbon sequestration research. Together with a proper conversion of gray value into wood density and species-specific allometric relationships, HXRCT has the clear potential to provide a better quantification of intra-annual dynamics of biomass production and carbon sequestration through high-resolution 3D mapping of carbon within the forming growth ring and among different tree species.

HXRCT is a growing and wide spreading technology in the field of tree- and wood sciences. Both forthcoming technological improvements of scanning devices and protocol optimization (batch preparation and actual scanning of samples), as well as increased computer capacity, will offer potential for the further development of a high-throughput analysis of xylogenesis with HXRCT. Such development will allow to (1) cover a broad range of biomes and climates by increasing drastically the number of sampled individuals and to (2) integrate different tree components (i.e., root, stem, and branches) that may present different phenological patterns and react differently to similar climatic conditions. Thus, the HRXCT approach, developed here, offers a promising avenue to investigate seasonal plant-climate growth relationships spanning a wide range of wood structures, species, biomes, and climatic conditions, and to generate a substantial amount of data that could trigger a new generation of vegetations models fully considering the seasonal wood growth dynamics.

## DATA AVAILABILITY STATEMENT

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

## AUTHOR CONTRIBUTIONS

RL, JVdB, and MC designed the study. BM, MC, and RL collected the samples. JG and PP processed the sample for microtomy. RL and JG conducted the analysis of the histologic slices. RL conducted the X-ray computed tomography analysis and the analysis of the data. RL, JVdB, MC, and HB wrote the first draft of the manuscript. RL, JVdB, MC, and HB wrote the revised manuscript. BM, JG, and PP made significant contributions to the manuscript. All authors gave final approval for publication.

## FUNDING

This research was funded by the BOF Special Research Fund for RL (BOF Starting Grant JVdB, BOFSTG2018000701), the UGCT Center of Expertise (BOF.EXP.2017.0007), the Research Foundation Flanders (FWO grant G009720N), and the LEAF-FALL starting ERC grant (714916).

## ACKNOWLEDGMENTS

We acknowledge the institutions that gave permission to conduct research in the study areas: Agency for Nature and Forest
of the Flemish Government (ANB), the Military Defense of
Belgium, and City of Brasschaat. We would like to thank Stijn
Willen for the technical support and Inge Dox and Sebastien
Leys for their help in collecting and processing the samples
for microtomy.

REFERENCES

Andrianantenaina, A. N., Rathgeber, C. B. K., Pérez-de-Lis, G., Cuny, H., and Ruelle, J. (2019). Quantifying intra-annual dynamics of carbon sequestration in the forming wood: a novel histologic approach. *Ann. Forest Sci.* 76:62. doi: 10.1007/s13595-019-0846-7

Barthelme, S. (2019). *imager: Image Processing Library Based on 'CImg'. R package version 0.41.42.

Bastin, J. F., Fayolle, A., Tarelkin, Y., et al. (2015). Wood specific gravity variations and biomass of central african tree species: the simple choice of the outer wood. *PLoS ONE* 10:e0142146. doi: 10.1371/journal.pone.0142146

Brabant, L., Vlassenbroeck, J., De Witte, Y., et al. (2011). Three-dimensional analysis of high-resolution X-ray computed tomography data with Morpho+. *Microsc. Microanal.* 17, 252–263. doi: 10.1017/S1431927610094389

Bray, D. (2000). "Critical point drying of biological specimens for scanning electron microscopy," in *Supercritical Fluid Methods and Protocols*, eds J. R. Williams, A. A. Clifford (Totowa, NJ: Humana Press), 235–243. doi: 10.1385/1-59259-030-6:235

Brodersen, C. R., and Roddy, A. B. (2016). New frontiers in the three-dimensional visualization of plant structure and function. *Am. J. Bot.* 103, 184–188. doi: 10.3732/ajb.1500532

Camarero, J. J., Olano, J. M., and Parras, A. (2010). Plastic bimodal xylogenesis in conifers from continental Mediterranean climates. *New Phytol.* 185, 471–480. doi: 10.1111/j.1469-8137.2009.03073.x

Campioli, M., Vincke, C., Jonard, M., Kint, V., Demarée, G., and Ponette, Q. (2012). Current status and predicted impact of climate change on forest production and biogeochemistry in the temperate oceanic European zone: review and prospects for Belgium as a case study. *J. For. Res.* 17, 1–18. doi: 10.1007/s10310-011-0255-8

Cufar, K., Prislan, P., de Luis, M., and Gričar, J. (2008). Tree-ring variation, wood formation and phenology of beech (Fagus sylvatica) from a representative site in Slovenia, SE Central Europe. *Trees* 22, 749–758. doi: 10.1007/s00468-008-0235-6

Cuny, H. E., Rathgeber, C. B., Kiesse, T. S., Hartmann, F. P., Barbeito, I., and Fournier, M. (2013). Generalized additive models reveal the intrinsic complexity of wood formation dynamics. *J. Exp. Bot.* 64, 1983–1994. doi: 10.1093/jxb/ert057

Cuny, H. E., Rathgeber, C. B. K., Frank, D., Fonti, P., and Fournier, M. (2014). Kinetics of tracheid development explain conifer tree-ring structure. *New Phytol.* 203, 1231–1241. doi: 10.1111/nph.12871

Cuny, H. E., Rathgeber, C. B. K., Frank, D., Fonti, P., Mäkinen, H., Prislan, P., et al. (2015). Woody biomass production lags stem-girth increase by over one month in coniferous forests. *Nat. Plants* 1:15160. doi: 10.1038/nplants.2015.160

De Baerdemaeker, N. J. F., Stock, M., Van den Bulcke, J., De Baets, B., Van Hoorebeke, L., and Steppe, K. (2019). X-ray microtomography and linear discriminant analysis enable detection of embolism-related acoustic emissions. *Plant Methods* 15:153. doi: 10.1186/s13007-019-0543-4

De Mil, T., Vannoppen, A., Beeckman, H., Van Acker, J., and Van den Bulcke, J. (2016). A field-to-desktop toolchain for X-ray CT densitometry enables tree ring analysis. *Ann. Bot.* 117, 1187–1196. doi: 10.1093/aob/mcw063

Delpierre, N., Vitasse, Y., Chuine, I., Guillemot, J., Bazot, S., Rutishauser, T., et al. (2016). Temperate and boreal forest tree phenology: from organ-scale processes to terrestrial ecosystem models. *Ann. For. Sci.* 73, 5–25. doi: 10.1007/s13595-015-0477-6

Deslauriers, A., Rossi, S., Anfodillo, T., and Saracino, A. (2008). Cambial phenology, wood formation and temperature thresholds in two contrasting years at high altitude in southern Italy. *Tree Physiol.* 28, 863–871. doi: 10.1093/treephys/28.6.863

SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.698640/full#supplementary-material

Dierick, M., Van Loo, D., Masschaele, B., Bulcke, J., Acker, J., Cnudde, V., et al. (2014). Recent micro-CT scanner developments at UGCT. *Nucl. Instrum. Methods Phys. Res. B.* 324, 35–40. doi: 10.1016/j.nimb.2013.10.051

Dox, I., Gričar, J., Marchand, L. J., Leys, S., Zuccarini, P., Géron, C., et al. (2020). Timeline of autumn phenology in temperate deciduous trees. *Tree Physiol.* 40, 1001–1013. doi: 10.1093/treephys/tpaa058

Dox, I., Prislan, P., Gričar, J., Mariën, B., Delpierre, N., Flores, O., et al. (2021). Drought elicits contrasting responses on the autumn dynamics of wood formation in late successional deciduous tree species. *Tree Physiol.* 41, 1171–1185. doi: 10.1093/treephys/tpaa175

Fatichi, S., Leuzinger, S., and Körner, C. (2014). Moving beyond photosynthesis: from carbon source to sink-driven vegetation modeling. *New Phytol.* 201, 1086–1095. doi: 10.1111/nph.12614

Fatichi, S., Pappas, C., Zscheischler, J., and Leuzinger, S. (2019). Modelling carbon sources and sinks in terrestrial vegetation. *New Phytol* 221, 652–668. doi: 10.1111/nph.15451

Friend, A. D., Eckes-Shephard, A. H., Fonti, P., Rademacher, T. T., Rathgeber, C. B. K., Richardson, A. D., et al. (2019). On the need to consider wood formation processes in global vegetation models and a suggested approach. *Ann. For. Sci.* 76:49. doi: 10.1007/s13595-019-0819-x

Gričar, J., Zavadlav, S., Jyske, T., Lavrič, M., Laakso, T., Hafner, P., et al. (2018). Effect of soil water availability on intra-annual xylem and phloem formation and non-structural carbohydrate pools in stem of Quercus pubescens. *Tree Physiology* 39: 222–233. doi: 10.1093/treephys/tpy101

Gutiérrez, Y., Ott, D., Töpperwien, M., Salditt, T., and Scherber, C. (2018). X-ray computed tomography and its potential in ecological research: a review of studies and optimization of specimen preparation. *Ecol. Evol.* 8, 7717–7732. doi: 10.1002/ece3.4149

Hubau, W., Lewis, S. L., Phillips, O. L., Affum-Baffoe, K., Beeckman, H., Cuní-Sanchez, A., et al. (2020). Asynchronous carbon sink saturation in African and Amazonian tropical forests. *Nature* 579, 80–87. doi: 10.1038/s41586-020-2035-0

Knipfer, T., Eustis, A., Brodersen, C., Walker, A. M., and McElrone, A. J. (2015). Grapevine species from varied native habitats exhibit differences in embolism formation/repair associated with leaf gas exchange and root pressure. *Plant Cell Environ.* 38, 1503–1513. doi: 10.1111/pce.12497

Legendre, P. (2018). *Imodel2: Model II Regression. R Package Version 1.7-3.*

Lehnebach, R., Bossu, J., Va, S., Morel, H., Amusant, N., Nicolini, É., et al. (2019). Wood density variations of Legume trees in French Guiana along the shade tolerance continuum: heartwood effects on radial patterns and Gradients. *Forests* 10:80. doi: 10.3390/f10020080

Leroux, O., Leroux, F., Bellefroid, E., Claeys, M., Couvreur, M., Borgonie, G., et al. (2009). A new preparation method to study fresh plant structures with X-ray computed tomography. *J. Microsc.* 233, 1–4. doi: 10.1111/j.1365-2818.2008.03088.x

Marchand, L. J., Dox, I., Gricar, J., Prislan, P., Leys, S., Van den Bulcke, J., et al. (2020). Inter-individual variability in spring phenology of temperate deciduous trees depends on species, tree size and previous year autumn phenology. *Agric. For. Meteorol.* 290:108031. doi: 10.1016/j.agrformet.2020.108031

Mariën, B., Dox, I., De Boeck, H. J., Willems, P., Leys, S., Papadimitriou, D., et al. (2021). Does drought advance the onset of autumn leaf senescence in temperate deciduous forest trees? *Biogeosciences* 18, 3309–3330. doi: 10.5194/bg-18-3309-2021

Ooms, J. (2020). *magick: Advanced Graphics and Image-Processing in R. R Package Version 2.3.*

Paganin, D., Mayo, S. C., Gureyev, T. E., Miller, P. R., and Wilkins, S. W. (2002). Simultaneous phase and amplitude extraction from a
single defocused image of a homogeneous object. *J. Microsc.* 206, 33–40.
doi: 10.1046/j.1365-2818.2002.01010.x

Prislan, P., Cufar, K., and De Luis, M., Gričar, J. (2018). Precipitation is not limiting for xylem formation dynamics and vessel development in European beech from two temperate forest sites. *Tree Physiol.* 38, 186–197.
doi: 10.1093/treephys/tpx167

Prislan, P., Gričar, J., and Cufar, K. (2014). *Wood Sample Preparation for Microscopic Analysis - Protocol of the sTReESS COST Action FP1106.*
Available online at: http://www.streess-cost.eu/images/stories/Documents/protocol_wood_sample_preparation_for_microscopic_analysis.pdf

Prislan, P., Gričar, J., de Luis, M., and Smith, K. T., Cufar, K. (2013).
Phenological variation in xylem and phloem formation in *Fagus
sylvatica* from two contrasting sites. *Agric. For. Meteorol.* 180, 142–151.
doi: 10.1016/j.agrformet.2013.06.001

Pya, N., and Wood, S. N. (2015). Shape constrained additive models. *Stat. Comput.* 25, 543–559. doi: 10.1007/s.11222-013-9448-7

R Development Core Team (2019). *R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing.

Rasband, W. S. (2012). *Image*. Bethesda, MD: U. S. National Institutes of Health.

Rathgeber, C. (2017). Conifer tree-ring density interannual variability—anatomical, physiological and environmental determinants. *New Phytol.* 216, 621–625. doi: 10.1111/nph.14763

Rossi, S., Anfodillo, T., and Menardi, R. (2006). TrepHor: a new tool for sampling microcores from tree stems. *IAWA J.* 27:89. doi: 10.1163/22941932-900 00139

Rossi, S., Deslauriers, A., and Morin, H. (2003). Application of the Gompertz equation for the study of xylem cell development. *Dendrochronologia* 21, 33–39.
doi: 10.1078/1125-7865-00034

Rossi, S., Rathgeber, C. B. K., and Deslauriers, A. (2009). Comparing needle and shoot phenology with xylem development on three conifer species in Italy. *Ann. For. Sci.* 66, 206. doi: 10.1051/forest/2008088

Simpson, G. L. (2020). gratia: Graceful 'ggplot'-Based Graphics and Other Functions for GAMs Fitted Using 'mgcv'. *R Package Version* 0.3.0.

Stuppy, W. H., Maisano, J. A., Colbert, M. W., Rudall, P. J., and Rowe, T.
B. (2003). Three-dimensional analysis of plant structure using high-
resolution X-ray computed tomography. *Trends Plant Sci.* 8, 2–6.
doi: 10.1016/S1360-1385(02)00004-3

Van den Bulcke, J., Boone, M., Dhaene, J., Loo, D., Hoorebeke, L., Wyffels, F., et al. (2019). Advanced X-ray CT scanning can boost tree-ring research for earth-system sciences. *Ann. Bot.* 124, 837–847. doi: 10.1093/aob/mcz126

Vlassenbroeck, J., Dierick, M., Masschaele, B., Cnudde, V., Van Hoorebeke, L., and Jacobs, P. (2007). Software tools for quantification of X-ray microtomography at the UGCT. *Nucl. Instrum. Methods Phys. Res. B.* 580, 442–445. doi: 10.1016/j.nima.2007.05.073

Wood, S. N. (2017). *Generalized Additive Models: An Introduction with R*. Boca Raton, FL: Chapman and Hall/CRC. doi: 10.1201/9781315370279

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

**Publisher's Note:** All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Copyright © 2021 Lehnebach, Campioli, Gričar, Prislan, Mariën, Beeckman and Van den Bulcke. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.