Climatic responses and variability in bark anatomical traits of 23 Picea species

In woody plants, bark is an important protective tissue which can participate in photosynthesis, manage water loss, and transport assimilates. Studying the bark anatomical traits can provide insight into plant environmental adaptation strategies. However, a systematic understanding of the variability in bark anatomical traits and their drivers is lacking in woody plants. In this study, the bark anatomical traits of 23 Picea species were determined in a common garden experiment. We analyzed interspecific differences and interpreted the patterns in bark anatomical traits in relation to phylogenetic relationships and climatic factors of each species according to its global distribution. The results showed that there were interspecific differences in bark anatomical traits of Picea species. Phloem thickness was positively correlated with parenchyma cell size, possibly related to the roles of parenchyma cells in the radial transport of assimilates. Sieve cell size was negatively correlated with the radial diameter of resin ducts, and differences in sieve cells were possibly related to the formation and expansion of resin ducts. There were no significant phylogenetic signals for any bark anatomical trait, except the tangential diameter of resin ducts. Phloem thickness and parenchyma cell size were affected by temperature-related factors of their native range, while sieve cell size was influenced by precipitation-related factors. Bark anatomical traits were not significantly different under wet and dry climates. This study makes an important contribution to our understanding of variability in bark anatomical traits among Picea species and their ecological adaptations.

In woody plants, bark is an important protective tissue which can participate in photosynthesis, manage water loss, and transport assimilates. Studying the bark anatomical traits can provide insight into plant environmental adaptation strategies. However, a systematic understanding of the variability in bark anatomical traits and their drivers is lacking in woody plants. In this study, the bark anatomical traits of 23 Picea species were determined in a common garden experiment. We analyzed interspecific differences and interpreted the patterns in bark anatomical traits in relation to phylogenetic relationships and climatic factors of each species according to its global distribution. The results showed that there were interspecific differences in bark anatomical traits of Picea species. Phloem thickness was positively correlated with parenchyma cell size, possibly related to the roles of parenchyma cells in the radial transport of assimilates. Sieve cell size was negatively correlated with the radial diameter of resin ducts, and differences in sieve cells were possibly related to the formation and expansion of resin ducts. There were no significant phylogenetic signals for any bark anatomical trait, except the tangential diameter of resin ducts. Phloem thickness and parenchyma cell size were affected by temperature-related factors of their native range, while sieve cell size was influenced by precipitation-related factors. Bark anatomical traits were not significantly different under wet and dry climates. This study makes an important contribution to our understanding of variability in bark anatomical traits among Picea species and their ecological adaptations.

