Wood Density Variations of Legume Trees in French Guiana along the Shade Tolerance Continuum: Heartwood Effects on Radial Patterns and Gradients

: Increasing or decreasing wood density (WD) from pith to bark is commonly observed in tropical tree species. The different types of WD radial variations, long been considered to depict the diversity of growth and mechanical strategies among forest guilds (heliophilic vs. shade-tolerant), were never analyzed in the light of heartwood (HW) formation. Yet, the additional mass of chemical extractives associated to HW formation increases WD and might affect both WD radial gradient (i.e., the slope of the relation between WD and radial distance) and pattern (i.e., linear or nonlinear variation). We studied 16 legumes species from French Guiana representing a wide diversity of growth strategies and positions on the shade-tolerance continuum. Using WD measurements and available HW extractives content values, we computed WD corrected by the extractive content and analyzed the effect of HW on WD radial gradients and patterns. We also related WD variations to demographic variables, such as sapling growth and mortality rates. Regardless of the position along the shade-tolerance continuum, correcting WD gradients reveals only increasing gradients. We determined three types of corrected WD patterns: (1) the upward curvilinear pattern is a speciﬁc feature of heliophilic species, whereas (2) the linear and (3) the downward curvilinear patterns are observed in both mid- and late-successional species. In addition, we found that saplings growth and mortality rates are better correlated with the corrected WD at stem center than with the uncorrected value: taking into account the effect of HW extractives on WD radial variations provides unbiased interpretation of biomass accumulation and tree mechanical strategies. Rather than a speciﬁc feature of heliophilic species, the increasing WD gradient is a shared strategy regardless of the shade tolerance habit. Finally, our study stresses to consider the occurrence of HW when using WD.


Introduction
Basic wood density (WD), which is defined as the ratio of oven-dried wood mass to fresh volume [1] is an effective biophysical property that is commonly used by wood scientists and technologists. The consideration of WD has recently been extended to other disciplines and is now widely used as a functional trait in the field of functional ecology [2][3][4][5]. The success of WD can be largely attributed to its ability to relate the carbon investment per unit volume of stem and to its integrative power of diverse properties, characteristics, traits and strategies of the materials, and organisms under study. and outer wood than in intermediate wood [34]. This pattern is difficult to reconcile with the Woodcock and Shier [29] hypothesis.
Secondary changes in wood were never considered in the analysis of WD radial gradients and patterns. Indeed, rather than originating from a particular mechanical strategy, both non-monotonic patterns of variations in WD and decreasing gradients can be explained by the formation of heartwood. A significant proportion of tropical species have nonfunctional often colored inner wood called heartwood (HW). Unlike functional sapwood (SW), HW no longer contains living parenchyma cells [36,37]. Before their death, these cells participate in the synthesis of secondary metabolites (i.e., extractives) [38] that impregnate cell walls and intercellular spaces, thus participating to tree defense against biotic stress. By their impregnation, extractives increase the mass without changing the apparent volume of wood. Thus, HW density (HWD) results from (1) the amount of material originating from xylogenesis (i.e., lignin and cellulose) and (2) the amount of extractives originating from HW formation (i.e., duraminization). Royer et al. [39] recorded a decrease of between 10% and 16% of WD after chemical extraction in Eperua falcata Aubl., a tree species with high extractive content [40]. As the presence of HW extractives is a common characteristic of intermediate to late successional species that evolved wood protection strategies [41], it is very likely that the high inner wood density observed among these species is partly related to the presence of HW.
Although HW extractives have been suggested to have an effect on the WD gradient by some authors [26,33], to date no study has specifically addressed the effect of HW on WD or assessed potential misinterpretation when inferring a species' functional strategies through radial variations in WD. Here, we report on detailed analysis of WD variations in 16 Guyanese legume species that are well distributed along the shade-tolerance continuum. Combining a fine-scale analysis of WD radial variations with both demographic and HW chemistry data the present study intends to (1) assess the specific effect of HW extractives on WD radial gradients and patterns, (2) to relate species successional status and vital rates (e.g., growth and mortality) to both WD of the raw-wood and WD free of extractives, and finally, (3) to reinterpret WD variations in terms of biomechanical and biomass accumulation strategies taking into account the effect of HW.

Study Site
All the sampling was performed in the Paracou experimental site (5 • 27 N, 52 • 92 W) in French Guiana, a region characterized by a tropical climate with two dry seasons from mid-August to mid-November and during March or April. This site is a "terra firme" rain forest belonging to the Caesalpiniaceae facies [42], a type of forest typical of French Guiana [43]. The growth of all trees with a diameter at breast height (DBH) of more than 10 cm has been fully and regularly monitored in permanent plots since 1984. Tree sampling took place in unmonitored areas of the Paracou experimental site.

Study Species
We focused on Leguminosae tree species (16 species in total) because of the following reasons. Leguminosae species are (1) a major component of Amazonian tree flora [44], (2) very abundant in French Guiana, and (3) widely distributed along the shade-tolerance continuum [42]. Thus, we expected to cover a wide range of average species-level WD as well as notable diversity in WD gradients and patterns. Second, this family also contains many species used for timber in French Guiana, meaning it is of quite high commercial interest.

