We compiled wood data from 794 Chinese tree species (Yang et al., 2001, 2009). Most tree samples were taken from the Eastern Monsoonal climate zone in China, where the temperate, subtropical, and tropical natural forests are distributed (Wu, 1980). The anatomical data was mainly collected from publications by The Chinese Academy of Forestry, which hosts the largest wood collection in China, and the wood anatomical traits were measured following their methods. In brief, most angiosperm tree species were sampled from mature trees with a DBH over 20 cm in natural forests. The sampled disks were collected at a height of 1.3m. These disks were cut, from pith to bark, into six equal parts according to the equidistance method and anatomical samples were taken from the middle of the outermost part (the closest to the bark). The entire span of one growth ring was sectioned with a microtome. Wood anatomical traits were measured using a light microscope and analyzed with an image analysis system (Q570). For vessel diameter and vessel density, 100 vessels were measured in earlywood and latewood of each sample. Other important anatomical traits measured were the proportion of cross sectional area occupied by vessels, fibers (including tracheid), rays, and axial parenchyma. To measure tissue fractions of cross sections, images were taken with a Q570 image analyser. The percentage area of each cell type, however, was manually determined for 1 mm² (Yang et al., 2001). All species names were verified using the Plant List¹ to correct for synonyms, while synonyms were removed. In total we obtained wood anatomy data for 780 woody angiosperms species (Appendix 1). The anatomical data were deposited in TRY² (Kattge et al., 2011).

As the sampling locations were not recorded for most of the species in the database, climate data for the natural distribution range in China were retrieved for each species. Among the 780 species in our database, climate information for 700 species over their range were extracted from the Atlas of Woody Plants in China (Fang et al., 2009). This book provided mean values for thirteen climatic variables across the range of each species on county level, as well as the species’ distribution map at the country level. The climate variables included three groups. (1) The thermal indices: mean annual temperature (MAT, °C), annual biotic temperature (ABT, °C), potential evapotranspiration (PET, mm), Kira’s warmth index (WI, °C ⋅ month) and coldness index (CI, °C ⋅ month), mean temperature of the warmest month (MTWM, °C), and mean temperature of the coldest month (MTCM, °C). (2) The humid/arid indices: mean annual precipitation (MAP, mm), precipitation in the warmest quarter of the year (PWQ, mm), and precipitation in the coldest quarter of the year (PCQ, mm). (3) Integrative indices included annual actual evapotranspiration (AET, mm) and moisture index (Im), and vegetation net primary production (NPP, g.a^-1.M^-2) estimated using the CASA model (for details, see Fang et al., 2009). In addition, maximum and minimum latitude and longitude for each species were also extracted from the distribution map to calculate the latitude and longitude midpoints of each species’ range (Appendix 1). Elevation could not be integrated in our analyses because we largely lack altitudinal gradients in our dataset. However, since most of the species of our dataset were distributed in the eastern part of China where altitudinal gradients are largely limited, it is reasonable to assume that altitudinal gradients play a less significant role than latitudinal gradients, especially under temperate climate in China. Only 19 species with a major distribution in western part of China were included in our dataset (e.g., Quercus semecarpifolia, Rhododendron delavayi, Sorbus tianschanica, Populus euphratica).

We build a phylogenetic tree for 780 species in this study using the R package “brranching” (Chamberlain, 2016). The stored super-tree was “R20120829” and the branch lengths of the phylo-tree were calculated using the BLADJ algorithm in phylocom (Webb et al., 2008). The tissue fractions (refereed as CFs afterward) were square root-transformed to increase normality in the data, before statistical analyses were performed with R.3.2 (R Core Team, 2015). Based on the phylo-tree obtained (Appendix 2), phylogenetic signals of quantitative traits were calculated as Pagel’s lambda using the “phylosignal” function in the R package “phytools” (Revell, 2012), phylogenetic independent contrasts for CFs were calculated using the AOT algorithm in phylocom (Webb et al., 2008), and the association between PICs for each cell type was analyzed using the “corr.test” function with the R package “psych” (Revelle, 2017).

