Variation in Bark Allocation and Rugosity Across Seven Co-occurring Southeastern US Tree Species

Bark is a complex multifunctional structure of woody plants that varies widely among species. Thick bark is a primary trait that can protect trees from heat generated in surface fires. Outer bark on species that allocate resources to thick bark also tends to be rugose, with bark being thickest at the ridges and thinnest in the furrows. Tree diameter or wood diameter is often used as a predictor for bark thickness but little attention has been made on other factors that might affect bark development and allocation. Here we test multiple mixed effect models to evaluate additional factors (height growth rate, measure height) that correlate with bark allocation and present a method to quantify bark rugosity. We focused on seven co-occurring native tree species in the Tallahatchie Experimental Forest in north Mississippi. Approximately ten saplings of Carya tomentosa, Nyssa sylvatica, Prunus serotina, Pinus echinata, Pinus taeda, Quercus marilandica, and Quercus falcata were destructively sampled for stem analyses. Outer bark thickness (OBT) ranged from 0.01 to 0.77 cm with the thickest maximum outer bark occurring on P. taeda (0.77 cm) and the thinnest maximum outer bark occurring on P. serotina (0.17 cm). Our outer bark allocation models suggest that some individuals with rapid height growth allocate less to outer bark in C. tomentosa, N. sylvatica, P. taeda, and P. serotina, but not for P. echinata or either oak species. All species except for C. tomentosa and N. sylvatica showed evidence for outer bark taper, allocating more outer bark at the base of the bole. Inner bark also was tapered in Carya and the oaks. Bark rugosity varied among species from 0.00 (very smooth) to 0.17 (very rugose) with P. Serotina and C. tomentosa having the smoothest bark. OBT was the best fixed effect for all species. Aside from providing data for several important yet understudied species, our rugosity measures offer promise for incorporating into fluid dynamics fire behavior models.