Selected Trees
For this study, we only sampled trees with a DBH of between 10 and 15 cm (total height of between 8 and 25 m) ( Table 1). As the canopy position does not seem to affect WD gradients [45], sampling small trees should allow a trade-off between the time required for sampling and measurements and the observation of significant variations in WD. Moreover, trees in these DBH classes are establishing in the canopy (i.e., occupying the subcanopy strata), therefore the crown expansion has already started. Hence, we expected to find clear and significant variations in WD originating from a shift in allocation from height to diameter growth in this size range. Moreover, as interindividual variation in WD is minimal compared to ontogenetic and interspecific variations [34,41,46], we sampled only 2 or 3 individuals per species in order to perform a finer scale pith-to-bark analysis of WD (see below).

Tree Logging, Collection, and Preparation of Wood Discs
Tree DBH was measured before the trees were felled. Once the tree was on the ground, tree height was measured using a measuring tape placed along the tree. We collected one wood disc at the base of the trunk of each tree (i.e., 0.5 m from the ground or 0.8 m when butt swell was observed). All the wood discs were extracted with a chainsaw and immediately referenced and sealed in plastics bags until WD was measured in the laboratory. Within the 24 h following sampling in the field, a bench planer was used to reduce the thickness of the wood disc to 2-3 cm and obtain a plane surface to enable bark/wood and SW/HW boundaries to be distinguished precisely.

Assignment of Radial Positions and Measurement of Wood Density
One diametric (i.e., bark to bark through the pith) 2-3 cm wide wood sample, was extracted from each wood disc using a band saw. Each diametric sample was then tangentially cut into 0.5 cm segments from bark to bark, taking care to include pith in one of the segments. Bark segments were discarded. Each wood segment was numbered and assigned to a radial position (i.e., distance from the pith in cm) with the segment containing pith assigned to the central position (i.e., 0 cm). Negative positions were assigned to one side of the disc and positive positions to the other side. For each segment, we determined the nature of the tissue with the naked eye (i.e., bark, HW, SW, and pith) (Figure 1). When a segment contained both HW and SW, it was classified as a HW-SW mixed segment. Inner segments were attributed to HW when their color clearly contrasted with outer segments. The green volume of each segment was measured using the buoyancy (G) of the sample, i.e., G = weight in air (W a )−submerged weight (W fl ). The mass of the samples was then stabilized in the oven at 103 • C for 48 h and dry mass was recorded (W d ). Both measurements were made using a YDK03 density kit (Sartorius, Goettingen, Germany) and a CP224S 0.2 mg precision balance (Sartorius, Goettingen, Germany). WD was computed as the ratio of dry mass to green volume [1], i.e., W d /(W a − W fl ) or W d /G.

Comparisons of the Inner, Outer, and Overall Wood Density at Species Level
For each species, we established the inner and outer WD distributions (i.e., in WD and out WD) by pooling the 6 innermost and 6 outermost WD values of each tree (3 per radius), respectively ( Figure 1). In the present study, the term 'inner-wood' may refer to SW or HW, depending on the presence or absence of a distinct visual difference between HW and SW, whereas 'outer-wood' refers exclusively to the outermost SW. The species overall WD distribution was established by pooling all WD values. We then compared the overall WD distributions among species using a Kruskal-Wallis test followed by a post-hoc Tukey's range test. Within each species, out WD and inner inn WD distributions were also compared using a Kruskal-Wallis test. All statistical analyses in this study were performed using R statistical software [48]. We also computed the species overall WD (WD mean ) as the arithmetic mean of all WD values measured on both radii, after removing samples containing the pith. For a list of variables used in this study and their abbreviations, see Table 2. WD, HWD, SWD Wood density, heartwood density and sapwood density (g cm −3 ) WD mean , corWD mean Observed and corrected mean wood density (g cm −3 ) WD ctr , corWD ctr Observed and corrected wood density at the stem center (g cm −3 ) WD grad , corWD grad Observed and corrected wood density radial gradient (g cm −3 m −1 ) inn WD, out WD Inner and outer wood density distribution (g cm −3 ) RGR95, sRGR95 95th percentile of the radial growth rate distribution of mature tree and sapling (cm year −1 ) sMR Sapling mortality rate (Proportion of dead saplings per year)

Decoupling the Effect of Xylogenesis from that of Heartwood Formation on WD
In order to separate the effects of wood formation and HW extractives production on WD, we estimated the WD produced by young trees (i.e., HW-free inner wood), by (1) carrying out a WD correction based on preexisting or newly acquired chemical data as well as (2) the collection and measurement of WD in young individuals of a subset of 3 species.

Estimation of Wood Density Corrected by Extractive Content
We computed WD corrected by extractive content (corWD) of HW samples for species showing HW. As the deposition of extractives within the cell wall during HW formation increases the wood mass, corWD was computed as the ratio of corrected dry mass (cW d ) to green volume (G). cW d was computed using extractive content (extractive mass per wood dry mass) at the species level originating from several kinds of sources (Table S1). The main source was the CIRAD (French agricultural research and international cooperation organization) wood chemistry database (CIRAD, unpublished), which contains 35,823 entries from 671 taxa (577 identified to the species level). We also used a few published studies that addressed specific questions regarding the HW natural durability and extractive content in the species selected for the present study (Table S1). As extractive content data were unavailable for Zygia tetragona Barneby & J.W. Grimes, we proceeded to conduct chemical extraction in water and methanol from HW powder using the same methodology as the CIRAD wood chemistry database (methodology of the Centre Technique Forestier Tropical (CTFT)) (Appendix A, for details on methods and results) to measure extractives content and used this value to compute corWD for this species. We also computed the species overall corrected WD (corWD mean ) as the arithmetic mean by pooling all corWD values from HW and all WD values from SW values, after removing samples corresponding to bark and pith.