To describe different aspects of xylem architecture, we calculated the F and S vessel metrics according to Zanne et al. (2010). F, which is the product of the mean vessel size (A, mm²) times the vessel number per unit area (N, mm^-2), measures the wood fraction that is occupied by the vessel space, which is interpreted as directly measured vessel tissue fraction in this paper. Increases in F should be correlated with lower mechanical strength (Zanne et al., 2010). S is the ratio of the same anatomical traits (A/N; mm⁴), and measures the variation in vessel composition. A higher S reveals a greater contribution of large vessels to water conduction in a given area, and therefore indicates increased water transport efficiency, but brings a heightened risk of frost-induced embolism. These two metrics represent orthogonal axes of variation (Zanne et al., 2010). All indices were log₁₀-tranformed to increase their normality for statistical analyses.

We conducted two types of phylogenetic analyses to account for the possible influence of phylogeny on the results that may explain the relationship between CFs and climatic factors. Firstly, the phylogenetic generalized least square method (PGLS) was used to build bi-variate and multi-variate models in a phylogenetic context, with the “gls” function in the R package “nlme” (Pinheiro et al., 2016), in which each tissue fraction was expressed as a function of climatic variables. In this method, the influence of phylogeny on a trait was included in the models’ correlation structure as a simple Brownian motion, thus various climatic variables for a trait could be compared within the same correlation structure. As there were 13 variables as potential predictors of each CF, we needed to reduce the number of predictors to minimize collinearity. We did a three-step procedure to estimate the combined effect of climate and traits on tissue fractions (Zheng et al., 2017). As most of the climatic variables are closely inter-correlated (Appendix 3), we started by reducing the initial pool of 13 climatic variables to a smaller number based on the results of the bivariate analysis. Firstly, the variables with an R² < 0.02 in explaining CFs were excluded (see Appendix 4). Next, within each of three groups of climate indices, we used an AIC-based backward elimination to select the variables that best explained CFs variation. Lastly, all climatic predictors that survived the first two steps were included to build the full models, and AIC was used to select the most parsimonious model. As there are too many possible interactions that would overwhelm our available degrees of freedom in the full models, we did not include any interactions. Determination coefficients (R_phy²) for each model were calculated using the maximum likelihood (ML) estimation. For our phylogenetic analysis, R_phy² was defined as the proportion of variation explained by the linear model, taking the correlation among observations into account, divided by that of the null model (intercept-only) with the same correlation structure in PGLS. Therefore, R_phy² can be regarded as a pseudo-R².

For the second phylogenetic test, we used a mixed-model analysis of covariance (MACO) to partition contributions from phylogeny and climate (He et al., 2009). In contrast to the commonly used method of phylogenetically independent contrasts or phylogenetically general least square, this method has the advantage that the effect of phylogeny on CFs can be quantified as the percentage of variations explained by ancient lineages formed at different times in the evolutionary history. Using this method, we created our phylogenetic groups by “cutting” the tree at a given divergence time, since taxonomic groupings may adequately represent the pattern of evolutionary relationships between species, but not represent the rate of evolution. We created two contrasts, a first-order grouping cut at 59 million years before present (n = 85 groups, Appendix 2), and a second-order grouping cut at 28 million years before present (n = 292 groups). We chose these ages because they resulted in roughly the same number of groups, and thus degrees of freedom, as the families and genera used in the taxonomic analysis. In brief, this method was exactly analogous to the decomposition of the sums of squares in an analysis of variance (ANOVA), and the total sums of squares was partitioned into component sums of products for the above relationships. F-ratios of mean products ( = sum of products/degree of freedom) were used to test for significance. Using an ordinary rather than REML mixed-model analysis of covariance, we could obtain the sum of products that were additive, allowing us to express the influence of phylogeny and climate on CFs by the percent of the total sum of the products they explained (for details, see He et al., 2009). Consequently, the relative effects of phylogenetic and climatic factors on CFs could be compared (He et al., 2009).