INTRODUCTION
The development of thick bark in woody plants has widely been cited as a trait selected for by frequent fire regimes as it protects trees from the heat generated in surface fires (Spalt and Reifsnyder, 1962;Hare, 1965;Vines, 1968;Pausas, 2015). Tree species that occur in frequently burned ecosystems (e.g., savannas and woodlands) generally develop thicker bark than those found in less flammable environments (Hoffmann et al., 2003;Lawes et al., 2011;Schafer et al., 2015). This thick bark is one of a suite of traits that can enable "fire tolerant" species to persist and grow into the canopy in ecosystems with a frequent fire regime, while thin-barked "fire intolerant" species are usually top killed (Varner et al., 2016).
Bark includes all tissues from the outer surface of the bole and branches to the vascular cambium. Outer bark, or rhytidome, consists of the collection of dead cells beyond the phellogen. Outer bark appearance differs among species depending on factors such as the position of the first periderm, the development of any additional periderms, and the composition and structure of phloem cells (Borger, 1973). Outer bark is suggested to have multiple functions including insulation from heat of fire, protection from herbivory, and structural support (Vines, 1968;Rosell, 2019). Inner bark is the living portion of bark that originates from the vascular cambium and extends to the phellogen (cork cambium). Inner bark contains secondary phloem, the cortex (for some species), and the phelloderm and is responsible for the transport of photosynthates, storage of non-structural carbohydrates and water, and the closing of wounds (Romero and Bolker, 2008;Romero et al., 2009;Rosell and Olson, 2014;Rosell, 2019). Although inner bark is not generally considered a primary fire defense, its high density and water content may provide protection from fire (van Mantgem and Schwartz, 2003;Schafer et al., 2015), however, several studies suggest that moisture content in bark decreases heat resistance due to high heat conductivity (Vines, 1968;Odhiambo et al., 2014).
Although outer bark is considered the primary defense against fire, most studies on post-fire tree mortality report bark thickness based on the combined thickness of outer and inner bark, often not distinguishing between the two structures (e.g., Harmon, 1984;Hengst and Dawson, 1994;Lawes et al., 2011;and many others). In larger trees, total bark thickness is often measured at the stem base or at breast height (approximately 137 cm) with a gauge that is incapable of differentiating outer and inner barks (e.g., Harmon, 1984;Hengst and Dawson, 1994;Lawes et al., 2011). Other studies utilize cross sections of species to measure bark with calipers (Hammond et al., 2015;Shearman et al., 2018). Total bark thickness allometric equations developed using the relationship of bark thickness and stem diameter generally fit very well as linear, logarithmic, or quadratic functions (Hengst and Dawson, 1994;Jackson et al., 1999;Lawes et al., 2013;Rosell, 2016), however, the correlation is largely due to the relationship between inner bark and stem diameter as outer bark does not correlate as closely (Rosell, 2016). The poorer correlation between outer bark and stem diameter suggests that other factors may affect the allocation of resources to outer bark. Thus, studies that fail to distinguish between outer and inner bark or that utilize allometric equations only fit to total bark thickness may miss important species differences that help explain patterns in mortality studies.
A common pattern is that outer bark is rarely uniform around the bole of a tree, especially in species that develop thick outer bark at relatively young ages. These species develop rough or "rugose" bark which can be impressively thick along the "ridges" and quite thin in the "furrows" (Figure 1). Previous studies on bark thickness suggest that the thin bark of furrows to be FIGURE 1 | Turkey oak (Quercus laevis) develops thick rugose bark even in small diameter stems. more prone to injury due to the lack of thermal protection (Smith and Sutherland, 1999). However, Fahnestock andHare (1964) andO'Brien et al. (2018) demonstrated that temperature in the furrows is lower than that of the ridges during a surface fire. Whether this temperature difference between ridges and furrows can be significant enough to prevent damage to the cambium is uncertain. Despite the acknowledged variation in thickness for species with rugose bark, bark thickness is usually reported as the average thickness of multiple measurements taken around the stem. Because pyrophytic species in frequently burned woodlands and savannas have been found to have thicker bark than mesophytic species in less frequently burned forests (Odhiambo et al., 2014;Pausas, 2015;Schafer et al., 2015), it follows that they are likely to have more rugose bark as well.
Differences in the allometric relationships between bark thickness and diameter suggest that the rate of bark allocation changes throughout the life histories of some species (Jackson et al., 1999). Relative bark thickness (RBT, the proportion of bark relative to stem diameter) also changes along the height of the bole, with pyrophytic species in frequent fire regimes having a more dramatic decrease in RBT than mesophytic species (Graves et al., 2014;Hammond et al., 2015;Shearman et al., 2018). Surface fires in these ecosystems are typically characterized as having low flame heights, with trees experiencing maximum temperatures within 1 m above the ground (Fahnestock and Hare, 1964;Vines, 1968). By developing thick basal bark that tapers with height, trees that are adapted to frequent fire regimes can optimize resource allocation where it is needed most. Failing to incorporate bark taper when estimating bark thickness could result in poor models of post-fire tree mortality. For example, fire modeling programs such as FOFEM that still rely on linear models of bark thickness could lead to increased error in the estimation of bark thickness and subsequent estimations of tree survival (Zeibig-Kichas et al., 2016). Recent studies examining bark thickness as a fire adapted trait in frequent fire regimes have recognized that differences in bark allocation among species should be most apparent in saplings, as it would be necessary to quickly develop thick bark during the fire-free period to recruit into the canopy (Jackson et al., 1999;Hammond et al., 2015;Schafer et al., 2015;Shearman et al., 2018).
To identify factors, such as height growth rate or measure height, that might correlate with the development of bark, we studied bark allocation among seven species in northern Mississippi. We identify factors that correlate with inner and outer bark allocation in each species. We also present a method of characterizing bark rugosity and compare rugosity measurements among species. We hypothesized that pyrophytic species would have thicker, more rugose bark than mesophytic species. We also hypothesized that rugosity would be a function of bark thickness.