Selected Young Trees of Three Species
Another possible way to analyze the effect of HW extractives on WD variations would be to compare the WD radial profiles of individuals in which HW formation has not yet started with trees where HW is present and clearly demarcated. However, the timing or size that triggers the HW formation process depends on the species and/or the growth environment [36,49], so there is no general rule to enable the selection of trees of sufficient size (i.e., showing significant variations in WD) that do not contain HW. Nonetheless, we decided to compare the sapwood density (SWD) of the trees described above with SWD collected from younger/smaller trees that are likely to not contain HW, or at the most, only very small amounts.
Three species-Bocoa prouacensis Aubl., Eperua grandiflora (Aubl.) Benth., and Zygia tetragona-were compared using 3 to 4 individuals per species whose DBH ranged from 2.5 to 7.5. As explained above, we measured the WD radial profiles of these smaller trees at 0.5 m in order to compare them with radial patterns of bigger trees measured at the same height.

Successional Status and Demographic Variables
Successional status and demographic data for each species sampled were gathered from different kinds of sources. For the successional status of each species, we used the Favrichon [47] classification, established on the Paracou experimental site, which defined 3 classes of species along the shade-tolerance continuum: heliophilic, semi-tolerant, and shade-tolerant species (Table 1). Successional groups were determined by Favrichon [47] through mutltivariate analysis and classifications with mean diameter, proportion of recruits (tree with DBH >10 cm), and annual growth rate for both natural and perturbed forest at the species level as input variables. In order to link the different WD variables to demographic variables of both mature trees and saplings, we used data available from different census protocols. For mature trees, we used data from the GUYAFOR network database (Guyafor, DataBase of the French Guiana Permanent Plot Network, CIRAD-CNRS-ONF) in order to compute the maximal radial growth rate (in cm/year) and the maximum DBH for each species (DBHmax). The maximal radial growth rate was computed as the 95th percentile of the radial growth rate distribution (RGR95) containing the RGR of each tree averaged over a 30-year period. Data were restricted to the Paracou experimental site in order to avoid potential biases due to differences between sites. For juvenile trees, as the GUYAFOR census protocol considers only individual trees with a DBH greater than 10 cm, the variables for the saplings were computed using two data inventories provided by Molino and Sabatier [50] and Vincent et al. [51]. Between 1995 and 1997, Molino and Sabatier [50] conducted a census, measured and identified all trees with a DBH ranging from 2 to 10 cm, over ten 20 × 250 m transects at Paracou. Vincent et al. [51] recorded mortality events along the same transects, and measured the DBH of the saplings that remained in 2002. From these two measurements, we computed the sapling maximal radial growth rate (sRGR95) at the species level in the same way as for the mature trees of the GUYAFOR network. Molino and Sabatier [50] and Vincent et al. [51] data inventories, were also used to compute the sapling mortality rate (sMR) as the proportion of the number of dead saplings divided by the mean period between the two censuses. To avoid bias due to small sample size, species for which RGR95, sRGR95, and sMR was computed using fewer than 20 individuals were discarded. As a result, among the 16 species used for WD measurements; the sRGR95 and sMR values of 6 species were discarded and the RGR95 values of 2 species were discarded. For a list of variables used in this study and their abbreviations, see Table 2.

Modeling the Radial Pattern
The WD radial variations of the SW and HW in the sampled individuals were modeled at the species-level using multilevel linear mixed effect models [52]. These models enable integration of random factor effects that may significantly impact estimation of the model parameters and its significance. As variations in WD along the pith-to-bark radius are a common feature of tropical trees, we used distance from the pith as a fixed effect and selected 2 nested random factor levels to allow radial variations in WD to be modeled at different scales in the same species. These 2 nested random factors levels were (1) the individual and (2) the measured pith-to-bark radius within an individual.
As curvilinear patterns exist [20,34], we extended our modeling approach to include quadratic expressions.
In this way, we specified both linear and quadratic full multilevel mixed-effects models, for whom all terms have a random effect at all nested levels considered. Let WD ijk be the WD value of the kth 0.5 cm interval from the pith of the jth pith-to-bark radius within the ith individual, our full species-level models, in the linear (Equation (1)) and quadratic (Equation (2)) forms, respectively, are expressed as follows.
with x corresponding to the distance from the pith and ε ijk corresponding to the within group error. We also specified two mean models (i.e., no radial distance effect), that include (Equation (3)) or do not include (Equation (4)) an individual random effect, in order to highlights species with no significant radial gradient:

Model Selection Procedure
To evaluate the most suitable model, both quadratic (Equation (1)) and linear models (Equation (2)) were derived in 2 sets of reduced models varying by the inclusion of the different random error terms (7 models per form including the full model of each form). Next, all these models and the two mean models (Equations (3) and (4)) were fitted to the species-level data. All the candidate models were compared with the corrected Akaike Information Criterion (AICc) [53], available in the 'AICcmodavg' R package (2.1-1) [54], penalizing models with more parameters in order to keep the most parsimonious and suitable model for each species. We then computed differences between the lowest AICc and the AICc of each model. We only considered fits with the lowest AICc and fits for which AICc differences did not exceed 2, as candidate models. For each candidate model, we checked (1) the validity of linear-mixed model assumptions by graphical examination of residuals [52] and (2) the significance of the fixed effects. Finally, we selected the model with the lowest parameters numbers that did not violate model assumptions and with significant fixed effects. Log-likelihood ratio tests were performed to compare nested models (i.e., models that differ only by the inclusion of one or several random effects) and to check the significance of the inclusion of random effects.