Xylem tissue fractions varied greatly among the 700 angiosperm woody species. For vessel fraction, values range from 8.19 to 66.50%; for fiber fraction, values range from 63.08 to 92.74%; for ray parenchyma fraction, values range from 22.87 to 65.17%; for axial parenchyma fraction, values range from 0 to 54.31%; for total parenchyma fraction (sum of ray parenchyma and axial parenchyma fraction), values range from 26.23 to 70.34%. Frequency of values for axial parenchyma fraction is strongly right skewed while other four fractions were less skewed (Figure 2).

All five xylem tissue fractions, i.e., vessel, fiber, ray, axial parenchyma, and total parenchyma (ray and axial parenchyma combined), showed significant phylogenetic signals (P < 0.001). Among the five tissue fractions, axial parenchyma had the highest value of lambda (0.895) followed by vessels (0.852), total parenchyma (0.760), fibers (0.732), and ray parenchyma (0.673).

Tissue fractions of vessels, fibers, ray, and total parenchyma showed significant latitudinal patterns with and without taking phylogeny into account (Figure 3); and the AICs were lower for PGLS models than corresponding linear models (Appendix 4). Vessel fraction increased (0.67%∼44.22%) with rising latitude (N18.6°∼N72.1°), while tissue fractions of fibers, ray, and total parenchyma were higher at lower latitudes. No significant linear relationship was detected between axial parenchyma fraction and latitude (P = 0.718 and 0.053 for the phylogenetic and linear model, respectively). The other two variables describing vessel composition (S) and vessel lumen area (F) showed opposite latitudinal pattern, with S decreasing while F increasing with a rising latitude (Figure 3).

In the bi-variate models for tissue fractions and climatic variables, the majority of the 13 climatic indices correlated positively with vessels but negatively with fibers, ray, and total parenchyma fractions; however, the correlation between axial parenchyma fractions and any climatic index was not significant (Figure 4 and Appendix 4). MAT (mean annual temperature) and MTCM (mean temperature of the coldest month) had higher R_phy² than other temperature indices. AET (annual actual evapotranspiration) was the most influential index in integrative climatic indices, while three precipitation indices had less explanatory power compared to temperature indices for all tissue fractions in term of R_phy² (Appendix 4).

Among the top multivariate PGLS models with different tissue fractions as a function of climatic variables, MAT, MTCM, and AET were the only three climatic variables entered, while precipitation indices were excluded (Table 1). MAT was included in the top multivariate models for vessels, fibers, ray, and total parenchyma fractions, but not for axial parenchyma. AET was included in the models for both vessel and fiber tissue fractions, and MTCM for vessel, ray parenchyma, and total parenchyma tissue fractions. Climatic factors explained a larger proportion of vessel fraction variances (R_phy² = 0.401) than those of other tissue fractions.

Decomposition of covariance showed that phylogeny and climate together could explain about 59.7–72.8% of the variation for different tissue fractions (Figure 5). However, the explanatory power of climatic parameters only was very small, while phylogeny could explain about 50.8–72.2% variance of the tissue fractions (Appendix 5). For vessel fraction, phylogeny alone accounted for 50.8% of the variance, with climate alone accounting for 5.3% of the variance, and phylogeny and climate sharing 11.1% of the variance.

Opposing trends between vessel and non-vessel tissue fractions along latitudinal gradients were found among angiosperm species across China, which supports our hypothesis. These spatial pattern of xylem tissue fractions also suggested that xylem space allocation shows at least partly adaptation of the plant hydraulic system in woody plants. Therefore, the trade-off between non-vessel and vessel tissue fractions reflects the hydraulic properties of the vessel tissue and the additional xylem tasks such as hydraulic capacity (Jupa et al., 2016), storage of non-structural carbohydrates by parenchyma (and living fibers) (Plavcová et al., 2016), and protection against pathogens (Morris et al., 2016a). On one hand, trees at high latitudes generally possess small vessel diameter to reduce the risk of frost-induced hydraulic failure (Zanne et al., 2014). Increased vessel density, however, and relatively high vessel tissue fractions could mediate reduction of water conductivity for trees at high latitudes. Since hydraulic conductivity scales to the fourth power of the vessel diameter, but linearly with vessel density according to the Hagen–Poiseuille law, it is possible that increased vessel density alone is not enough to compensate the effect of narrow vessel diameters. Moreover, vessel diameter is well known to be negatively correlated with vessel density. On the other hand, reduced fiber and ray fractions in trees located at higher latitudes might be associated with reduced growth rate or plant size compared to trees in tropical regions. However, changes in fiber and ray tissue fractions would not imperil xylem’s multiple functions that are essential for a plant’s life history under various climatic conditions, as diversified xylem space allocation in angiosperms could be viewed as a strategy to adapt to different environments (Bittencourt et al., 2016). For instance, wide variance of tissue fractions occurred in a rather narrow range of wood density for 69 Australian tree species (Zieminska et al., 2015). No clear latitudinal pattern of axial parenchyma was found, which is probably due to a lack of tropical tree species in our data. Levels of axial parenchyma in tropical trees are generally higher than those of temperate trees, and on a global scale axial parenchyma fraction declines along latitude (Morris et al., 2016b).

