Understanding the Formation of Heartwood in Larch Using Synchrotron Infrared Imaging Combined With Multivariate Analysis and Atomic Force Microscope Infrared Spectroscopy

Formation of extractive-rich heartwood is a process in live trees that make them and the wood obtained from them more resistant to fungal degradation. Despite the importance of this natural mechanism, little is known about the deposition pathways and cellular level distribution of extractives. Here we follow heartwood formation in Larix gmelinii var. Japonica by use of synchrotron infrared images analyzed by the unmixing method Multivariate Curve Resolution – Alternating Least Squares (MCR-ALS). A subset of the specimens was also analyzed using atomic force microscopy infrared spectroscopy. The main spectral changes observed in the transition zone when going from sapwood to heartwood was a decrease in the intensity of a peak at approximately 1660 cm-1 and an increase in a peak at approximately 1640 cm-1. There are several possible interpretations of this observation. One possibility that is supported by the MCR-ALS unmixing is that heartwood formation in larch is a type II or Juglans-type of heartwood formation, where phenolic precursors to extractives accumulate in the sapwood rays. They are then oxidized and/or condensed in the transition zone and spread to the neighboring cells in the heartwood.


INTRODUCTION
Heartwood (HW) formation is the final step in the life cycle of ray cells. Before cell death, ray cells undergo metabolic changes in the transition zone between sapwood (SW) and HW, resulting in increased synthesis of secondary metabolic compounds called extractives. The extractives have a significant effect on the properties of wood, most notably regarding its resistance to fungal decay and other forms of biological attack (Hillis, 1987;Hinterstoisser et al., 2000;Schultz and Nicholas, 2000;Taylor et al., 2002). In order to understand how extractives contribute to wood durability, many studies have focused on the chemical interplay between extractives and decay agents (Valette et al., 2017). However, there is an increasing understanding that the distribution of extractives within cell walls also plays an important role, though hitherto under-investigated (Kampe and Magel, 2013).
The process of HW formation has been studied for decades and is known to be associated with parenchyma cell death, disappearance of storage material, and increase in extractives content (Hillis, 1985;Hillis, 1987). Kampe and Magel (2013) described two possible HW formation mechanisms: Type I or Robinia-Type proposes the accumulation of the phenolic extractives in the transition zone without any indication of phenolic precursors in the aging SW (Nair et al., 1981;Magel et al., 1994;Bergström et al., 1999); Type II or Juglans-Type of HW formation suggests a gradual accumulation of phenolic precursors in the aging SW tissues. In type II, HW extractives are formed in the transition zone by primary and secondary reactions, such as oxidation and hydrolysis of precursor substances (Dellus et al., 1997;Burtin et al., 1998;Mayer et al., 2006). Once the extractives are formed, they are released into the lumina of neighboring cells and cell walls (Déjardin et al., 2010;Kampe and Magel, 2013). When inside the cell walls, a few studies suggest that at least some extractives are covalently bound to the structural cell wall polymers through enzymatic activity (Monties, 1991).
The genus Larix species (larch) are an important European resource for durable wood (Hillis, 1987). The extractives in larch belong to the molecular families of terpenoids, flavonoids, lignans, fatty acids, and galactans (Zule et al., 2015;Zule et al., 2017). Like all conifers, larch trees contain high amounts of oleoresin, produced by specialized epithelial cells surrounding resin canals and ray parenchyma cells. The resin is composed of fatty acids and esters thereof, as well as various subgroups of terpenoids, and is distributed throughout the HW and SW through the network of resin canals and ray cells (Hillis, 1987). It was shown for Pinus sylvestris that in HW, the composition of resin is enriched with phenolic compounds, presumably produced by ray parenchyma cells during HW formation (Felhofer et al., 2018).
Within the large family of terpenoids, diterpenoid acids (called resin acids) constitute the largest part of the oleoresin (Higuchi, 1997), and they have been repeatedly shown to have fungicidal properties (Keeling and Bohlmann, 2006), which is also the case for triterpenes, known as sterols (Burčová et al., 2018). Triglycerides and fatty acids may have a role in moisture regulation (Tomppo et al., 2011), which is important in the context of degradation by microorganisms. The principal phenolic compounds detectable in larch HW are flavonoids, the main compounds being taxifolin (C 15 H 12 O 7 ) and dihydrokaempferol (C 15 H 12 O 6 ). Minor amounts of lignans can also be found. Flavonoids have been attributed with fungicidal properties, but their main potential seems to be their ability to scavenge different types of radicals, as well as to reduce and chelate metals (Giwa and Swan, 1975;Cao et al., 1997;Babkin et al., 2001;Ivanova et al., 2012;Zule et al., 2017).
Larch HW is appreciated for its good mechanical properties, its color, and specially for its natural durability (Gierlinger et al., 2004). A strong relationship between extractives content and brown-rot decay resistance has been shown (Gierlinger et al., 2002;Windeisen et al., 2002). Nevertheless, very little is known about the formation and distribution of larch extractives within the xylem tissue at the cell and cell wall level. Their micro and nano-scale distribution is of importance (Taylor et al., 2002) since extractives are more effective against wood degradation within cell walls than in extracellular voids (Hillis, 1987). To investigate extractives in context with the microstructure, TOF-SIMS imaging has been applied in Cryptomeria japonica trees and showed that the extractives tend to accumulate near radial rays (Imai et al., 2005;Saito et al., 2008). Recently, the potential of Confocal Raman Microscopy to follow the extractive distribution in sapwood (SW) and heartwood (HW) of Scots pine (Pinus sylvestris, a moderately durable species) was shown (Belt et al., 2017;Felhofer et al., 2018). On the micro-level, pinosylvins were reported in the lumen, as well as in the compound middle lamella (CML), cell corner (CC), and pits of tracheid cells.
Synchrotron Radiation Fourier Transform Infrared (SR-FTIR) imaging is the ideal technique to study the extractive deposition patterns at the microscale during HW formation in larch because of the high brightness and high collimation of the beam and avoidance of the fluorescent problems experienced with Raman microspectroscopy. SR-FTIR imaging provides spatial and spectral information about the samples and, therefore, informs on the composition and location of the different sample constituents. Despite the relevant information contained in SR-FTIR images, the analysis of this kind of measurement is not straightforward because of the often large image sizes and the mixed signal components present in the spectra collected. To help in the signal-unmixing task, multivariate analysis tools like Multivariate Curve Resolution -Alternating Least Squares (MCR-ALS) are used. Indeed, MCR-ALS has already been proven to adapt particularly well to hyperspectral image analysis because of the easy introduction of external spectral and spatial information about the image in the analysis and the ability to work with both single and multiset (several) image structures (Tauler et al., 1995;de Juan et al., 2004;Piqueras et al., 2014). This approach is the main tool used in the current study to obtain distribution maps and spectral signatures of the wood sample (Felten et al., 2015;Piqueras et al., 2015).
To supplement the SR-FTIR images, atomic force microscopy infrared spectroscopy (AFM-IR) was used. AFM-IR combines atomic force microscopy (AFM) with pulsed IR laser ( Figure S1) to obtain localized mid-IR spectra (3600-900 cm -1 ) of regions as small as tens of nm in the horizontal plane and with a vertical resolution of~0.1 nm (Dazzi and Prater, 2017). Such resolution surpasses the resolution of optical IR instruments making AFM-IR a suitable technique to study nanoscale properties of wood materials. Within the study of plants, AFM-IR has been used to analyze the composition of thylakoids (Janik et al., 2013) and of epicuticular wax (Farber et al., 2019), and to understand how the structure and composition of the Populus nigra cell wall affects water transport within the xylem (Pereira et al., 2018). Further, it has been used to identify the products of the reaction between the cell wall of Pinus taeda and a phenol-formaldehyde resin (Wang et al., 2005). However, to our knowledge, no one has yet studied the nanoscale compositional variations between the cell wall, the middle lamella, and the rays of SW and HW.
The objective of this work was to obtain a detailed overview of the HW formation process in larch (Larch gmelinii i.e. var Japonica) at the micro and nanoscale by combining SR-FTIR, AFM-IR, and advanced chemometric tools.