Determination of the Species Corrected WD Radial Pattern
The corrected WD radial pattern of each species was determined following the modeling procedure. The determination of the corrected WD radial pattern is straightforward for species without HW (i.e., one single fit per radius), but is not for HW-species for which we modeled separately HWD and SWD. Thus, for HW-species we fitted simple linear and quadratic regressions over WD data compounded by SWD values and corrected HWD values and applied the same procedure as explained above. Three types of qualitative patterns were determined based on the mathematical form of the selected model: "linear" for linear form, "upward curvilinear" for quadratic form with β 1 < 0 and β 2 > 0, and "downward curvilinear" for quadratic form with β 1 > 0 and β 2 < 0.

Calculation of Observed and Corrected Radial Gradients
Based on the modeling of radial patterns, we computed the observed WD radial gradient of each species (i.e., the change in WD per radius unit) as where WD ctr equals the WD at the stem center (i.e., the estimation of β 0 from the whole stem or HW modeling of HW-free and HW species, respectively) and WD x refers to the predicted WD at the radial distance x in m. Here, we chose x = 0.06 m because this was the outermost value covered by all the species sampled in this study. The corrected radial gradient of each species was computed as with corWD ctr equal to the corrected WD at stem center. For species for which both mature trees and saplings were collected, we used the intercept estimated from sapling modeling as corWD ctr . When no saplings were collected, we used the intercept estimated from SWD modeling after checking the consistency between this parameter and the observed range of WD at the stem center corrected by extractive content. The corrected WD gradient (corWD grad ) was not computed for species whose extractive contents were not available.

Correlations among WD Variables and between WD Variables and Demographic Variables
We then proceeded to correlations among WD variables (i.e., WD mean , corWD mean , WD ctr , corWD ctr , WD grad , and corWD grad ) and between WD variables and demographic variables (i.e., RGR95, sRGR95, sMR, and DBHmax). After checking the normality of these variables, we used the Pearson correlation with the rcorr() function implemented in the 'Hmisc' package (4.1-1) [55]. The significance of each correlation coefficient was assessed by the p-value provided by the rcorr() function.

Linking Successional Status and Corrected WD Radial Pattern to WD and Demographic Variables
Significant differences between both WD and demographic variables and successional status and corrected WD radial pattern were assessed using the Kruskal-Wallis rank test followed by a post-hoc Tukey honest significant differences test when significant differences were revealed by the Kruskal-Wallis test. Kruskal-Wallis and Tukey tests were performed using 'base' (3.6.0) [48] and 'agricolae' (1.3-0) R packages [56], respectively.

Observed Overall, Inner and Outer Wood Density at the Species Level
Our sampling covered a wide range of WD values at both sample level (0.186 in P. nitida Miq. to 1.156 g cm −3 in B. prouacensis) and species level (0.304 in P. nitida to 1.055 g cm −3 in Z. tetragona). The observed range of WD was continuously covered by the 16 selected species (Figure 2). Species with low overall WD exhibited lower inn WD than out WD, and inversely, species with high overall WD exhibits higher inn WD than out WD. Interestingly, species in which inn WD > out WD all exhibited clearly demarcated HW, apart from V. Americana Aubl., for which no significant differences between inn WD and out WD were observed. However, species in which inn WD < out WD are both HW-free and HW species. T. guianensis (Benth.) Zarucchi & Herend. was the only HW-free species that exhibited no differences between inn WD and out WD (Figure 2).

Qualitative Assessment of WD Radial Patterns in Light of the Presence of Heartwood
Taking the presence of HW into account, the observation of variations in WD from the pith to the bark enabled three qualitative types of radial patterns to be distinguished (Figure 3) resulting in the above-mentioned differences between inn WD and out WD (Figure 2). Type 1 was represented by species with no HW in which WD increased from the pith to the bark (e.g., P. velutina Benoist, P. pendula (Willd.) Walp., P. nitida, I. stipularis DC., and E. schomburgkii (Benth.) Benth.). Type 2 corresponded to species with HW in which WD increased from the pith to the bark (e.g., R. speciosum (Benoist) Gazel ex Barneby and A. jupunba (Willd.) Britton & Killip). Type 3 corresponded to species with HW in which the HWD was higher than SWD (E. falcata, E.grandiflora, B. prouacensis, and Z. tetragona). Taking into account the presence of HW and SW highlights the effect of HW on WD radial changes ( Figure 3). Indeed, among Type 3 species, a sharp decrease in WD occurred at the HW-SW transition. The WD values of mixed type samples (i.e., samples made of SW and HW) were often intermediates to samples containing either only SW or only HW. In Type 2 species, the WD of mixed type samples was lower than that of the adjacent HW samples. Figure 3. Diversity of radial patterns of WD variations and the effect of heartwood. Wood density is plotted according to distance from the pith for 9 species sampled. Each column refers to a qualitative type of radial pattern (Type 1, HW-free species and increasing WD; Type 2, HW species and increasing WD; and Type 3, HW species and decreasing WD). Pith, sapwood, sapwood-heartwood, and heartwood samples are indicated with black, yellow, orange, and brown symbols, respectively. The vertical dotted lines indicate the center of each stem (i.e., the pith). Pictures displayed below each plot show the diametric sample considered and the absence/presence of heartwood.