Introduction
In woody plants, bark refers to all tissues outside the vascular cambium and is an important component of the stem (Evert, 2006). Structurally, it is divided into two distinct parts: the outer bark (OB) and the inner bark (IB). The OB consists of dead cells, generally the rhytidome, and serves functions such as mechanical support and protection against pathogens. The IB, on the other hand, consists of living tissues including the phloem, and is responsible for the storage and transport of water and photosynthetic assimilates (Rosell et al., 2017). Bark plays important roles in the basic physiological functions of woody plants, and the sizes of tissues and cell morphologies in the periderm and phloem are determined by several developmental pathways and affect plant functions (Srivastava, 1964;Paine et al., 2010;Rosell, 2019). Bark morphology and structure are closely related to the physiological and ecological processes of plants. For example, the cortex in the bark contains chloroplasts capable of photosynthesis, which convert carbon dioxide produced by mitochondrial respiration and flowing in the xylem into sugars (Pfanz, 2008). It also increases bark oxygen concentration and reduces plant stem hypoxia (Wittmann and Pfanz, 2018). Furthermore, the phellem and lenticels in the bark structure regulate the exchange of water, oxygen, and carbon dioxide between the stem and its environment (Lendzian, 2006). For example, the phellem cells have suberin, a waxy substance that makes them impermeable to gases and water. Loram-Lourenco et al. (2022) found that the water conductance of bark across species was related to the morphoanatomical characteristics of the outer bark (i.e., thickness, density, and lenticel investment), while these outer bark characteristics were related to stem transpiration and respiration. For example, the phellem on the bark surface, which consists of dead cells, effectively prevents excessive water loss from the plant (Leite and Pereira, 2017). Therefore, bark anatomical traits are important to clarify their multifunctionality, resource allocation trade-offs, and environmental adaptative mechanisms (Rosell, 2019).
Environmental adaptability is apparent in the ecological strategies of most plant organs and tissues, including bark anatomical traits (Wright et al., 2004;Freschet et al., 2021;Yang et al., 2022). Differences in bark thickness and tissue structure can develop under stress, and such changes are associated with plant resilience (Kopanina et al., 2022). The phylogenetic niche conservatism hypothesis suggests that more closely related species are more likely to have similar functional traits and that conserved and similar functional traits will exhibit stronger phylogenetic signals (Blomberg et al., 2003;Losos, 2008). Plant anatomical traits can exhibit a wide range of variation due to climate-driven effects, for example, xylem vessels of plants in arid regions are often characterized by narrow lumens and thick walls (Baas et al., 1983). Many previous studies have shown that bark structural characteristics are related to fire factors (Lawes et al., 2013;Pausas, 2015) and that species from fire-prone areas tend to have thicker bark that better protects phloem tissues from destruction. Thus, it is important to study how bark anatomical traits vary due to environmental factors in the context of phylogeny.
There are about 40 Picea species worldwide, they are widely distributed in boreal, temperate, and subtropical high-altitude regions of the northern hemisphere (Farjon, 2001), making up significant portions of forests. Picea species provide ecological benefits and their bark is well utilized as a forest by-product (Harkin, 1971). For example, spruce bark extract is a natural antioxidant and anti-inflammatory agent. According to the information recorded by the Global Biodiversity Information Facility (GBIF, http://www.gbif.org), this genus has been widely introduced to various regions of the world and has strong environmental adaptability. Previous studies on Picea species have focused on the morphological and anatomical traits of wood (Piermattei et al., 2020;Puchi et al., 2020), pollen (Jia et al., 2014), and needles . Importantly, the organ or tissue anatomical traits of Picea species have been shown to be closely related to their physiological and ecological functions. For example, certain needle anatomical traits determine photosynthetic performance , and xylem cell number and cell lumen area affect the hydraulic systems of trees (Sperry and Tyree, 1988;Puchi et al., 2020). The structure and function of the bark influence water transport and storage in the plant. For example, low density bark has a higher hydraulic conductivity (Loram-Lourenco et al., 2020). Indeed, the organ or tissue anatomical traits in Picea can provide new insight into interspecific relationships and the mechanisms by which biotic and abiotic factors affect them. However, studies on the interspecific variation in bark anatomical traits and the underlying driving mechanisms in species of this genus are lacking.
In this study, the bark anatomical traits of 23 Picea species native to North America, Europe, and Asia were examined in conjunction with information on climatic factors of the Picea species. This was done to address the following questions: (1) are there interspecific differences in bark anatomical traits among Picea species and are anatomical traits phylogenetically conserved among species in different habitats; (2) are there trade-offs in changes among bark anatomical traits; and (3)  . The area has an elevation of 1,550-2,100 m, average annual temperature of 7.2°C, average annual precipitation of 757 mm, and average relative humidity of 78%, and the growing conditions are similar for all Picea species in the nursery. In October 2020, 23 Picea species from the experimental nursery were sampled. All species were sown in 2008. After three years of cultivation, individual trees were transplanted into the nursery with 1.5 m of spacing between them. For each Picea species, three single plants of uniform size and normal growth were selected and bark samples 2 × 2 cm in size were taken at 30 cm from the base of each plant. Bark samples include all tissues peeled from the cambium to the surface of the bark.

Bark sample fixation
Referring to the method of fixation in needles and pollen of Picea (Jia et al., 2014;Wang et al., 2021), bark sections were prepared using the following 10 steps: (1) Sample fixation: Bark samples were fixed with FAA fixative (90 mL ethanol + 5 mL formaldehyde + 5 mL acetic acid) for 24 h. (2) Sample dehydration waxing: Cut 2-3 mm bark samples along the tangential direction and rinsed in running water for 30 min, placed in 15% ethanol for 2 h, and transferred to a dehydrator (DIAPATH, Donatello) for dehydration. (3) Sample embedding: Melted wax was poured into an embedding frame and before the wax solidified the tissue was removed from the dehydration box and placed into the frame according to the embedding surface. Finally, the samples were cooled on the -20°C freezing table to solidify the wax, which was then trimmed. (4) Sample sectioning: Trimmed wax blocks were placed in a paraffin slicer (Shanghai Leica Instruments Co., Ltd., China, RM2016) and sectioned at a thickness of 4 mm. (5) Sample dewaxing: Sections were dewaxed and hydrated using ethylene glycol monoethyl ether acetate, ethanol solution, and distilled water. (6) Safranin O staining: Sections were placed in safranin O staining solution for 2 min, and then washed briefly in distilled water to remove excess dye. (7) Decolorization: Sections were placed sequentially in a 50%, 70%, and 80% alcohol gradient for 3-8 s each, in order to wash away the excess safranin O staining solution. (8) Fast green staining: Sections were placed in fast green staining solution for 6-20 s and dehydrated in anhydrous ethanol. (9) Sample sealing: Sections were placed in xylene for 5 minutes and sealed with neutral balsam. (10) Microscopic examination: Bark samples were observed using an optical microscope (Nikon Eclipse E100, Japan) and imaged with a high-definition camera (Nikon DS-U3, Japan).

Data acquisition
We selected three bark samples for each species as biological replicates. Considering the differences between primary and secondary structures in the bark and the difficulty in differentiating them, 30 parenchyma cells and sieve cells were randomly selected from the scanned image of each sample to measure their radial and tangential widths. In this study, bark anatomical sections were measured using the CaseViewer software (3DHISTECH CaseViewer, Budapest, Hungary). The anatomical trait related to bark structure and function were selected for assessment (Angyalossy et al., 2016;Schweingruber et al., 2019;Rosner and Morris, 2022), including (Table 1): periderm thickness (PE), cortex thickness (CO), phloem thickness (PH), tangential diameter of resin duct (TR), radial diameter of resin duct (RR), phloem ray width (PR), tangential diameter of sieve cell (TS), radial diameter of sieve cell (RS), tangential diameter of parenchyma cell (TP), and radial diameter of parenchyma cell (RP). In order to reduce the effects of uneven growth of bark cells and tissues, the aspect ratios of some anatomical traits were also calculated, including tangential diameter of resin duct/radial diameter of resin duct (TR/RR), tangential diameter of sieve cell/radial diameter of sieve cell (TS/RS), and tangential diameter of parenchyma cell/radial diameter of parenchyma cell (TP/RP). A detailed bark anatomy is provided in Figure 1.
To further understand the responses and adaptations of bark anatomical traits of Picea species to climatic factors, we collected information regarding all original collection sites of the 23 Picea species by referencing Ouyang et al. (2021) and GBIF (http://www.gbif.org) ( Figure 2 and Table 2). Each species averaged all sites within the GBIF record range to represent the center position of the sampling site. The temperature and precipitation related climate factors were also obtained from WorldClim v2.0 (http://www.worldclim.org) in conjunction with the geographic information of the sources. These included the annual mean temperature, mean temperature of the wettest quarter, mean temperature of the driest quarter, annual precipitation, precipitation of the wettest quarter, and precipitation of the driest quarter. The global aridity index was obtained at the CGIAR consortium for spatial information (CGIAR-CSI, https:// cgiarcsi.community) (Zomer et al., 2022). In addition, based on mean annual precipitation, we classified areas with > 500 mm as moist and areas with < 500 mm as dry (Smith et al., 2008). The climate variables for each species were averaged across all source sites for the species using the 'raster' package in R 3.6.3 software (R Core Team, 2020).

Statistical analysis
The 'V.PhyloMaker' package (Jin and Qian, 2019) and the 'Plantlist' package (Zhang, 2018) in the R 3.6.3 software were used to elucidate the influence of phylogenetic relationships on Picea species attributes. The Pagel's l and Blomberg's K values (Harmon et al., 2008;Kembel et al., 2010) of bark anatomical traits  Cross section of bark anatomical traits, exemplified by Picea engelmannii. Abbreviations: ca, cambium; co, cortex; ib, inner bark; ob, outer bark; pc, parenchyma cell; pd, phelloderm; pe, periderm; pg, phellogen; ph, phloem; pl, phellem; pr, phloem ray; rd, resin duct; rh, rhytidome; sc, sieve cell; sp, suberized filling tissue or phellem cells with polyphenolic content; uf, unsuberized filling tissue; x, xylem. were calculated using phylogenetic trees to assess the phylogenetic signal of bark anatomical traits in Picea species. To further elucidate the relationships among bark anatomical traits, phylogenetic independent contrasts (PIC) were calculated for each bark anatomical trait using the 'ape' package (Felsenstein, 1985). General Linear Model (GLM) with ANOVA and Duncan's multiple range test were performed in SAS 9.4 (SAS Institute Inc., Raleigh, NC) to assess differences in bark anatomical traits among Picea species. Pearson's correlation analysis and principal component analysis (PCA) were used to assess the relationships between different bark anatomical traits and their associations with climatic factors of their native range Meng et al., 2022). Independent sample T-tests were used to analyze differences in bark anatomical traits under dry and moist conditions. The phylogenetic tree of this study was plotted online in ChiPlot (http://www.chiplot.online/). All statistical analysis was conducted in the R 3.6.3 software unless specified and P < 0.05 was considered statistically significant.

Bark anatomical traits and their interspecific differences
The bark of the 23 Picea species consisted of similar tissues (Figure 1 and Supplementary Figure 1). The outer bark is mainly consisted of a thick periderm and a thin rhytidome. The periderm had suberized filling tissue (A looser tissue arising outward from the phellogen in the lenticels) or phellem cells with polyphenolic content and also included unsuberized filling tissue. Between the cortex and the periderm was phellogen, with more neatly arranged cells. The cortex included resin ducts, sieve cells, and parenchyma cells, and the secondary resin ducts were mainly composed of about 2-3 layers of epithelial cells. The phloem was composed of multiple layers of parenchyma cells, sieve cells, and radially distributed phloem rays, however, the sieve cells in the secondary phloem collapsed (Figure 1). Among all the observed species, P. pungens had the largest RR (384.92 mm), while P. sitchensis had the largest PH (1614.60 mm). P. obovata had the smallest TP (32.59 mm) ( Table 3). The resin ducts size of Picea species ranges from about 132.32-803.80 mm, with large variation. ANOVA and Duncan's multiple range test showed that there were significant differences in bark anatomical traits among the 23 Picea species (Tables 3, 4). The principal component analysis showed that the first principal component was loaded mostly by PH and CO, while the second principal component was loaded mostly by TR, RR, and PR (Figure 3).

Correlations between bark anatomical traits
The present study showed that there were significant correlations between the bark anatomical traits of Picea species (Figure 4). There was a significant positive correlation between phloem thickness and parenchyma cell size (P < 0.05). There was a significant negative correlation between sieve cell size and RR (P < 0.05). The aspect ratio of parenchyma cells was significantly and positively correlated with sieve cell size (P < 0.05). There was a highly significant positive correlation between TS and PH (P < 0.001).

Phylogenetic signals of bark anatomical traits
The analysis of phylogenetic signals showed that TR had a weak phylogenetic signal (Blomberg's K = 0.47, P < 0.05). However, none of the other bark anatomical traits of Picea species had significant phylogenetic signals (Table 5 and Figure 5).

Relationships between bark anatomical traits and climatic factors
The bark anatomical traits of Picea species were influenced by climatic factors of their native range ( Figure 6). PR exhibited a positive correlation with the annual mean temperature (P < 0.05), while both PH and TP were positively correlated with the mean temperature of the driest quarter (P < 0.01). TP was positively correlated with annual precipitation (P < 0.05). However, there was a negative correlation between TS and precipitation of the driest quarter (P < 0.05). For all Picea species, TR was not correlated with any climatic factor (P > 0.05). These results clearly indicated that climatic factors of their native range had driving effects on bark anatomical traits. In addition, we found no significant differences (P > 0.05) in bark anatomical traits between dry and moist conditions (Supplementary Figure 2). This study suggested that the PH and TP of bark in Picea species are driven by temperaturerelated factors of their native range while TS is influenced by precipitation-related factors, but all traits are relatively insensitive to moist and dry environments at the global scale.

Discussion
In this study, we assessed the interspecific patterns in bark anatomical traits of Picea species and analyzed the phylogenetic relationships and climatic drivers underlying the differences in bark anatomical traits. The large variability in bark anatomical traits among the 23 Picea species suggested that species had a significant influence on bark anatomical traits. In addition, some of the bark anatomical traits exhibited synergistic relationships, there were no obvious trade-offs between traits. In general, bark anatomical traits showed only a slight influence from phylogeny, but strong influences from climatic factors of their native range.

Functions and interspecific variation in bark anatomical traits
It is generally accepted that tree phenotypes are shaped by a combination of genetic, developmental, and environmental factors  (Coleman et al., 1994). Our study revealed specific interspecific differences in anatomical traits of the outer and inner bark, suggesting that the bark performs different functions among species. Sieve cells mainly perform the axial transport of assimilates, whereas parenchyma cells perform the radial transport of assimilates (Deng et al., 2020). Different bark anatomical structures determine bark functional traits, for example, differences in the radial structures of the phloem lead to discontinuities in the radial distribution of in-situ water content and saturation osmotic potential (Rosner et al., 2001). Bark is multifunctional and involved in the transport of photosynthetic assimilates, transpiration of plants, mechanical support, and defense against fire, insects, and pathogens. Bark often has lenticels that regulate water loss under dry conditions (Lendzian, 2006). Studies have shown that bark transpiration under drought conditions can be the source of more than half of the water loss from the plant and that transpiration from bark is generally a passive process not associated with plant metabolism (Lintunen et al., 2021). This water dissipation function tends to be more closely related to the outer bark (Loram-Lourenco et al., 2022), while water storage and transport are more related to the structure and function of the inner bark (Loram-Lourenco et al., 2020). Picea species typically have lenticels in the outer bark associated with water loss (Rosner and Kartusch, 2003), but lenticel anatomy was not studied here. Cortex contains chloroplasts that perform photosynthesis, reduce CO 2 production by the stem, and prevent acidification of the cortex (Pfanz, 2008). Bark tissues have ecological strategies appropriate to the environment of different locations, with thicker bark occurring in fire-prone areas and thinner bark in tropical areas (Paine et al., 2010). Similarly, bark in areas with high insect infestation rates exhibited a denser composition and other induced defense strategies (Franceschi et al., 2005) associated with specific bark structures. In addition, the structure of the phloem in the bark changes with age. For example, sieve cells tend to accumulate and take on irregular shapes. Furthermore, the thickening of the cell walls of parenchyma cells can turn them into stone cells, thus halting their cellular activity (Schweingruber et al., 2019). All species in our study were 12 years old, which controlled for any error caused by differences in age. There were significant interspecific differences in bark anatomy, which is similar to the results of interspecific differences in needle anatomy of Picea . All growth conditions in our common garden experiment were consistent, so we concluded that genetic factors had an important influence. Correlation analysis and principal component analysis of bark anatomical traits revealed significant positive correlations among most anatomical traits, correlations which were consistent even when the phylogenetic independent contrasts were resolved (Supplementary Figure 3). This indicated there were strong correlations among different bark anatomical traits. Both phloem thickness and parenchyma cell size were significantly and positively correlated, which may have been related to the involvement of parenchyma cells in the radial transport of assimilates. Sieve cell size was significantly negatively correlated with the radial diameter of resin ducts, and the variation of sieve cells might have been related to the expansion of resin ducts (Esau, 1969). Previous studies have shown that the swelling of the sieve cells was accompanied by a decrease in the number of resin ducts, which had a positive effect on the resistance of the plant (Krokene et al., 2008). Our results demonstrated that bark anatomical traits were generally located along the same axis in the principal component analysis, which indicated that there were no trade-offs among them, but mostly synergistic or complementary relationships. Thus, this study provides new insights into the unique position in which the bark economics spectrum (BES) is located in the plant economics spectrum (PES) (Li et al., 2022).

Drivers of variation in bark anatomical traits of Picea species
Except for TR, there were no significant phylogenetic signals for bark anatomical traits in any of the 23 Picea species. This indicated that phylogeny has little influence on the variation in bark anatomical traits, which was similar to the observation that bark thickness traits of angiosperms were not significantly influenced by phylogeny (Rosell et al., 2014). Our analysis indicated that climatic factors of their native range have a strong influence on bark anatomical traits. Martıń-Sanz et al. (2019) found that dry conditions were not conducive to bark biomass partitioning and thus led to thinner bark, which is generally consistent with the positive correlation between the precipitation of the wettest quarter (PREWQ) and phloem thickness in our study. The mean temperature of the driest quarter had a strong effect on parenchyma cell size, which increased significantly with increasing temperature. There can also be seasonal differences in saturated osmotic pressure in the secondary phloem (Rosner et al., 2001), which may be related to the climate-driven nature of the phloem structure. Previous studies have revealed associations between parenchyma cell size and the partitioning of nonstructural carbohydrates, a mechanism of assimilate partitioning that regulates osmotic pressure in the phloem and is associated with plant resistance to embolism (Janssen et al., 2020). Poorter et al. (2014) concluded that there were no significant Principal component analysis of bark anatomical traits. For trait abbreviations, see Table 1. Nie et al. 10.3389/fpls.2023.1201553 Frontiers in Plant Science frontiersin.org differences in bark characteristics between dry and moist forests, however, our results suggested that bark anatomical traits were more strongly driven by temperature than moisture conditions. Previous studies have shown that resin ducts are associated with conifer defense functions, with wide resin ducts providing more resin and more effective protection to conifers (Mason et al., 2019). It has also been shown that the number of resin ducts is strongly correlated with temperature (Novak et al., 2013), which is consistent with our finding of a weak correlation between resin duct size and temperature. There was a weak negative correlation between resin ducts and the precipitation of the wettest quarter (PREWQ) in our study, suggesting that resin production may be reduced under less stressful conditions. Bark surface insect activity and microbial composition are often linked to bark anatomical traits. For example, the bark structure of inverted wood is associated with invertebrate and microbial composition within communities, and the greater the overall variation in bark traits, the greater the variation in faunal community composition and species richness on its surface (Zuo et al., 2016). The cortex and phloem are mostly living tissues consisting of multiple layers of parenchyma cells and sieve cells in which chemicals such as phenolics, terpenoid resins, and alkaloids provide a second layer of protection for the tree (Franceschi et al., 2005). In conclusion, bark exhibits corresponding physiological and ecological strategies when disturbed by biotic and abiotic factors.

Conclusions
Our common garden experiment revealed large interspecific differences in bark anatomical traits. Bark anatomical traits were not phylogenetically influenced but were influenced by climatic conditions of the plant origin. Interspecific variation in bark Correlations among bark anatomical traits of 23 Picea species. *: P < 0.05, **: P < 0.01, ***: P < 0.001. For trait abbreviations, see Table 1.  Phylogenetic relationships of Picea species and the variability in their bark anatomical traits. Different heatmap colors represent different bark anatomical trait sizes. All bark anatomical traits were normalized to 0-1. For trait abbreviations, see Table 1. Correlation between bark anatomical traits and climatic factors in 23 Picea species. PE, periderm thickness; CO, cortex thickness; PHm phloem thickness; TRm tangential diameter of resin duct; PR, width of phloem ray; TS, tangential diameter of sieve cell; TP, tangential diameter of parenchyma cell; AMT, annual mean temperature; MTWQ, mean temperature of the wettest quarter; MTDQ, mean temperature of the driest quarter; APRE, annual precipitation; PREWQ, precipitation of the wettest quarter; PREDQ, precipitation of the driest quarter; AI, aridity index. All climatic factors were acquired in their native range. *: P < 0.05, **: P < 0.01.
anatomical traits is an important driver of multifunctional variation in bark. It appeared that all bark anatomical traits contributed in the same dimension, influencing to the mechanical and chemical defenses of the plant. In addition, bark anatomical traits were driven by climatic factors at the seed source, which may indicate long-term climatic adaptation by the plants. Our study provided important new insight into the variability of bark anatomical traits in Picea, but this effort should be followed by controlled experiments to elucidate the plasticity of bark structures and use molecular techniques to reveal the functions and mechanisms of bark structural variation.

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 authors.

Author contributions
WN contributed to execute the experiment, data analysis, manuscript drafting. YD, YL, CT, YW, and YY conducted field data collection. WX, JL, JM, and SA provided experimental guidance. ZPJ, ZRJ, and JW designed the study. All authors contributed to the article and approved the submitted version.

Funding
This study was supported by the National Promotion Project of Forestry and Grassland Achievements (2020133110) and the National Natural Science Foundation of China (31500540).