Sample Preparation
For this imaging study, a sample of Larix gmelinii var japonica (Kurile larch) was taken at 1.3 m stem height. The tree was felled in October 2017 near Hørsholm, north of Copenhagen, Denmark. The wood was freeze dried to avoid possible artifacts from air drying. From an area including nine annual rings, nine tangential (D T1 -D T9 ) and nine cross-sections (D X1 -D X9 ) were cut without any embedding of the samples, in order to preserve the major content and distribution of extractives. Samples of nominally 10 µm thickness were obtained using a Leica microtome (Leica RM2255) (Figure 1). For transportation, the samples were placed on a glass slide and covered with a glass coverslip.

Synchrotron Infrared Imaging
All the tangential and cross-sections were imaged at the IR beamline MIRAS of the ALBA synchrotron (Cerdanyola del Vallés, Spain, proposal 2018022761). Before IR imaging, samples were inspected by light microscopy in order to select regions of interest (ROI) that included ray and tracheid cells and appeared to have a sample transparency that would allow SR-FTIR measurements in transmission mode. After area selection, the sample was carefully transferred onto a ZnSe disc of 1 mm thickness. The sample edges were fixed to the disc with tape to avoid sample movement during imaging. The IR measurements were acquired with a Bruker system (Hyperion 3000 microscope coupled to a Vertex 70 spectrometer) equipped with a liquidnitrogen cooled mercury cadmium telluride (MCT) detector.
Preliminary tests were carried out prior to the IR imaging. During these tests, possible effects of the high-energy IR beam on the wood material was studied by collecting punctual IR spectra at different exposure times. No adverse effects on the tissue or the spectra were observed for any of the exposure times tested. Therefore, we are confident that the energy of the IR beamline did not adversely influence our SR-FTIR measurements.
The IR images were acquired in transmission mode, using a 36x objective. The images were collected with a 3 x 3 µm 2 spatial resolution. All spectra were obtained in the infrared region (4000 −800 cm −1 ) with 64 co-added scans. Absorbance representation was used throughout. A total of nine SR-FTIR images were acquired in both tangential and cross directions.