Modeling of WD Radial Variation in Heartwood and Sapwood
Modeling the radial variations of WD reveals that HW-free species exhibit exclusively positive radial patterns ( Figure 4A-D,F,G,M). Three out of the seven species have an upward curvilinear pattern (i.e., the slope increases with increasing distance from the pith; Figure 4A-C) whereas three others have linear patterns (Figure 4D,F,G) and one has a downward curvilinear pattern (i.e., the slope decreases with increasing distance from the pith; Figure 4M). Among the species with HW, A. jupunba, D. guianensis Amshoff and R. speciosum only show a slight decrease in WD at the HW/SW boundary ( Figure 4E,H,L) resulting in an almost continuous downward curvilinear pattern. In these three species, HWD increases significantly, whereas SWD only varies significantly in A. jupunba and D. guianensis ( Figure 4E,L, respectively). In the others species with HW, a marked decrease in WD can be seen at the HW/SW boundary ( Figure 4I-K,N-P). In V. americana, E. falcata, and E. grandiflora ( Figure 4I-K), both SWD and HWD increase linearly. In B. prouacensis and Z. tetragona ( Figure 4O,P), there is a significant downward curvilinear pattern in HW, whereas SWD increases linearly in B. prouacensis and is constant in Z. tetragona. P. venosa (M.Vahl)Benth. is the only one species with a significant pattern inversion between SW and HW with decreasing and increasing HWD and SWD, respectively ( Figure 4N). The range of corrected inner WD and the estimation of WDctr by modeling SWD variations of bigger individuals (i.e., the value of the intercept β 0 , see Table S2) are in partial agreement for V. americana, P. venosa, E. falcata, and E. grandiflora, the estimation of β 0 is in the range of corrected inner WD, whereas for B. prouacensis, Z. tetragona, and D. guianensis, estimated β 0 is below the observed range (the description of the selected models is available in Table S2. The statistics of all the fitted models are given in Table S3).  Figure 2).

Effect of Extractives in Bocoa prouacensis, Zygia tetragona, and Eperua falcata
Interestingly, the juxtaposition of the SWD modeled patterns of saplings and larger trees produces remarkably continuous tendencies across the whole radius ( Figure 4N-P). Modeling of SWD variations in young trees reveals only positive linear patterns of variation in the three species considered. A new comparison between the ranges of corrected inn WD and the estimation of WD ctr by modeling SWD variations on saplings (i.e., the value of the intercept β 0 , see Table S4) shows that the estimation of β 0 lies far below the corrected inn WD range for these three species (Figure 4N-P).

Observed and Corrected Radial Gradient
The observed radial gradient (WD grad ) varies from −2.74 to 5.84 g cm −3 m −1 in P. venosa and P. nitida, respectively, with an average value of 1.82 g cm −3 m −1 ( Table S3). The computation of the corrected radial gradient (corWD grad ) only resulted in positive gradients with the lowest expressed by P. venosa (i.e., corWD grad = 0.75 g cm −3 m −1 ) and increase the average value to 2.75 g cm −3 m −1 .

Correlations among WD Variables and between WD Variables and Demographic Variables
Significant correlations among WD variables are observed (Table 3). Whereas uncorrected and corrected mean WD and WD at the stem center (i.e., WD mean , corWD mean , WD ctr , and corWD ctr ) are highly positively correlated (r > 0.95), these three variables are negatively correlated with the uncorrected WD gradient (WD grad ) suggesting that the higher the averaged WD variables, the lower the value of WD grad . However, although still significant, the correlations between corWD grad and other WD variables are weaker.
Regarding demographic variables, WD variables are significantly correlated with RGR95 and sRGR95. The higher the values of WD mean , corWD mean , WD ctr , and cWD ctr , the lower the values of both RGR95 and sRGR95, although negative correlations between sRGR95 with WD mean , corWD mean , WD grad , and cWD grad are marginally significant. However, no significant correlation is apparent between DBHmax and WD variables. Regarding sMR only a significant negative correlation with cWD ctr is apparent, while interestingly, sMR is not significantly correlated with WD ctr . Table 3. Pearson correlations among WD variables and between WD variables and demographic parameters. Stars indicate the level of significance for each Spearman coefficient (*** <0.001); WDmean and corWDmean, uncorrected and corrected mean WD; WDctr and corWDctr, uncorrected and corrected mean WD at the stem center; WDgrad and corWDgrad, uncorrected and corrected WD gradient; RGR95, 95th percentile of mature tree radial growth rate; sRGR95, 95th percentile of saplings radial growth rate; sMR, sapling mortality rate; DBH, maximal DBH value.

Links between WD Variables, WD Profile, and Successional Status
Both observed and corrected species-average WD, as well as WDctr, differed significantly according to the shade tolerance index ( Figure 5, upper panel). Values for heliophilic species were significantly lower than those of semi-tolerant and shade-tolerant species for these variables. Despite the systematically higher average values observed for shade-tolerant species (Table S5), we did not detect any significant differences between shade-tolerant species and semi-tolerant species ( Figure 5, upper panel). Neither WD grad nor corWD grad showed any significant differences between shade tolerance indexes ( Figure 5, upper panel).
The same tendencies emerged when distinguishing the different types of radial pattern using the WD variables ( Figure 5, lower panel). For this comparison, we classified our species in three groups of corrected radial pattern by visual examination of the radial tendencies ( Figure 4): (1) upward curvilinear pattern (i.e., P. nitida, P. velutina Benoist, and P. pendula); (2) linear (i.e., T. guianensis, E. schomburgkii (Benth.)Benth., I. stipularis, D. guianensis, E. falcata, E. grandiflora, V. americana, and P. venosa); and (3) downward curvilinear pattern (i.e., A. jupumba, R. speciosum, S. panacoco (Aubl.) Cowan, B. prouacensis, and Z. tetragona). Species with an upward curvilinear pattern showed the lowest value of both WDmean and corWDmean, as well as WDctr and corWDctr whereas species with a downward curvilinear pattern showed the highest. Again, despite systematically higher averages values observed for species with a curvilinear downward pattern (Table S5) we did not observe any significant differences between species with a curvilinear pattern and species with a linear pattern ( Figure 5, lower panel).
The same analysis performed with demographic variables instead of variables related to WD, led to significant differences in RGR95 and sMR between shade-tolerance indexes (Table S5). Heliophilic species had the highest RGR95 and sMR, whereas shade tolerant species had lower values for these variables; semi-tolerant species were between the two extremes. More surprisingly, we did not detect any significant differences between shade tolerance indexes according to sRGR95 or DBHmax. When conducting the same analysis according to WD radial patterns, we did not find any significant differences (Table S5).

