The samples were collected from the experimental nursery of outdoor in the Shaba experimental base of the Research Institute of Forestry of Xiaolong Mountain in Gansu Province, China (104°38’E, 34°07’N). 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.

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 μm. (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).

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

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 λ and Blomberg’s K values (Harmon et al., 2008; Kembel et al., 2010) of bark anatomical traits 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 (Ouyang et al., 2021; 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.

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 μm), while P. sitchensis had the largest PH (1614.60 μm). P. obovata had the smallest TP (32.59 μm) ( Table 3 ). The resin ducts size of Picea species ranges from about 132.32–803.80 μm, 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 ).

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

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

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 temperature-related 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.

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-Lourenço et al., 2022), while water storage and transport are more related to the structure and function of the inner bark (Loram-Lourenço 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₂ 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 (Wang et al., 2021). 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).

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ín-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 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.