Atomic Force Microscopy Infrared Spectroscopy
AFM-IR exploits the photothermally induced resonance effect to detect the absorption of IR radiation with the AFM tip (Dazzi et al., 2005). In short, the sample is irradiated with the IR source and mechanically expands to dissipate the absorbed energy ( Figure S1). The strongest expansion happens when the sample is irradiated with the IR wavelength that corresponds to the maximum absorption by the sample. Thus, by placing the AFM tip directly above the irradiated area, it is possible to detect the expansion of a sample by monitoring the deflection and oscillation of the AFM cantilever. This thermal expansion is directly proportional to the absorption coefficient of the excited area (Dazzi et al., 2012). By analyzing the many frequencies and amplitudes of the resonant oscillations of the AFM cantilever with Fourier transform, we extracted the useful information and reconstructed the IR spectra of various regions on the sample with a spatial resolution that is close to the size of an AFM tip.
We used a nanoIR from Analysis Instruments, Inc., to obtain mid-IR spectra of the ray, the middle lamella, and the secondary cell wall. Only two of the samples of the cross-sections (D X1 and D X8 ) were used for the AFM-IR investigation. We fixed the edges of a sample with adhesive tape to a glass slide and acquired AFM images in contact mode. First, we found a suitable region where we could see both the middle lamella and the tracheid cell wall, or the ray and the adjacent tracheid cell wall and imaged it. Then, we collected and averaged three-background IR spectra with the resolution of 4 cm -1 to account for the variations in power of the IR source. For the IR spectral acquisition, we chose the second mode of the cantilever vibration to record the signal. This mode had a frequency of about 190 kHz and we chose the frequency window to be ±25 kHz to account for the variations in thermoelastic properties between the cell walls and the middle lamella or the ray. The second mode was chosen to improve the signal to noise ratio of the cantilever amplitude. For AFM-IR, it is crucial to position the IR laser directly at the tip-sample contact to make sure that the recorded spectrum originates directly from the area at the sample above where the tip is positioned. To do so, we scanned around the tip-sample contact area with the IR laser to find the highest cantilever amplitude, which corresponds to the position where the IR laser has maximum power. Once the tip, sample, and laser positions were optimized, we collected the spectra at various positions of the sample with the energy resolution of 4 cm -1 and by co-averaging 256 scans. After IR spectral acquisition, the same area was inspected by using the build-in camera. No laser damage was observed on any of the samples. The AFM images were flattened to remove the tilt and the AFM-IR spectra were smoothed by a Savitzky-Golay filter (2 nd polynomial degree, 15 points window size) (Savitzky and Golay, 1964).

SR-FTIR DATA TREATMENT
The data treatment of SR-FTIR images consisted of two consecutive steps: (1) preprocessing of the image spectra to correct for scattering effects and (2) analysis by Multivariate Curve Resolution -Alternating Least Squares (MCR-ALS) (Tauler et al., 1995;de Juan and Tauler, 2006) to obtain pure spectra of the image constituents and their related distribution maps. The next subsections describe these steps.

Data Preprocessing
Due to the thickness and density of the samples, the infrared spectra of ray and tracheid cells were oversaturated in both the low and the high spectral wavelength range. Therefore, these spectral areas had to be excluded, and only the 1200-1750 cm -1 range was included in the analysis.
Infrared spectra are prone to artifacts because of Mie scattering associated with surface irregularities. Such artifacts may produce a broad oscillation in the baseline spectrum and can lead to distortions in both the position and intensity of absorption bands (Romeo et al., 2006). The raw data were corrected by the algorithm Asymmetric Least Squares (AsLS) (Eilers, 2004), which has been demonstrated to cope well with this type of scattering (Piqueras et al., 2013).