The rising vessel fraction toward higher latitudes came to a great extent at the cost of decreased fiber fractions, which was associated with a change in total parenchyma fractions. PIC correlation analysis indicated that a strong association existed between vessel fraction and fiber fraction (r_pic = -0.70), and between fiber fraction and total parenchyma fraction (r_pic = -0.63), and only weak association between total parenchyma fraction and vessel fraction (r_pic = -0.08) (Table 2). Therefore, a major trade-off between vessel and fiber fraction occurred along a latitudinal gradient. This is in agreement with previous studies, which also showed that strong trade-offs existed between vessel and fiber fractions within the xylem space (Fichtler and Worbes, 2012; Zheng and Martínez-Cabrera, 2013; Ziemińska et al., 2013). It has been suggested that structural trade-off between parenchyma and fiber fraction is due to the relatively high proportion of both tissues in wood (Ziemińska et al., 2015). Our study agrees with the idea that a strong trade-off exists between parenchyma and fiber fractions. Taken together, these results suggested that a changing environment along a latitudinal gradient in China has a significant effect on the structural trade-offs between the three main cell types (vessel, fiber and parenchyma). Moreover, the structural trade-offs are directly linked to functional trade-offs, although each cell type may fulfill multiple functions rather than a specific function.

To further explore the changes in vessel composition along a latitudinal gradient, we analyzed the trends of the vessel composition index (S, mm⁴) and vessel lumen area index (F, unitless) along a latitudinal gradient. Zanne et al. (2010) showed that S could explain a large proportion of structural variance within the hydraulic system, and that variation in conductivity is more strongly influenced by S than by F (Zanne et al., 2010). Results showed that S decreased, whereas F increased significantly with rising latitude (Figure 3). The latitudinal trend of F is consistent with that of vessel tissue fraction, and both metrics have the same meaning, although F is calculated by multiplying vessel diameter measurements with vessel density data, while vessel fraction is measured directly. The higher value of S at low latitudes indicates that the conducting area is comprised of a few large vessels, which are highly efficient for long distance water transport, while low values of S at high latitudes indicate that the conducting area is comprised of many small vessels, which are resistant to embolism risk under freeze-thaw conditions (Zanne et al., 2014). Therefore, the opposite latitudinal trends of S and F suggest a trade-off between vessel size and vessel number in the hydraulic system of woody angiosperm species across latitude.

Climate could have an effect on xylem space use as most climatic indices were significantly associated with vessel, fiber, ray, and total parenchyma fractions, except for axial parenchyma fraction (Appendix 4). Generally, thermal indices have better explanatory power than moisture indices. For instance, in two bi-variate models explaining variance in vessel tissue fraction by mean annual temperature (MAT) or by mean annual precipitation (MAP), R²_phy of the former model (0.367) was higher than that of latter (0.301). A global study on parenchyma patterns found that MAT had a stronger influence on total parenchyma fractions than MAP, but no significant relationship was found between axial parenchyma and these two climatic indices (Morris et al., 2016b), which is in line with our results. We found that climatic variables also influence vessel and fiber fraction (Appendix 4). Moreover, the moisture index (Im), an indicator for regional aridity, was only weakly associated with vessel fraction (R²_phy = 0.058), showing that an aridity gradient is not obvious across the latitudinal range in China. Most of the woody angiosperms species in this study were distributed in the eastern part of China, where a monsoonal climate dominates with a synchrony of thermal and moisture conditions. This is consistent with the fact that water was generally not a limiting factor for trees, even in regions that experience freezing events (Bjorklund et al., 2017).