MATERIALS AND METHODS
Our study site (34.506499, −89.498432) was located within the Tallahatchie Experimental Forest near Holly Springs, Mississippi. The upland forest soils of this part of the Greater Yazoo River Watershed are acidic sandy loams and silt loams on the ridge tops with acidic loamy sands on slopes (Cannon and Brewer, 2013). Prior to fire exclusion and before extensive logging, the site was dominated by pyrophytic species such as Quercus velutina, Quercus marilandica, Quercus stellata, Quercus falcata, and Pinus echinata (Surrette et al., 2008;Cannon and Brewer, 2013). In recent years, multiple disturbance events have occurred at the site, including an EF4-intensity tornado in 2008 and three biennial prescribed burns in the 8 years following (Cannon and Brewer, 2013;Brewer, 2015Brewer, , 2016. The initial tornado disturbance reduced canopy cover to approximately 45%, which recovered to 60% by 2012 (Brewer et al., 2012;Brewer, 2015).
We selected approximately ten saplings of seven native tree species: Carya tomentosa (n = 10), Nyssa sylvatica (n = 10), P. echinata (n = 11), Pinus taeda, Prunus serotina (n = 10), Q. marilandica (n = 10), and Q. falcata (n = 11). Saplings were subjectively selected to cover the widest range of diameters possible. We cut cross-sections along the stem every 10 cm from the base (0 cm) to 100 cm and every 20 cm from 100 to 200 cm. Each cross-section was air-dried and then sanded with up to 1000-grit sandpaper. All cross-sections were aged by counting rings under a dissecting microscope. Cross-sections at 0 (or 10 cm if the basal cross-section had substantial damage from fire) and 140 cm heights (measure height) along with sections from each age break (i.e., cross-sections that had fewer rings than the previous cross-section below) were digitally scanned on a flatbed scanner at 4,800 dpi (Epson Perfection V39, Suwa, Nagano, Japan). Total cross-section area, outer and inner bark area, and wood area were measured directly using ImageJ image analysis software (version 1.53c, Schneider et al., 2012) and tracing the actual cross-sections. Wood diameter and inner and outer bark thickness (OBT) were then calculated from the area measurements based on the relationship of the diameter of a circle to its area.

Bark Allocation
For each species, we fit a series of mixed effect models using both outer and inner bark cross-sectional area as the response variable. Fixed effects included wood cross-sectional area, measure height, and height growth rate (total height/age of the lowest crosssection). A unique tree number for each individual was included as a random effect to account for the lack of independence of multiple observations on the same tree. Bark area, wood area, and height growth rate were first log-transformed to improve heteroscedasticity of the residuals. Seven separate models were fit for each dependent variable using different combinations of fixed effects. Models were first evaluated using the corrected Akaike Information Criteria (AICc, Hurvich and Tsai, 1989). Top models identified by having the lowest AICc (including models with AICc < 2, Burnham and Anderson, 2004) were further evaluated and fixed effects were checked for significance (alpha < 0.05). The model with the lowest AICc and all significant fixed effects was chosen as the best model. If more than one model was considered the best and there was no significant difference between models (based on ANOVA), the model with the fewest terms was selected. Damaged cross-sections due to fire scars were removed from the analyses prior to fitting the models (n = 4 for C. tomentosa, 5 for P. taeda, 1 for P. serotina, and 3 for Q. marilandica).

Bark Rugosity
Rugosity, f r , is a measure of the variations in amplitude of a surface and is calculated as the ratio of the actual surface area (A r ) to the geometric surface area (A g ): For our bark rugosity (B r ) measurement, we calculated the ratio of the actual cross-sectional area at each measured height, A mht , to the area of the convex hull at that measured height, A conv(mht), which we then take the compliment to scale the ratio so that increasing value corresponds to increasing rugosity (Figure 2): For each species, mixed effect models were fit using bark rugosity as a response variable and tree individual as a random effect. Fixed effects included species, wood diameter, outer and inner bark thickness, measure height of the cross-section, and height growth rate. A total of eight models were fit, varying the fixed FIGURE 2 | Bark rugosity of a cross-section, in this case a sample of Quercus marilandica at 10 cm height, can be measured as the compliment of the ratio of the cross-sectional area (area of the actual sample) to the area of its convex hull (red line). effects for each model. Model selection was made using the same methods as above.
To test significance of rugosity and OBT among different species, we used ANOVA to compare the basal cross-sections of each species. We also compared diameters of basal cross-sections among species with ANOVA to see if any differences in rugosity or bark thickness were also observed in their diameters. We then fit a mixed effect model to all the species and all cross-sections combined using bark rugosity as the response variable and OBT as a fixed effect. Again, tree individual was considered a random effect. Differences in slopes between species were identified by examining all pairwise comparisons among species.
All models in all analyses were fit in R (version 3.3.1, R Core Team, 2019) using the lmer function in the lme4 package (Bates et al., 2015). P-values for fixed effects were calculated using the lmerTest package (Kuznetsova et al., 2017). Pseudo-R 2 values were calculated with the r.squaredGLMM function in the MuMIn package (Barton, 2019). Multiple comparisons were made using the emmeans package (Lenth, 2021). Raw data are available in the Supplementary Material.
Outer bark investment was affected by different factors for different species although wood cross-sectional area was a positive significant effect for all models. Height growth rate was a negative effect (i.e., faster height growth had lower bark allocation) for C. tomentosa, N. sylvatica, P. taeda, and P. serotina, but not for P. echinata or either oak (Q. falcata, and Q. marilandica). Measure height also had a negative effect for all species except for C. tomentosa and N. sylvatica ( Table 4). All outer bark investment models had high explanatory   power, with fixed effects explaining 83-96% of the variance and combined fixed and random effects explaining 89-97% of the variance ( Table 4).
The best inner bark investment models all included wood cross-sectional area as a significant fixed effect. Measure height had a negative effect on C. tomentosa, Q. falcata, and Q. marilandica, but was not significant for the other species (Table 4). Fixed effects explained a slightly higher amount of the variance in the inner bark models than the outer bark models, explaining 93-98% of the variance. Fixed and random effects combined explained 97-99% of the variance in the data ( Table 4).
All species had the same structure for the best bark rugosity model. OBT was the only significant fixed effect in the best models. The models suggest a higher slope for Q. marilandica (0.28), N. sylvatica (0.23), and Q. falcata (0.20), than the other species (Table 5). Pairwise comparisons on the full model of all species found that Q. marilandica had significantly higher slopes than C. tomentosa, P. echinata, P. taeda, and P. serotina, but not N. sylvatica or Q. falcata. N. sylvatica, and Q. falcata also had significantly higher slopes than P. taeda but not with any other species (Table 6). Species and OBT as fixed effects in the full model explained 77% of the variance with the random effect of tree individual explaining an additional 6%.

DISCUSSION
Previous studies on bark in relation to fire have overwhelmingly found that bark thickness is a highly significant predictor in post-fire tree mortality and that bark thickness increases with stem diameter (e.g., Hare, 1965;Lawes et al., 2011;Pausas, 2015). Several studies also acknowledge the contributions of other bark properties such as density and moisture content (e.g., Pinard and Huffman, 1997;van Mantgem and Schwartz, 2003;Brando et al., 2012;Odhiambo et al., 2014;Nolan et al., 2020). Comparatively fewer studies have examined bark surface structure in relation to fire, each utilizing different methods of measurement from categorical to quantitative (Barlow et al., 2003;Bauer et al., 2010;Nolan et al., 2020). Our study provides a quantitative method that can be used to standardize bark rugosity measurements.
Allocation of resources to produce thick bark presumably comes at a cost to other functions. For example, thick bark reduces the ability to photosynthesize along the stem and branches . Our results suggest that height growth may also be a compromise to bark allocation for some species. All species in our study, with the exception of P. taeda, resprout readily after topkill, and therefore most of the individuals in our samples were likely resprouts, which could explain the significance of height growth rate in our models. In fast growing resprouts of P. serotina, for example, individuals that had the fastest height growth allocated less to outer bark than slower growing individuals. Height growth was not significant in most of the species that are considered fire tolerant (Q. marilandica, Q. falcata, and P. echinata). This pattern suggests the possibility of different priorities in resprouting between fire tolerant and fire intolerant species where the latter species prioritize height growth to compete with light availability over investment for fire defense (Hammond et al., 2015). Height growth was also significant in the outer bark models for P. taeda and for C. tomentosa, which some may consider to be pyrophytic species (Varner et al., unpublished). Pinus taeda is fire-tolerant in larger sized individuals, however, saplings can be fire-sensitive and are often killed in frequent fire regimes (Williams, 1998;Stewart et al., 2015;Robertson et al., 2021). C. tomentosa is often lumped in with other Carya species as a member of the oak-hickory forest type that was maintained under intermediate, low-intensity historical fire regimes (Abrams, 1992;Frost, 1998). However, Monk et al. (1990) suggests that the "oak-hickory" forest type is not supported by the data and should be designated as "mixed oak" or simply "oak forests." The outer bark of C. tomentosa saplings was extremely thin in our dataset despite having the thickest average inner bark. Additionally, C. tomentosa is easily topkilled by fire (Smith, 1990;Coladonato, 1992), which calls in to question C. tomentosa as a pyrophyte (but see Varner et al. in this special issue on its flammability).
Our study found an effect of measure height on OBT in the oaks and pines, which indicates that outer bark tapers more than wood in these species. This result is consistent with the findings of others that suggest that pyrophytic species develop thicker bark at the base of the bole as a fire protection strategy (Graves et al., 2014;Hammond et al., 2015;Shearman et al., 2018;Kidd and Varner, 2019). However, we also found the effect of measure height in P. serotina, which is not consistent with this hypothesis. Shearman et al. (2018) found evidence for bark taper in C. tomentosa although they only measured total bark thickness in their study. Here, we found an inner bark taper with height for C. tomentosa (as well as the oaks) which likely accounts for the discrepancies between our study and Shearman et al. (2018).
Thick bark at the base of trees, particularly fire-sensitive saplings, is intriguing but the evidence of this phenomenon as a fire-adapted trait is still lacking. Midgley and Lawes (2016) suggest that the decline in relative bark thickness with height is too small to be of ecological significance. However, studies that have identified an effect of height on bark allocation have thus far been on relatively small stems. Graves et al. (2014) sampled trees up to 12 cm DBH, Hammond et al. (2015) and Kidd and Varner (2019) samples were all under 5 cm, Shearman et al. (2018) samples were all under 10 cm DBH. Future studies need to examine this relationship in a larger range of diameters before discounting it as insignificant for fire survival. Our lack of understanding on the mechanism for this effect is also worth investigation. Finally, it is uncertain as to whether producing thick basal bark synergizes with other traits such as rugosity to better protect stems during surface fires.
Bark rugosity is largely a consequence of developing thick bark. As thick bark develops, cracking or flaking of bark must occur as the radial growth of wood exerts tangential stress on bark tissues. In pines and other species with flaky, scaly bark, bark develops as new, discontinuous periderm layers form beneath the previous periderms in overlapping arcs (Borger, 1973). These discontinuous arcs of periderm layers allow bark flakes, sheets, or platelike strips to be shed due to growth and weathering. Bark in furrowed species, like the oaks in our study, form due to interlocking sclerified tissues in the outer bark that adheres to the periderm cells resulting in bark that does not shed or flake (Borger, 1973). In our models of bark rugosity, OBT explained the most variability in the oaks and less variability in the pines, possibly due to shedding of bark flakes in the pine species. The strong relationship between rugosity and OBT in the oaks suggest that rugosity could be useful in developing more accurate bark thickness models in these species. The poor performance of our model for P. serotina was due to that species having the smoothest (i.e., least rugose) bark in the data set.
The survival benefit of rugose bark is unclear because of its correlation with OBT in the species we studied. Bauer et al. (2010) found that including bark surface structure (measured as the difference in maximum and minimum thickness divided by the mean) slightly improved their models of thermal conductivity, but stated further study was needed to understand the significance for fire resitance. Fahnestock and Hare (1964)  hypothesized that the reduction in temperature in the furrows during a fire may offset the vulnerability of having thin bark in these areas. Rugose bark could potentially alter the fluid dynamics at the fire-bark surface interface due to deeper boundary layers (Dickinson and Johnson, 2001). As a flaming front passes a tree bole, leeward vortices can develop, increasing leeward flames and causing increased stem heating on the leeward side (Gutsell and Johnson, 1996). The formation of these vortices depends on the Reynolds number, a dimensionless number based on upstream wind velocity, tree diameter, and the kinematic viscosity of the air (Gutsell and Johnson, 1996). However, this model of the fluid dynamics of stem heating assumes the tree is a smooth cylinder and the influence of surface rugosity on the fluid dynamics still needs to be explored. Butler and Dickinson (2010) justified the use of a smooth surfaced model of heat transfer through bark by reasoning that the higher temperatures observed on the ridges of a stem (Kayll, 1963;Fahnestock and Hare, 1964;O'Brien et al., 2018) would have a similar rate of heat transfer as in the fissures that experience lower temperatures. Jones et al. (2004) also reasoned that because tree mortality due to cambial necrosis would require girdling of the tree and thus damage under the thick ridges as well as the fissures, that a one-dimensional model considering heat transfer through the thick ridges alone would be sufficient in predicting mortality. However, Barlow et al. (2003) found that bark texture (measured categorically as "rough, " "medium, " or "smooth") was an important determination of tree mortality with surviving trees in burnt forests having significantly rougher bark in the smaller (0.2-0.6 cm) bark thickness classes. Barlow et al.'s (2003) finding that bark texture was only significant in trees with thinner bark might suggest that rugosity might benefit smaller diameter trees and that larger trees can withstand fire due to having thick bark alone. Therefore, future heat transfer models might benefit from including a rugosity term as proposed in our study. Although our study found significant differences in both OBT and bark rugosity among species, our results are limited to one study site and a limited range of species and sapling sizes. These results largely follow the pattern of fire-tolerant "pyrophytes" (e.g., the oaks and pines) having thicker and more rugose bark than the fire-sensitive "mesophytes" (e.g., P. serotina), however, other environmental factors could influence bark development. For example, Glitzenstein and Harcombe (1979) found differences in Q. falcata bark roughness, with drier sites (measured as percent of sand content in soil) having more furrowed bark than mesic sites. Howard (1977)  also comments that tree vigor, age, and height on the tree affect bark thickness, texture, and depth of fissures, citing Guttenberg's (1951) anecdotal evidence of low vigor resulting in rougher bark. These reports suggest a plastic response of bark thickness and rugosity, with drier, more fire-prone environments eliciting more allocation to bark within species than mesic sites. Clearly more studies are needed to evaluate this hypothesis. On a global scale, fire regime is the main environmental variable explaining bark thickness across species (Rosell, 2016;Pausas, 2017). Based on the relationship between bark thickness and rugosity, it follows that fire regime likely explains patterns in rugosity across species as well. Our study, along with the limited literature on bark structure and fire, questions whether tree survival and recruitment in frequent surface fire regimes may be more than just insulation from thick bark. Rugosity may be just as important in smaller diameter stems if it alters fire behavior around the bole.
In the absence of fire, the tradeoff of height growth and bark allocation may be one mechanism that allows fire-sensitive trees to dominate during the mesophication process (Nowacki and Abrams, 2008). Future climate scenarios project increases in both likelihood of severe droughts and in fire potential throughout large portions of the United States with projected increases in fire season of 2-3 months in the southern United States (Liu et al., 2013;Mitchell et al., 2014). These changes in climate may allow us to test hypotheses about the importance of thick rugose bark in pyrophytic landscapes through long-term studies on bark thickness and rugosity in changing ecosystems.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
TS and JV designed the study. TS collected the data and conducted the analyses, and wrote the manuscript. Both authors contributed to the article and approved the submitted version.

FUNDING
We acknowledged funding from the Joint Fire Science Program under project JFSP 13-1-04-49.