Hyperspectral Image Resolution
The goal of hyperspectral image resolution and, consequently, of the multivariate curve resolution alternating least squares (MCR-ALS) algorithm, is the decomposition of the original raw image data into distribution maps and pure spectra of the constituents present in the imaged sample (Tauler et al., 1995;Jaumot et al., 2005;de Juan and Tauler, 2006;Tauler et al., 2009). In the matrix form, hyperspectral images can be well described by a bilinear model based on the Beer-Lambert law [Eq.(1)], where the D matrix contains the original raw spectra, which are decomposed into a set of concentration profiles (C matrix) and corresponding pure spectra (S T matrix) of the constituents present in the image. Every row of the S T matrix corresponds to the pure spectrum of an image constituent, while every column of the C matrix of concentration profiles corresponds to the related pixel-to-pixel variation of its chemical concentration. It should be pointed out that each column of the C matrix can be refolded appropriately in order to recover the original two-dimensional spatial image structure and then pure distribution maps are obtained.
MCR-ALS is a flexible method that allows analyzing a single image individually or several images simultaneously. To obtain a complete and reliable description of the HW formation, the simultaneous analysis of the recorded images along the nineannual rings obtained from larch sample was performed. Thus, two multisets were built and analyzed separately-formed by nine and eight images collected in the cross-sectional (D X ) and tangential directions (D T ), respectively. See Figure 2 for a visualization of the image multiset structure. In these multisets, the spectral dimension of all images is the same, while the image dimensions may differ between images because the images are unfolded before being merged into a single matrix (Figure 2). Due to extreme cases of over-saturation, the first image of the cross-sections (D T1 ) and the last image of the tangential sections (D X9 ) had to be excluded from the multiset structures D T and D X , i.e: A multiset structure also follows the bilinear model based on Beer-Lambert law [see Eq. (1)]. In this example, image multiset analysis by MCR-ALS provides a single matrix S T of pure spectra, identical for all the images analyzed, and a C matrix formed by as many submatrices as the number of images included in the data set. Every column of each C submatrix can be refolded conveniently to recover the distribution map of each constituent present in the different images of the data set ( Figure 2).
MCR-ALS multiset analysis was performed on both multiset structures D X and D T following the MCR-ALS steps described in the literature (Jaumot et al., 2005). The first step consisted of determining the number of components involved during HW formation by singular value decomposition of the whole preprocessed D X and D T matrices (Golub and Reinsch, 1970). Five contributions were needed to describe the variation in both multisets. Then, initial estimates of pure spectra were obtained with a method based on SIMPLISMA (Windig and Guilment, 1991). The spectral estimates and the original multiset are used to perform an iterative alternating least squares optimization of matrices C and S T under constraints. To obtain unmixed resolved profiles that are chemically meaningful, the constraints used in the resolution of both multiset structures were nonnegativity in both, the concentration and the spectral profiles (Bro and de Jong, 1997), and normalization of pure spectra in the S T matrix (using 2-norm, i.e., the Euclidean norm).
After a preliminary MCR-ALS analysis of both multisets (D X and D T ) under non-negativity constraints, we realized by inspection of the distribution maps obtained that there were components absent in some sub-images of both multisets. As a consequence, a new MCR-ALS analysis was performed to obtain more accurate solutions by imposing the additional constraint of correspondence of species, which encodes information on the presence/absence of constituents in the concerned images of the multisets (Tauler et al., 2009), i.e. when a certain constituent is absent in one image, the related concentration profile is null. The lack of fit was 6.37% and 7.44% for D X and D T , respectively, which is satisfactory for FTIR measurements (Felten et al., 2015).
It is important to note that, when resolving images of biological samples, each resolved contribution (component) may refer to a mixture of chemical compounds of defined composition (polysaccharides, sugars, polymers, fatty acids, flavonoids…) that are present in a particular location of the sample and often represent a distinct biological element, e.g., a tissue or a cell part. This means that an MCR contribution is not necessarily a pure chemical compound. In imaging, the component discriminating ability of MCR will also depend on FIGURE 2 | MCR-ALS analysis for a multiset structure of images, where x and y are spatial pixels and l represents the wavenumbers of the spectra. Dx n is the augmented data matrix; C is the concentration profiles; and S T is the pure spectra matrix.
how fine the spatial resolution of the imaging technique used is, e.g., if two components are differently distributed at nanoscale level, MCR will not resolve them if imaging is performed at microscale level.