Discussion
Through the combination of fine-scale analysis of WD variation, demographic data, and HW extractive content value, our study (1) underlines the risks of functional misinterpretation of WD radial variation when WD is not corrected by HW extractives content value, (2) proposes a classification of corrected WD along the shade-tolerance continuum that enables avoiding this bias, and (3) highlights the functional significance of corrected WD value at the stem center.

The Effect of HW Extractives Impregnation on Wood Density Radial Gradients and Possible Functional Misinterpretation Due to Uncorrected Wood Density
We differentiate three qualitative types of WD patterns regarding the occurrence of HW and the differences between inn WD and out WD (Figure 3). Type 1 is HW-free and exhibits higher out WD than inn WD. However this trend is also observed in Type 2 species with HW. While the non-consideration of HW in these two types will lead to the same qualitative evaluation of WD variation (i.e., greater under the bark than at the center of the stem), from a quantitative point of view, the radial gradient in Type 2 species will be biased by the production of HW and underestimated. However, since the HW formation process increases the mass of the wood through extractives synthesis, the failure to take the presence of HW into account only affects the magnitude of the gradient but not its direction in that case. This statement does not hold for Type 3 species with both HW and a greater WD at the stem center than under the bark. This result confirms the effect of HW on WD radial variations already suspected by several authors [26,33]. Here, prior knowledge of WD value before HW formation is needed or can be estimated trough WD correction or by sampling young individuals in which HW formation has not started yet. Such information is indispensable to differentiate both the direction and the magnitude of the radial gradient attributable to xylogenesis alone from the gradient attributable to both xylogenesis and HW formation and to accurately relate WD radial pattern to tree growth strategy.