Temperature had an important influence on trade-offs between xylem tissue fractions along altitudinal gradients. The top multivariate models, which include MAT, MTCM and AET, support this point (Table 1). Temperature was a critical factor in limiting the production and differentiation of xylem cells in cold climates (Rossi et al., 2012). For instance, mean annual temperature and coldest temperatures were correlated with maximum plant height globally (R² = 0.170∼0.195) (Moles et al., 2009). Our results showed that a correlation between MCTM and vessel and parenchyma fractions was higher than those between MAT and vessel and parenchyma fractions in relevant bi-variate models (Appendix 4), suggesting that winter temperature may influence xylem development, and thus xylem structure and functions. According to the growth limitation hypothesis (Körner, 1998), there is a minimum temperature that permits sufficient production of new cells and differentiation of xylem tissue in higher plants. An average daily temperature of 6 to 8°C is commonly considered to be a thermal boundary layer below which metabolic activities are inhibited (Rossi et al., 2012). In our dataset, c. 300 species had MTCM values lower than 8°C, thus explaining the close association between MTCM and vessel fraction as well as fiber, ray and total parenchyma fractions, which is in agreement with the growth limiting hypothesis. Furthermore, recent studies on wood formation showed that temperature could regulate the timing of cambial reactivation and xylem differentiation for temperate trees, thus temperature could regulate the important xylem and phloem functions of transport, storage, and mechanical support (reviewed by Begum et al., 2018). In addition, a decrease in temperature might be the critical factor in the control of latewood formation and the cessation of cambial activity (Begum et al., 2018).

In this paper, xylem tissue fraction demonstrates a clear pattern along a large latitudinal gradient, where phylogeny alone could explain approximately half of all variance in the traits studied. Our results thus suggest that xylem anatomy is conservative, being largely determined by evolutionary history rather than ecological conditions. We found it unexpected that climate alone has a respectively low explanatory power over tissue fraction variance given that bi-variated models indicated a significant correlation between climatic indices and tissue fractions in general. Commonly, climate had a strong influence on plant performance (such as size, phenology, productivity, and some hydraulic traits) on a global scale. For instance, theoretical hydraulic conductivity (K_th) of 28 Eucalyptus species across Australia was significantly associated with aridity index (Pfautsch et al., 2016), while xylem vulnerability to embolism of Fagus sylvatica across Europe was associated with mean annual temperature and aridity (Schuldt et al., 2016). One possible reason for the low explanatory power of climatic indices on xylem tissue fractions could be that the climatic indices used in this paper represented a species’ average climatic condition over ones range rather than actual sample sites, while trees’ performance (including hydraulic function and xylem structure) might differ between their geographical center and margin (Stojnic et al., 2018).

Another possible reason for the low explanatory power of climate on xylem tissue fractions might lie in the modularity of angiosperms. Higher plant species grow through repeated increases in their modules (e.g., leaves, branches), so the effect of climate on the functions at the individual component level is likely compromised by the tree’s overall architecture. In commenting on why only a weak trade-off between xylem safety and efficiency was detected at the branch level in a global study by Gleason et al. (2016b), Brodersen (2016) suggested that hydraulic variables of stems and branches might not be enough to represent the whole tree’s hydraulic functioning, and integration of hydraulic data from along the entire hydraulic pathway would be useful. The same concept of modularity has also been described under different terms, such as hydraulic vulnerability segmentation (Johnson et al., 2016), hydraulic sectoriality (Zanne et al., 2006), and hydraulic redundancy (Ewers et al., 2007). Future studies should integrate structural anatomy of xylem and its multiple functions at different plant levels, such as at the cellular level (i.e., the role of pit membranes in water transport; Jansen et al., 2018), the tissue level (i.e., connection between different cell types, especially parenchyma and vessels), the branch level (i.e., autonomy of branch functions), and at the whole-plant level (i.e., the xylem continuum from minor roots to minor leaf veins).