Exploratory Analysis of the SR-FTIR Images
In order to identify the main spectral variations during the HW formation of Kurile larch, an exploratory analysis was done. A small area of the ray and tracheid cell wall was selected for each of the images collected in the cross-sectional direction ( Figure  3A). Figure 3B shows the average spectra of the ray area selected for each of the collected images (the average spectra of the tracheid cell wall area can be found in the supplementary information ( Figure S2). The main spectral features when going from the SW to the HW of the ring area are the emergence of a new band at 1640 cm -1 and the decrease in intensity of the band at 1660 cm -1 . These spectral features appear in the sample D X5 , which we therefore determined to be the transition zone between SW and HW ( Figure 3C). According to the literature, the band at 1640 cm -1 could correspond to adsorbed water (Popescu et al., 2007) or to a carbonyl stretching (C = O) as a consequence of the presence of para substituted ketone or aryl aldehydes (Lawther et al., 1996;Shi et al., 2012). Since the OH band at 3360 cm -1 appears not to covary with the band at 1640 cm -1 , we find it unlikely that the latter is related to water adsorption.

Resolution of Image Multiset Structures. MCR-ALS on Complete SR-FTIR Images
MCR-ALS multiset analysis on complete images provides the biological spectral signatures and distribution maps needed to integrally describe the ray and tracheid cells across the transition zone during the HW formation of larch. Although, as mentioned before, there is not necessarily a one-to-one correspondence between the MCR-ALS contribution and the individual chemical FIGURE 3 | (A) Representation of the area selection of the lumen, cell wall, and ray. (B) Average spectra of the ray area selected for each of the cross-section images collected across the heartwood formation zone. (C) Zoom of the spectral range from 1200 cm -1 to 1800 cm -1 of the average spectra of the ray area selected for each of the cross-section images. The spectra gradually change from blue (sapwood) to red (heartwood) color.
compounds, each contribution can be associated with some particular kind(s) of plant cell region, because of the morphology of the distribution maps and the spectroscopic features found in the resolved spectra. The resolution results for the multiset of the cross sections (D X ) are shown in Figure 4. The multiset analysis of the tangential section images (D T ) gave similar MCR-ALS results as D X and are shown in Figure S3. The light microscopy images, corresponding to the imaged ray and the two surrounding tracheid cells areas are shown in the left of Figure 4A. The resolved distribution maps of each image are displayed in the right side of Figure 4A and the related resolved spectra in Figure 4B.
MCR-ALS was able to resolve the main plant cell constituents of the wooden tissue and to give further details on chemical changes occurring in the transition from SW to HW. Component I is characterized by signals in the region between 1317-1370 cm -1 , which are mainly associated with cellulose (Colom et al., 2003;Popescu et al., 2010) (see Table 1). The corresponding distribution maps show higher intensity in what we identify as the secondary cell wall (S2) of the tracheids, known to be thick in latewood and rich in cellulose (Fengel and Wegener, 1989). This component was distributed evenly across all the growth rings, as is also seen clearly in Figure S3A of the tangential sections.
The most prominent bands of lignin at 1505 and 1610 cm -1 are associated with C = C stretching of the aromatic ring modes (Colom et al., 2003;Weiland and Guyonnet, 2003;Popescu et al., 2010;Gorsás et al., 2011). They appear in components I, III, IV, and V, but show higher intensity and thus higher lignin contribution in components IV and V. Components IV and V appear to be situated in the same anatomical segments, i.e. in the cell corners (CC) and in the compound middle lamella (CML; middle lamella + adjacent primary walls), which is consistent with previous studies showing high lignin concentration in these locations (Fengel and Wegener, 1989). The main spectral difference between component IV and V is the same spectral variation that was observed during the exploratory analysis of SR-FTIR images: Component IV shows the characteristic band at 1660 cm -1 , which is attributed to the ethylenic C = C (in coniferyl alcohol/sinapyl alcohol units) and C = O (in coniferaldehyde/ sinapaldehyde) bond stretches of lignin (Umesh and Rajai, 2010;Umesh et al., 2011;Bock and Gierlinger, 2019). This component was prevalent in the CC and CML of SW tracheid cells but disappears in the HW. The other feature observed in the preliminary analysis, the band at 1640 cm -1 , only appears in component V, which is distributed in the CC, CML, and part of the ray area of HW tracheids, as is also observed in Figure S3A. According to the literature, the IR band at 1640 cm -1 is assigned to a carbonyl stretching due to para substituted ketone or aryl aldehydes (Lawther et al., 1996;Shi et al., 2012). This band was also assigned to hydrogen bonding to the carbonyl group, as reported elsewhere (Pomar et al., 2002;Agarwal and Reiner, 2009;Bock and Gierlinger, 2019). Because component IV is present from D X1 to D X5 (SW region), and component V appears from D X5 to D X8 (HW region), it can be deducted that this fifth annual ring represents the transition zone, where the process of HW formation starts.
The resolved IR spectrum for component III shows a distinct, broad band at 1700-1736 cm -1 and is assumed to be formed by the overlapping of the C = O stretch vibration of acetyl or carboxylic acid (COOH) groups (Faix, 1991;Schwanninger et al., 2004) and the C = O stretch vibration of unconjugated ketones, carbonyls, and esters groups. This broader band centered at 1728 cm -1 could suggest the presence of resin acids in the rays since C = O stretching belonging to the -COOH group absorb around 1700 cm -1 (Traoré et al., 2018). Besides, bands are observed around 1600, 1637, and 1660 cm -1 ; the same ones as found in components IV and V, but weaker. Candidates for para substituted ketones are flavonoids such as taxifolin and dihydrokaempferol (Ruddick and Xie, 1994;Kocábová et al., 2016), known to be abundant in larch (Nisula, 2018). This component is represented in all the images in the ray, as well as tracheid lumina and S3 layer.
As mentioned before, resolved IR spectra reflect a mixture of different kinds of biomolecules. For example, the IR spectrum related to component IV consists of a mixture of mainly lignin and hemicelluloses, known by the presence of the band at 1738 cm -1 , which is attributed to ester carbonyl groups prominent in hemicelluloses (Stewart et al., 1995;Gorsás et al., 2011) (see Table 1 for the different IR bands of wood assembled from the literature). Finally, since the samples were measured in the dry state for the SR-FTIR experiment, the IR spectrum of component II associated with part of the lumen is not shown, since it does not contain any biologically relevant information. The spectrum corresponds to the IR absorbance of the ZnSe slide used in the measurements.

AFM-IR Spectra Analysis
The AFM-IR spectra collected in the tracheid cell in the crosssectional sections of SW (D X1 ) and HW (D X8 ) are shown in Figure 5, together with the AFM images. Where possible, we acquired the AFM image of the middle lamella and cell wall in a single image ( Figure 5A). However, in some places, the middle lamella (and the ray) contain topographical features that exceed a few micrometers. Such complex topography is difficult to image in contact mode because of the limited vertical range of the AFM scanner. In addition, the imaging of such complex topographical features blunts the AFM tip rapidly. Hence, we took images of various dimensions, depending on the topography of the area, in order to minimize the mechanical strain exerted on the AFM tip, but to still be able to get a view of both the cell wall and the middle lamella in the same image ( Figure 5B). Numbered crosses on the AFM images indicate the location of AFM-IR spectra. In contrast to the SR-FTIR spectra, the spectral region between 1000-1200 cm -1 is not saturated in the AFM-IR measurements, so we were able to obtain information from that region as well. Although the IR spectra of CML (point nr. 2) and S2 layer (point nr. 1) are very similar to each other in the SW ( Figure 5A), we can observe different intensities around 1504-1510 cm -1 , associated with C = C stretching of the aromatic ring modes of lignin. This band is more intense in the CML, which was also seen in the SR-FTIR data and is consistent with literature (Fengel and Wegener, 1989). The C = O stretching mode is slightly shifted to lower wavelengths (1724 cm -1 ) in point nr. 3, which is characteristic of pectins and/or hemicelluloses (Gorsás et al., 2011). Cellulose bands dominate the IR spectra in the S2 layer (see Table 1). In the case of the HW tracheids ( Figure 5B), no spectral differences between the CML (point nr. 2) and S2 layer (point nr. 1) were observed. By comparing AFM-IR spectra of Figures 5A (SW) and B (HW), we see the presence of the band at 1648 cm -1 , while the band at 1660 cm -1 is missing in the HW tracheid cell wall; the opposite is the case for SW. These are the same main spectroscopic features that were found in the analysis of SR-FTIR images. Finally, a new band at 1108 cm -1 appears in HW ( Figure 5B), linked with COH in plane deformation of celluloses and hemicelluloses  and/or with aromatic C-H in plane deformation of lignin. Figure 6 shows the AFM images and AFM-IR spectra of the ray region and the S2 cell wall of an adjacent tracheid in SW ( Figure 6A) and HW ( Figure 6B). Imaging of the ray and the cell wall in contact mode was particularly difficult because of their different material properties. This is why the image is blurry and  Schwanninger et al., 2004;b Fackler et al., 2010;c Ghosh et al., 2015;d Popescu et al., 2010;e Popescu et al., 2007;f Zu et al., 2012;g Faix, 1991;h Traoré et al., 2018;i Lawther et al., 1996;j Kocábová et al., 2016;k Bock and Gierlinger, 2019. it appears as if the ray is smeared over the cell wall ( Figure 6). Such topography and relationship between the ray and the cell wall are unrealistic and simply an artefact of imaging in contact mode. This artefact does not affect the acquisition of AFM-IR spectra or its spectral features because the spectra are acquired after imaging was finished and at the frequency characteristic for the tip-substrate system. Lignin and cellulose bands are more intense in the S2 cell wall (points nr. 3, 8, and 9) of SW compared to the ray area (points nr. 1, 2, 4, 5, 6, and 7). A characteristic peak occurs at 1076 cm -1 inside the ray, assigned to C-O bands in primary and secondary alcoholic groups (Compounds and Hergert, 1960). A low intensity of the band at 1728 cm -1 is seen in the S2 cell wall. The ray region spectrum of HW (points nr.1 and 2) ( Figure 6B) reveals characteristics peaks at 1692, 1756, and 1764 cm -1 , which represent the C = O stretching in conjugated ketones (Popescu et al., 2007) and carboxylic acid groups , alkyl esters (including the methyl ester of fatty acids), and in g-lactone (Lievens et al., 2011), respectively. The band at 1164 cm -1 is typical for 5,7-dihydroxysubstituted flavonoids (Zu et al., 2012) and indicates the presence of taxifolin. It is also important to highlight the presence of the band at 1072 cm -1 , associated with C-O deformation in primary and secondary alcohol groups of galactosyl subunits (Ghosh et al., 2015). In the S2 layer (point nr. 3), the AFM-IR spectrum shows more intense cellulose/hemicellulose bands, whereas the band at 1660 cm -1 is shifted to higher wavenumbers compared to the S2 layer in SW, likely due to the formation of C = O conjugated ketones. Nevertheless, the lignin band is again more intense and a shoulder at 1640-1660 cm -1 appears.

DISCUSSION
The formation of HW is linked to the occurrence of nonstructural substances called extractives, which play an important role in the resistance of wood to fungal decay (Hinterstoisser et al., 2000;Schultz and Nicholas, 2000). By combining high resolution SR-FTIR with the powerful unmixing algorithm MCR-ALS, we were able to identify a component associated mainly with phenolic compounds and likely with deposition of resin acids (component III from D X and  Figures 4 and S2). Component III shows prominent IR bands at 1637, 1658, and 1728 cm -1 in the ray, tracheid lumen, and S3 layer. The band around 1637 cm -1 was assigned to the ketone bond in taxifolin (Ruddick and Xie, 1994;Miklečić et al., 2012;Kocábová et al., 2016;Liu et al., 2018), one of the most abundant phenolic compounds in larch wood (Giwa and Swan, 1975;Babkin et al., 2001;Ivanova et al., 2012;Zule et al., 2017) (see reference spectrum of taxifolin in Figure S4 for comparison). The other peaks seen in the taxifolin spectrum are not seen due to overlap with the spectra of the structural wood cell wall biopolymers. The distinct, broad band 1700-1736 cm -1 centered at 1728 cm -1 , is assigned to the C = O vibration of carboxylic acid groups in resin acids (Faix, 1991;Schwanninger et al., 2004). It is important to emphasize that component III appears to be more common in the ray than in tracheid cells in the SW and vice versa in the HW (see Figures 4B and S3). As described by Hillis (1987), parenchyma cells die when HW forms and the polyphenols diffuse into cell walls. The relatively higher concentration of taxifolin and resin acid mixtures in the HW cell walls likely explains its natural durability.
With the MCR-ALS multiset analysis of tangential and crosssections, two MCR-ALS contributions were found to be directly linked to the process of HW formation in Kurile larch (components IV and V, Figures 4 and S3). Both of them mainly showed lignin bands but had different spatio-temporal distributions, as well as different spectral features. Component IV was located in the CC and CML of tracheid cells in the SW, although in the transition zone it was almost exclusively present in the ray area. The other lignin contribution (Component V) was distributed in CC, CML, and in the ray cells of HW. The main spectroscopic difference between component IV and V was the appearance of the band at 1640 cm -1 and the simultaneous disappearance of the band at 1660 cm -1 .The disappearance of the 1660 cm -1 band might be explained by a decrease in the intensity or as a drastic shift in the frequency to 1640 cm -1 .
The decrease of the intensity of band at 1655-1660 cm -1 could be explained in terms of condensation reactions (Yamauchi et al., 2005) of coniferyl alcohol in lignin. By condensation, coniferyl alcohol loses the ethylenic bond which then no longer contributes to the vibration at 1660 cm -1 . However, this process cannot explain the appearance of the 1640 cm -1 band. Possible coupling of the condensation reaction with the integration of new coniferyl aldehyde moieties into the lignin structure and their H-bonding to unreacted alcohols would decrease the frequency of the 1660 cm -1 band to 1640 cm -1 , as described by Bock and Gierlinger (2019) and Agarwal and Reiner (2009). Another possibility for the band shift is the oxidation of coniferyl alcohol to its aldehyde and subsequent H-bonding to the carbonyl groups during the aging of the wood cells. Lastly, the appearance of taxifolin and other flavonoids may explain the appearance of the 1640 cm -1 band, as they contain carbonyl vibrations at 1640 cm -1 . Since flavonoids are rich in OH-groups, they are likely to interact with lignin and cause the shift of the 1660 cm -1 band frequency.
From the MCR-ALS distribution maps, interpretation of the deposition pathways of the components involved during HW formation of larch was achieved. It is interesting to observe that distribution of extractives (component III from D X and D T , see Figures 4A and S3A) is represented in all the images and can be found in the ray region, lumen, and S3 layer of tracheid cells. It may suggest that precursor molecules are present before and after the actual transition from SW to HW. This pattern is observed in Juglans-Type II HW formation. Not much is known about this mechanism, but it has been described for deciduous trees, as well as conifers, such as Prunus, Platycarya, Eucalyptus, and Pseudotsuga (Kampe and Magel, 2013). If we follow their deposition in the distribution maps of Figures 4A and S3A, it seems that extractives are accumulated in the ray in SW and after ray cell death, they spread to the surrounding wood tissues. We can also see how component V emerges and component IV diminishes during the HW formation process. The similar IR spectra and spatial deposition indicate that component IV corresponds to a set of precursor molecules of component V. Thus, Type II seems the reasonable mechanism of larch HW formation.
AFM-IR spectroscopy revealed variations in the composition between the tracheid and the ray regions in SW and HW ( Figures 5A, B). The main spectral differences were found at 1108, 1456, 1648, and 1660 cm -1 . In the HW tracheid cell walls, the band at 1108 cm -1 , assigned to COH in plane deformation of celluloses and hemicelluloses  and/or to aromatic C-H in plane deformation of lignin, increases. At the same time, the band at 1456 cm -1 , related to C-H bending of methoxyl groups  becomes broader. According to the literature, the IR band at 1108 cm -1 could appear because of cross-linking reactions of -OH groups of cellulose/hemicellulose with phenolic compounds at the cell wall level (Wang et al., 2016), resulting in an increase in the hardness of the wood cell walls. The most important difference is again the appearance of the band at 1648 cm -1 and the absence of the band at 1660 cm -1 in the HW tracheid cell walls. As mentioned earlier, the band at 1655 cm -1 could be mainly present because of C = O and C = C groups in coniferyl aldehyde and coniferyl alcohol structures of lignin. Hence, it is likely that the reduction of band intensity may be attributable mainly to condensation reactions of lignin molecules or oxidative alteration of lignin.
It is well known that lignin condensation makes lignocellulosic biomass more recalcitrant, mainly due to limiting the accessibility to the polysaccharides in the cell wall . Additionally, generation of new carbonyl groups by oxidation, as discussed above, may increase non-productive binding of cellulases or enzyme inhibition via chelation of metal co-factors , thereby starving the fungus. The attachment of other, potentially fungitoxic or antioxidant phenolics, i.e. flavonoids, to these newly formed reaction sites are other possibilities that lignin modification would allow for. These steric or chemical inhibitory effects are important for understanding the durability of heartwood (Valette et al., 2017).
The AFM-IR spectra of the ray also show differences between SW and HW ( Figures 6A, B).
For example, the appearance of the band around 1680-1692 cm -1 suggests the presence of conjugated ketones and carboxylic acids, such as resin acids, in the HW ray and its surrounding cell wall tracheid. The presence of taxifolin inside the ray is supported by the presence of the band at 1164 cm -1 (a characteristic band for 5,7dihydroxysubstituted flavonoids) (Zu et al., 2012).
We can also observe an intense C = O stretching vibration at 1764 cm -1 and around 1750 cm -1 (see Figure 6B), indicating the presence of g-lactones and alkyl esters, including the methyl esters of fatty acids (Lievens et al., 2011). Prominent bands that appear at 1072 and 1728 cm -1 are assigned to C-O deformation in primary and secondary alcohol groups of galactosyl-and carbonyl of carboxylic groups.
The interpretation of the spectra is complex since mixtures of chemical compounds are present throughout the plant tissue, i.e spectra of pure chemical constituents are rarely possible to obtain, neither by use of high spatial resolution (as in AFM-IR), nor by use of MCR-ALS modelling. Consequently, comparing to a spectral data base with the most common components present in plant cell walls might help to further identify the chemical compounds, but it would most likely not lead to conclusive results. AFM-IR is presented in this study as a potent technique to further characterize plant cell wall components because of its higher spatial resolution than SR-IR imaging, and because spectra could be obtained for a broader range of wavenumbers. However, measurements of a more comprehensive sample set would be necessary in order to not simply illustrate the technique but obtain representative results.

CONCLUSIONS
MCR-ALS multiset analysis on sets of SR-FTIR images collected across the HW formation zone of Kurile larch provided a cellular level description of the components involved in HW formation. In particular, the IR resolved spectral signatures and comparison with IR reference spectra from the literature allowed us to identify taxifolin, one of the most abundant extractive in larch, in rays as well as in the lumen and S3 cell wall layer of adjacent tracheids. Moreover, refolding of the concentration profiles to the original image formats allowed us to see that one initial phenolic lignan contribution (component IV) was present in the SW, while a second somewhat similar contribution (component V) emerged in the transition zone and continued in the HW. Our interpretation of this result is that component IV is a set of precursor molecules for component V. Such a pattern is characteristic for Type II heartwood formation, also called Juglans-type. The main spectroscopic difference between component IV and V was the appearance of the band at 1640 cm -1 and the simultaneous disappearance of the band at 1660 cm -1 . We hypothesize that the disappearance of the 1660 cm -1 band may be attributable mainly to condensation reactions of lignin/lignan molecules or oxidative alteration of lignin. Lignin condensation reactions are known to make lignin more recalcitrant. Generation of new carbonyl groups by oxidation of coniferyl alcohol to coniferyl aldehyde could also help explain both the peak shift and the resistance against fungal attack of Kurile larch HW.
AFM-IR has been proven to be a powerful technique to study the nanoscale compositional variations between the cell wall, CML, and the ray of SW and HW of larch. The AFM-IR results confirmed the trends observed in the SR-FTIR image analysis and provided more detail about the plant cell wall composition as spectra were obtained for a broader spectral range. Conjugated ketones and carboxylic acids accompanied with the presence of g-lactone and alkyl ester were also found in the HW rays. Finally, AFM-IR spectra proposed the existence of cross-linked reactions of cellulose/ hemicelluloses with phenol compounds at the cell wall level in HW.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
All the authors discussed the results and contributed to the final manuscript. SP, SF, RO, AG-S, TK, AJ and LT carried out the SR-FTIR experiment and SJ performed the AFM-IR measurements.

FUNDING
The research was funded by VILLUM FONDEN through project 12404. RO, AG-S and AJ also acknowledge funding from the Spanish government through project CTQ2015-66254-C2-2-P.

ACKNOWLEDGMENTS
We are grateful to ALBA synchrotron facilities for enabling the SR-FTIR measurements, especially to Imma Martínez Rovira and Ibrahem Yousef for helping in the SR-FTIR measurements. SJ wishes to thank Tue Hassenkam and the Villum foundation "Experiment" for support for AFM-IR measurements. We acknowledge Illustrations Elaborated for help with production of AFM-IR Figure.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019. 01701/full#supplementary-material FIGURE S1 | Schematic diagram of AFM-IR setup. IR laser irradiates the sample which then expands. The expansion of the sample is detected by the AFM tip which deflects the whole AFM cantilever. The deflection is then monitored by tracking the movements of the AFM laser shone on the cantilever.

FIGURE S2 | (A)
Representation of the area selection of the lumen, cell wall, and ray. (B) Average spectra of the cell wall selected for each of the cross section images collected across the heartwood formation zone. (C) Zoom of the spectral range 1200 cm -1 to 1800 cm -1 of the average spectra of the cell wall selected for each of the cross section. The spectra gradually change from blue (sapwood) to red (heartwood) color.
FIGURE S3 | MCR-ALS results of the multiset structure formed by a series of tangential section images. (A) distribution maps of components involved in the heartwood formation of Kurile larch. Each line of maps represents the resolved maps of all constituents for a particular sample. Each column of maps represents the distribution map of a particular chemical constituent in all samples analyzed. Distribution maps use a gradual color scale where yellow color refers to large concentration values and blue color to small values. (B) Related pure spectra.