On the Corrected Radial Gradient, Patterns and Their Link with Overall, Inner/Outer WD, and the Shade-Tolerance Continuum
After correction of the WD using the concentration of extractives, our analysis revealed that all the species we sampled exhibited a positive WD radial gradient ranging from very strong (e.g., Parkia genus representatives) to very slight (e.g., B. prouacensis and Z. tetragona) ( Figure 5). These findings contrast with the bulk of studies on WD radial variations [26,29,[32][33][34][35], which reported both negative and positive gradients. We presume that this difference is mainly due to the failure to take into account the presence of HW that was previously hypothesized by several authors [26,33]. Another important effect could be related to the size and the ontogenetic stage reached by the individuals sampled, this aspect will be discussed in the last section of the discussion.
Our modeling of the corrected WD variation according to distance from the pith, revealed three different types of radial patterns: (1) an upward curvilinear pattern in species with the lowest overall WD as well as the lowest WD at the stem center (WD ctr ); (2) a linear pattern in species with an intermediate to high WD (according to observed and corrected overall WD and WD ctr ), but that occurs in a wide range of WD gradients; and (3) a downward curvilinear pattern mainly observed in species with very dense wood and in those with the least radial variations.
Despite the nonsignificant differences in corrected WD radial gradients between the different types of radial patterns (Figure 5), the highest WD gradients were observed among species with the upward curvilinear pattern. This recently described pattern of WD variation [20,34,57], in fact appears to be relatively common.
The striking agreement between the distributions of the radial pattern on one hand and the index of tolerance on the other hand, according to WD variables ( Figure 5), led us to characterize the relationships between these two criteria. In the 16 species in which we clearly identified the type of radial pattern: (1) the upward curvilinear pattern (i.e., three species) was only represented by heliophilic species; (2) the linear pattern was mostly represented by species with early to intermediate successional characteristics (i.e., two heliophilic and five semi-tolerant), V. americana being the only shade-tolerant species with a linear pattern; (3) the downward curvilinear pattern was represented by species with intermediate to late successional characteristics (i.e., two semi-tolerant and three shade-tolerant species) (Figures 4 and 5). Therefore, the gradual change in the corrected WD radial patterns, i.e., from upward to downward curvilinear via a linear pattern illustrates a gradual change in species-successional characteristics. Such findings provide a coherent framework for the inference of biomass allocation and mechanical strategies of establishing trees along the shade tolerance continuum ( Figure 6). Dashed lines indicate the observed WD range, gray band and solid line indicate the corrected WD range and the median trajectory, respectively. Note that for heliophilic species, the corrected and observed range are equal as the heliophilic species sampled are HW-free.

Do Low WD Gradients Only Illustrate Late-Successional Species?
Although we did not find any negative WD radial gradients, our study corroborates a common scheme of decreasing WD radial gradient from early to late-successional species [25,29]. Here, the following elements could explain this observation. Saplings of early-successional species exhibit high height growth rate [57][58][59]. This high investment in height growth should be achieved by the production of low WD at stem center. Thus, differential carbon allocation strategy between height growth and stem thickening represents a competitive advantage during the early stages. As height growth rate decreases from early to late successional species [59], we expect that such a difference in carbon investment between height growth and stem thickening during the early stages should be less pronounced in mid and late-successional species. This hypothesis agreed with the decreasing WD radial gradient and the increasing WD at stem center observed from early to late-successional species. However, it is important to note some deviations from this common scheme. Indeed, not only late-successional species have an almost constant WD from pith to bark. In our sample, T. guianensis, exhibited a significant but quite low WD gradient (i.e., 0.9 g cm −3 .m −1 ). Therefore we can speculate that the evolutionary history of this species favored a strategy that made it possible to achieve sufficient mechanical strength, at least by young individuals, with the same carbon investment, through the production of a larger but less dense stem [60,61]. However, this alternative strategy to the common pith-to-bark increase in WD, would be expected to increase the maintenance cost associated with the higher amount of bark and wood tissues [61]. The reduction in maintenance cost achieved by reducing the stem area is probably the main evolutionary driver of the common occurrence of a WD gradient. Reciprocally, given the negative correlation between WD and water storage [8,62], very low and almost invariable WD can be seen as an adaptive strategy of succulent-stem species, such as the African Baobab (Adansonia spp., Malvaceae) [63], at the cost of higher maintenance requirements. Interestingly, other representatives of the family Malvaceae in the Neotropics contain almost constant and very low density wood [24,25]. Given the above-mentioned functional trade-off between WD and water storage ability, it is important to emphasize the possible effect of wood multifunctionality as a driver of WD gradient magnitude.

Limitations to the Correction of WD by Extractive Content
The use of species-level extractive content to compute the corrected WD proved to be very useful to describe and understand WD variation exclusively produced by xylogenesis. Nevertheless, the partial agreement between the range of corrected inn WD and the estimation of WD ctr by modeling SWD variations revealed some limitations concerning this method.
First, from a mathematical point of view, it is important to mention that the parameters extrapolated from the polynomial model are not based on the true behavior of the data, thus more confidence should be placed in the interpretation of corrected WD ctr estimated from linear modeling than from polynomial modeling. For instance, polynomial modeling of SWD variation in D. guianensis-the most exploited timber species in French Guiana-produced a too low estimate of WD ctr , i.e., 0.12 g cm −3 , to be considered as a plausible value. A linear model fitted on the whole D. guianensis SW data provided a plausible WD ctr value of 0.645 g cm −3 . Except for this consideration, the disagreement between the ranges of corrected inn WD and extrapolated corrected WD ctr might result from intraspecific variability, which might not be sufficiently covered given the insufficient number of individuals sampled in the CIRAD wood chemistry database (Table S1). Moreover, HW extractives are known to vary radially within the same individual [64][65][66][67][68], which is supported by the observed radial variation of wood natural durability within the same individual [69], but might also vary positively or negatively from pith to bark depending on the species considered [64,70] or on the type of extractive [65][66][67]. Given the agreement between the range of corrected inn WD and the estimation of corWD ctr observed in P. venosa, we hypothesize that the decrease in the HWD from the pith to the HW/SW boundary might be due to the decrease in total extractive content. Finally, SW also contains extra material such as nonstructural carbohydrates, gums, resins, silicate, or precursors of extractives [37,49,71,72] that may impact SWD. Although it is generally acknowledged that the concentration of extra materials is generally lower in SW than in HW [70], some species may have high extractive content in the SW [36,72]. In this sense SWD should also be corrected by extra-material content.

WD Variations and Demographic Variables
Our results confirm the negative relationship between the radial growth rate and averaged WD of both saplings and mature trees [2,14,73]. Correlations involving sRGR95 and averaged WD are weaker and marginally significant compared with correlations involving RGR95 except for both corrected and uncorrected WD at the stem center (corWD ctr and WD ctr ). Interestingly, we found a better correlation between sRGR95 and corWD ctr than between sRGR95 and WD ctr . As it seems very likely that the HW formation process has not yet started in saplings, this result could be explained by the fact that it is more relevant to compare the growth rate of an individual and the wood it produces at the same time. It is generally acknowledged that high WD density at the stem center enables saplings to resist falling objects thanks to their high mechanical strength [3,6] and high natural durability [4,30]. Thus, in the same line of thought as above, it should be more relevant to link the sapling mortality rate (sMR) to the corWD ctr . Again, both the significant correlation between sMR and corWD ctr and the nonsignificant correlation between sMR and WD ctr stress the relevance of corWD ctr in predicting sapling mortality rate. This result agrees with the recent study of Osazuwa-Peters et al. [46], reporting a nonsignificant correlation between WD measured close to the pith on mature trees (i.e., uncorrected WD) and sapling mortality rate in 20 species from Barro Colorado Island, Panama. Our results show that the interspecific differences in corWD ctr are significantly explained by the observed differences in sMR. Thus, the mechanical characteristics required to survive in the understory are achieved by a relatively dense SW. The stiffness of the wood material is mainly provided by the cellulose produced during the cell wall maturation process [74], in this sense the deposition of secondary metabolites during HW formation should not affect the stiffness of the material, even though the opposite results obtained by some authors make the effect of extractives on mechanical properties a controversial topic [75]. Moreover, an increase of mass around the neutral axis of the stem is not of particular interest as bending stresses are exerted on the outer wood. Cellulose microfibril angle and fiber organization are also correlated with the stiffness of wood [76][77][78], underlining that WD should not be considered as the unique driver of wood stiffness.
Since HW is absent or poorly developed in sapling stems, we would also expect that the protective function of the HW is not effective but can be provided by extractives (or the secondary metabolites) produced by parenchyma cells and stored in the SW [49,79,80]. It is noteworthy that HW natural durability is mainly afforded by extractives content rather than WD alone [81] and that there is not any mechanistic evidence of the role of WD in pest resistance. The positive relationship between pest resistance and WD [4,30] might result from the occurrence of high extractives content in species with dense wood.
These considerations highlight the strong potential of the corrected WD value at the stem center (corWD ctr ) for the inference of sapling functional strategies in a retrospective context. The corWD ctr should be used as an unbiased proxy of both radial and mortality rates of young individuals instead of the uncorrected WD ctr . However, this statement only holds for species containing HW, i.e., semi-tolerant to shade-tolerant species. The simple uncorrected WD ctr can be used for HW-free early successional species. A significant positive correlation was also observed between uncorrected WD gradient (WD grad ) and RGR95, confirming the results of Hietz et al. [33]. Moreover, the nonsignificant differences in the corWD grad between shade-tolerance indexes confirm that a wide range of corWD grad can be observed in the same successional class [32] and that corWD ctr is a better proxy of the species position along the shade tolerance continuum [33].

What Can We Expect for Wood Density Gradients of Canopy Trees?
Whereas our study reveals that establishing trees present only increasing or constant WD values from the pith to the bark, it is legitimate to question if the same WD gradients and patterns would be valid for mature canopy trees. In the study of WD radial variation on 304 species in Panama and Ecuador, Hietz et al. [33] observed WD radial gradients ranging from −2 to +2.5 g cm −3 m −1 . Regarding the higher specific diversity sampled by Hietz et al. [33], we would expect a smaller range of WD gradients in our study. However, our range of WD grad is larger, i.e., −2.74 to 5.84 g cm −3 m −1 . We assume that this difference originates from the interaction between the size of the sampled tree and the widespread occurrence of a nonlinear WD radial pattern. For species with a downward curvilinear increase, the WD tends to asymptote with an increase in tree diameter. Thus, WD grad measured on bigger individuals will be lower than WD grad measured on smaller trees. Nevertheless, this cannot hold for species with an upward curvilinear or linear pattern. For instance, the WD radial gradient of P. nitida and S. panacoco, i.e., 5.84 and 2.78 g cm −3 m −1 , respectively, would produce WD values under the bark of 4.11 and 1.55 g cm −3 , respectively, for trees measuring 40 cm DBH, which are beyond the range of possible WD values given that the cell-wall density is approximately 1.5 g cm −3 [82]. This highlights the limit of our sampling in the inference of WD patterns for bigger individuals. We might, however, expect that in P. nitida and S. panacoco, the WD value would level off in bigger trees, in turn shifting the radial pattern observed in smaller individuals, as observed in Bagassa guianensis [68]. Thus, bigger individuals of S. panacoco should exhibit a downward curvilinear increase, whereas bigger P. nitida trees should exhibit a sigmoidal radial pattern, i.e., the combination of an upward and downward curvilinear increase not described in our study. This last consideration and the absence of any significant correlation between the maximum diameter of the species (DBHmax) and WD variables of establishing trees, highlight that these variables cannot be linked to, or be predictive of, the maximum size of the species, due to a size-effect already mentioned by Plourde et al. [35].

Conclusions
Whereas it was previously considered that late-successional species produce relatively denser wood at the stem center than at the periphery, in order to meet mechanical requirements and to resist pests, our study on the effect of heartwood formation on WD radial gradients in Amazonian leguminous tree species showed an exclusively significant positive pith-to-bark WD gradient.
This finding underlines the importance of considering the effects of extractives on WD to avoid biased interpretations when inferring species' mechanical strategies. We also point out that even if the early successional species have the strongest positive WD gradients, semi-tolerant species may also exhibit considerable positive variations. Thus, the particular allocation strategy enabling fast height growth followed by the mechanical reinforcement of the structure through the increasing WD from pith to bark is a strategy shared by all species in our sampling, regardless of their shade tolerance characteristics or position along the growth-mortality trade-off, rather than a specific feature of heliophilic species. From a practical point of view, the way that corrected WD varies radially, i.e., the corrected WD radial pattern, could be a relevant indicator of species shade tolerance indexes as a complement to both mean WD and WD gradient. Finally, our study stresses the importance of acknowledging and considering the effect of heartwood when inferring mechanical or other wood-related functions, especially for saplings, when using WD as a proxy in a retrospective approach. A step forward in the understanding of WD variations could be achieved through the integration of both WD and extractives content variation in trees.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/10/2/80/s1, Table S1: Extractive content values used to compute wood density corrected by the mass of extractives, Table S2: Description of the selected models, Table S3: Description of the fitted model. For each species the table provides the coefficients estimate, Table S4: Comparisons of wood density at stem center (WD ctr ) values obtained with the different approaches, Table S5: ANOVA results between shade-tolerance indexes or corrected WD radial pattern and wood density variables and demographic variables. A 1 g (M 0 ) sample of powder was mixed in 25 mL of methanol-water solution (80:20 v/v) for 72 h. The solution was then filtered with a filter paper and concentrated using a Rotavapor. The dry weight of the sample (M 2 ) allowed us to determine the extractive yield (YieldE) as YieldE = M 2 *100/M 0' , Three replicates were performed per individual and we used the mean value (See Table S1).