Inﬂuence of Colder Temperature on the Axial and Radial Parenchyma Fraction of Quercus ciliaris C.C.Huang & Y.T.Chang Wood and Its Relationship with Carbohydrate Reserve (NSC)

: Parenchyma in the secondary xylem comprises the main tissue for the storage of nonstructural carbohydrates (NSC) in woody plants. Across species, the amount of parenchyma depends on the general environment of the distribution area and determines to a large extent the NSC storage. However, little information is available on the relationship between parenchyma fractions, NSC storage, and the environmental inﬂuences within individual species. This information is crucial to assessing the adaptive capacities of tree populations in the context of increasing the frequency and severity of stress-inducing events. In this study, parenchyma fractions and NSC concentrations of the secondary xylem in trunks of a subtropical evergreen oak ( Quercus ciliaris C.C.Huang & Y.T.Chang) were quantiﬁed along an elevational gradient from 700 m to 1200 m a.s.l. in eastern China. Air temperatures within the distribution area correlated with altitude were recorded. The results showed that the total parenchyma fractions did not covary with the colder temperatures. However, axial parenchyma fractions were lower with a colder climate, while the fractions of multiseriate rays and total ray parenchyma were higher. Higher concentrations of starch and NSC were signiﬁcantly associated with larger axial parenchyma fractions. The sugar concentration displayed no signiﬁcant relationship with parenchyma fractions. These ﬁndings suggest that the total parenchyma fractions in secondary xylem do not increase in response to a colder climate, while colder temperatures drive changes in the composition of parenchyma for Q. ciliaris .


Introduction
It is well known that the majority of a tree's biomass is contained within the secondary xylem of the stems [1][2][3]. These woody tissues are generally composed of vessels (and, in some cases, tracheids), fibers, and parenchyma. Among these, parenchyma is linked with multiple plant functions, such as the transport and storage of non-structural carbohydrates (NSC, including starch and soluble sugar) and water, vessel embolism repair, among others [4]. Parenchyma can be further classified as either axial or ray parenchyma [5]. Although the relationship between the properties and functions of woody tissue components, such as vessels and fibers are well studied [6][7][8][9], the relationships between the properties and functions of parenchyma are much less well understood [10].
The patterns associated with variations in parenchyma fractions give insight into the functional consequences and ecological significance of low vs. high fractions. For This study was conducted in Tianma National Nature Reserve, located at the junction of Anhui, Henan, and Hubei provinces (115 • 20 -115 • 50 E, 30 • 10 -31 • 20 N), in eastern China. This region is characterized by a subtropical moist monsoon climate. The mean annual temperature (MAT) is 13.8 • C, with an extreme maximum temperature of 38.1 • C (July) and a minimum of −23 • C (January). The mean annual precipitation is 1489 mm, most of which occurs from May to September. This region consists of mixed evergreen and deciduous broad-leaved forests [39] with a transition zone buffering the two forest types [37,40].
In terms of the studied elevational gradient, the MAT varied from 10.7 • C to 12.8 • C during the year of sampling, with the coldest month mean temperature (CMMT) varying from −0.3 • C to −2.9 • C and the warmest month mean temperature (WMMT) varying from 21.3 • C to 23 that affects the distribution of Q. ciliaris. In this study, the fractions of the secondary xylem parenchyma and NSC concentrations in the trunk of Q. ciliaris along an elevational gradient were measured to examine the hypotheses that within the single species, (1) parenchyma fractions would decrease towards colder temperatures, and (2) parenchyma fractions would be positively correlated with NSC concentrations.

Study Area
This study was conducted in Tianma National Nature Reserve, located at the junction of Anhui, Henan, and Hubei provinces (115°20′-115°50′ E, 30°10′-31°20′ N), in eastern China. This region is characterized by a subtropical moist monsoon climate. The mean annual temperature (MAT) is 13.8 °C, with an extreme maximum temperature of 38.1 °C (July) and a minimum of −23 °C (January). The mean annual precipitation is 1489 mm, most of which occurs from May to September. This region consists of mixed evergreen and deciduous broad-leaved forests [39] with a transition zone buffering the two forest types [37,40].
In terms of the studied elevational gradient, the MAT varied from 10.7 °C to 12.8 °C during the year of sampling, with the coldest month mean temperature (CMMT) varying from −0.3 °C to −2.9 °C and the warmest month mean temperature (WMMT) varying from 21.3 °C to 23.8 °C (Figure 1). The upper elevational limit (ca. 1200 m a.s.l.) of Q. ciliaris in the study area corresponded to the CMMT of −2.9 °C, the WMMT of 21.3 °C, and the MAT of 10.7 °C (Figure 1). Figure 1. The relationships of temperature variables with elevation. The significant Pearson correlation coefficients were presented (p < 0.05). MAT, mean annual temperature; CMMT, the coldest month mean temperature; WMMT, the warmest month mean temperature.

Sample Collection and Determination
Quercus ciliaris, a dominant evergreen broad-leaved canopy tree species, reaches both the northern [36] and upper elevational distributional limits of subtropical evergreen broad-leaved forests [37,41,42] in the studied area. Tree core samples were collected with a 5.15 mm increment borer (Haglöf, Långsele, Sweden) from three healthy adults at 740 m, 837 m, 922 m, 1016 m, 1120 m, 1148 m, and 1200 m of elevation in October 2018 when no frost nights occurred. According to our monitoring temperature data (see Section 2.3), the mean temperature in October 2017 varied from 10.7 °C to 12.7 °C within the studied elevational gradient. Samples were collected from approximately 1 m above the ground. A total of 21 individuals were sampled. Two cores were collected from each tree, one of which was treated in a microwave oven within 12 h of sampling in order to denature enzymes [43]. Next, the samples were dried to constant mass at 65 °C. The outer 5 cm of secondary xylem was used for NSC measurement following an improved phenol-concentrated sulfuric acid method [44].
The second core from each tree was stored in FAA reagent (formalin: glacial acetic acid: ethanol) [45] until it was sectioned for anatomical observations. As the wood of Q. The relationships of temperature variables with elevation. The significant Pearson correlation coefficients were presented (p < 0.05). MAT, mean annual temperature; CMMT, the coldest month mean temperature; WMMT, the warmest month mean temperature.

Sample Collection and Determination
Quercus ciliaris, a dominant evergreen broad-leaved canopy tree species, reaches both the northern [36] and upper elevational distributional limits of subtropical evergreen broad-leaved forests [37,41,42] in the studied area. Tree core samples were collected with a 5.15 mm increment borer (Haglöf, Långsele, Sweden) from three healthy adults at 740 m, 837 m, 922 m, 1016 m, 1120 m, 1148 m, and 1200 m of elevation in October 2018 when no frost nights occurred. According to our monitoring temperature data (see Section 2.3), the mean temperature in October 2017 varied from 10.7 • C to 12.7 • C within the studied elevational gradient. Samples were collected from approximately 1 m above the ground. A total of 21 individuals were sampled. Two cores were collected from each tree, one of which was treated in a microwave oven within 12 h of sampling in order to denature enzymes [43]. Next, the samples were dried to constant mass at 65 • C. The outer 5 cm of secondary xylem was used for NSC measurement following an improved phenol-concentrated sulfuric acid method [44].
The second core from each tree was stored in FAA reagent (formalin: glacial acetic acid: ethanol) [45] until it was sectioned for anatomical observations. As the wood of Q. ciliaris is very hard, the cores were boiled for 2 h to soften the wood for sectioning. Next, transverse sections of entire tree cores approximately 20 µm thick were prepared using a microtome (WSL core-microtome, WSL, Birmensdorf, Switzerland) as described in Gärtner et al. 2015 [46]. The sections were rehydrated in distilled water and then immersed in a mixture of 0.35% w/v safranin (in 50% ethanol) and 0.65% w/v alcian blue (in distilled water) for 2 min, thoroughly washed with distilled water, and finally they were gradually dehydrated through an ethanol series (55%, 75%, 95%, 100%) before being mounted on a slide in Neo-Mount (Merck Millipore, Darmstadt, Germany) [31] (Figure 2a). One tangential section ( Figure 2b) and one radial section ( Figure 2c) were also prepared to help the analysis of the transverse sections [47]. Sequences of digital images were captured using a microscope (Olympus SZ61, Olympus Corporation, Tokyo, Japan), which were then stitched into a composite image of each full tree core using PTGui Pro (version 11.12, New House Internet Services BV, Rotterdam, The Netherlands). The fractions of axial and ray parenchyma were measured in the transverse sections. Since the biseriate or triseriate rays (two or three cells in width) were hardly visible, ray parenchyma was classified into the following two categories: thin rays, including uniseriate and, if any, biseriate or triseriate rays, and multiseriate rays (more than three cells wide). ciliaris is very hard, the cores were boiled for 2 h to soften the wood for sectioning. Next, transverse sections of entire tree cores approximately 20 μm thick were prepared using a microtome (WSL core-microtome, WSL, Birmensdorf, Switzerland) as described in Gärtner et al. 2015 [46]. The sections were rehydrated in distilled water and then immersed in a mixture of 0.35% w/v safranin (in 50% ethanol) and 0.65% w/v alcian blue (in distilled water) for 2 min, thoroughly washed with distilled water, and finally they were gradually dehydrated through an ethanol series (55%, 75%, 95%, 100%) before being mounted on a slide in Neo-Mount (Merck Millipore, Darmstadt, Germany) [31] (Figure 2a). One tangential section ( Figure 2b) and one radial section ( Figure 2c) were also prepared to help the analysis of the transverse sections [47]. Sequences of digital images were captured using a microscope (Olympus SZ61, Olympus Corporation, Tokyo, Japan), which were then stitched into a composite image of each full tree core using PTGui Pro (version 11.12, New House Internet Services BV, Rotterdam, The Netherlands). The fractions of axial and ray parenchyma were measured in the transverse sections. Since the biseriate or triseriate rays (two or three cells in width) were hardly visible, ray parenchyma was classified into the following two categories: thin rays, including uniseriate and, if any, biseriate or triseriate rays, and multiseriate rays (more than three cells wide).  As the parenchyma fraction of the trunk is known to vary radially [48], the parenchyma fraction for each individual is an averaged fraction calculated over 2 mm × 1 mm rectangular measurements taken from the outside, middle, and inside of the outer 5 cm of the composite core image. For a given rectangular measurement area, the fraction of parenchyma was defined as the area of parenchyma divided by the total area. Proportions of axial and ray parenchyma were quantified by manually outlining their areas in different colors in Photoshop software (version 19.1.9, Adobe Systems Inc., San Jose, CA, USA) and then measuring them in ImageJ (version 1.50b, National Institutes of Health, Bethesda, MD, USA). The scatter plots of parenchyma and NSC variables vs. elevation were shown in Figure 3.
As the parenchyma fraction of the trunk is known to vary radially [48], the parenchyma fraction for each individual is an averaged fraction calculated over 2 mm × 1 mm rectangular measurements taken from the outside, middle, and inside of the outer 5 cm of the composite core image. For a given rectangular measurement area, the fraction of parenchyma was defined as the area of parenchyma divided by the total area. Proportions of axial and ray parenchyma were quantified by manually outlining their areas in different colors in Photoshop software (version 19.1.9, Adobe Systems Inc., San Jose, CA, USA) and then measuring them in ImageJ (version 1.50b, National Institutes of Health, Bethesda, MD, USA). The scatter plots of parenchyma and NSC variables vs. elevation were shown in Figure 3. . Scatter plots of parenchyma fraction and non-structural carbohydrate concentration (yaxis) vs. Elevation (x-axis). The significant Pearson correlation coefficients were presented (p < 0.05). AP, axial parenchyma; RP, ray parenchyma (multiseriate rays + thin rays); RAP, AP + RP; NSC, total non-structural carbohydrate (soluble sugar + starch).
To visualize the distribution of starch within the samples, transverse sections were stained with Lugol's solution for 2 min and then thoroughly washed in distilled water for immediate observation with a microscope (Olympus SZ61, Olympus Corporation, Tokyo, Japan) ( Figure 4). Elevation (x-axis). The significant Pearson correlation coefficients were presented (p < 0.05). AP, axial parenchyma; RP, ray parenchyma (multiseriate rays + thin rays); RAP, AP + RP; NSC, total non-structural carbohydrate (soluble sugar + starch).
To visualize the distribution of starch within the samples, transverse sections were stained with Lugol's solution for 2 min and then thoroughly washed in distilled water for immediate observation with a microscope (Olympus SZ61, Olympus Corporation, Tokyo, Japan) ( Figure 4).  The detailed procedure for measuring NSC concentrations is described as because NSC concentration of the trunk is known to vary radially [49], all the ou of the dried tree cores were crushed into a fine powder. The soluble sugars were e by weighing about 40 mg of powder into a centrifuge tube and recording the weig ing 10 mL of 80% ethanol into tubes and incubating the closed tubes in a water ba min at 80 °C followed by centrifuging at 4000 r/min (2200 × g) for 15 min. The sup was transferred to 100 mL volumetric flasks. 5 mL of 80% ethanol was then add tubes. The processes of water bath, centrifuging, and supernatant collection w repeated. The total volume of the combined supernatant was adjusted to 100 mL measurement.
Starch remaining in the undissolved pellet after drying the residuals in the tu hydrolyzed in a boiling 3% HCl solution for 3 h. Moreover, after filtering the solu a new batch of 100 mL volumetric flasks, the total volume of the solution was ad 100 mL to be measured later.
In a glass tube, 1 mL of the extracting solution and 1 mL 28% phenol solut phenol adding 80% ethanol to 300 g) were added. Immediately afterwards, 5 mL trated sulfuric acid was added. The tube was shaken for 1 min and rested for a min before absorbance was measured at 490 nm (UV-9600). The standard solu sucrose. The sum of soluble sugar and starch concentration was referred to as "t structural carbohydrates" (NSC).

Temperature Measurement
The air temperature of each sampling site was recorded hourly for one year 2017-September 2018) using a HOBO Pro data logger (U23-001 Pro v2, Onset C Corporation, Bourne, MA) mounted approximately 1.5 m above the ground. As nary observations on the population density of Q. ciliaris in the study area show decrease in population density approaching the elevational distributional li mainly driven by colder climate [38], as well as the upper distributional limits of e broad-leaved forests in East Asia correspond with mean temperature of the colde [50], the temperature time series data were summarized to calculate the coldes mean temperature (CMMT) as the temperature index (Figure 1).

Data Analysis
Since the residuals of linear mixed models did not comply with normal dist The detailed procedure for measuring NSC concentrations is described as follows: because NSC concentration of the trunk is known to vary radially [49], all the outer 5 cm of the dried tree cores were crushed into a fine powder. The soluble sugars were extracted by weighing about 40 mg of powder into a centrifuge tube and recording the weight, adding 10 mL of 80% ethanol into tubes and incubating the closed tubes in a water bath for 60 min at 80 • C followed by centrifuging at 4000 r/min (2200× g) for 15 min. The supernatant was transferred to 100 mL volumetric flasks. 5 mL of 80% ethanol was then added to the tubes. The processes of water bath, centrifuging, and supernatant collection were then repeated. The total volume of the combined supernatant was adjusted to 100 mL for later measurement.
Starch remaining in the undissolved pellet after drying the residuals in the tubes was hydrolyzed in a boiling 3% HCl solution for 3 h. Moreover, after filtering the solution into a new batch of 100 mL volumetric flasks, the total volume of the solution was adjusted to 100 mL to be measured later.
In a glass tube, 1 mL of the extracting solution and 1 mL 28% phenol solution (84 g phenol adding 80% ethanol to 300 g) were added. Immediately afterwards, 5 mL concentrated sulfuric acid was added. The tube was shaken for 1 min and rested for at least 15 min before absorbance was measured at 490 nm (UV-9600). The standard solution was sucrose. The sum of soluble sugar and starch concentration was referred to as "total non-structural carbohydrates" (NSC).

Temperature Measurement
The air temperature of each sampling site was recorded hourly for one year (October 2017-September 2018) using a HOBO Pro data logger (U23-001 Pro v2, Onset Computer Corporation, Bourne, MA) mounted approximately 1.5 m above the ground. As preliminary observations on the population density of Q. ciliaris in the study area showed that a decrease in population density approaching the elevational distributional limit was mainly driven by colder climate [38], as well as the upper distributional limits of evergreen broadleaved forests in East Asia correspond with mean temperature of the coldest month [50], the temperature time series data were summarized to calculate the coldest month mean temperature (CMMT) as the temperature index ( Figure 1).

Data Analysis
Since the residuals of linear mixed models did not comply with normal distribution, we employed generalized linear mixed models utilizing a beta error structure and a logit link function to test the effects of temperature on the fractions of parenchyma using temperature as a fixed effect and elevation as a random effect. Elevation was included as a random effect in each linear mixed model since individuals occurring at the same elevation may be more similar and could therefore violate the assumption of independence. Ordinary linear regression was used to test the influence of parenchyma fraction on NSC. To explore the multiple correlation patterns in parenchyma and NSC, we also conducted a principal component analysis (PCA). All statistical analyses were implemented in R software (version 4.0.2) [51]. Generalized linear mixed models were written using the glmmTMB function in the glmmTMB package (version 1.1.2.3) [52]. The simulate Residuals function in the DHARMa package [53] was used for residual diagnostics of the generalized linear mixed models. Ordinary linear regression was performed with the lm function in the base package [51]. PCA was performed using the rda function in the vegan package (version 2.5.6) [54].

Variations in Parenchyma Fractions along a Temperature Gradient
The axial parenchyma fractions decreased significantly with a colder CMMT (Table 1). Meanwhile, the fractions of multiseriate ray parenchyma and total ray parenchyma increased significantly with a colder CMMT ( Table 1). The fractions of thin ray parenchyma and total parenchyma displayed a non-significant linear relationship with CMMT (Table 1).

Relationships between Parenchyma Fractions and Carbohydrate Storage
The results of the linear models indicated that NSC concentrations were not significantly associated with the total parenchyma fraction. However, a significantly positive correlation was identified between NSC and the axial parenchyma fraction, along with a significantly negative linear relationship between the fraction of ray parenchyma or multiseriate ray and NSC concentration ( Figure 5). Neither the fraction of any kind of ray parenchyma nor the total parenchyma exhibited linear associations with starch concentration. However, starch concentration was positively associated with the axial parenchyma fraction ( Figure 5). No significant linear relationships were identified for soluble sugar concentration and any kind of parenchyma ( Figure 5).
The results from the PCA show that the first and second axes accounted for 72.20% of the total variation ( Figure 6). The PCA confirmed that fractions of axial parenchyma, starch concentration, and NSC concentration displayed positive correlations with each other. The PCA showed that NSC was also negatively correlated with multiseriate rays. Fractions of total parenchyma were positively associated with ray parenchyma and its two components (multiseriate and thin ray). However, no correlation was identified between the fraction of ray parenchyma and the fraction of axial parenchyma. x FOR PEER REVIEW 8 of 14 Figure 5. Linear models test the effects of fraction of axial parenchyma, ray parenchyma, multiseriate ray parenchyma, thin ray parenchyma, and total parenchyma on the concentration of soluble sugar, starch, and total non-structural carbohydrates (NSC). The blue lines indicate that the relationships were significant (p < 0.05).
The results from the PCA show that the first and second axes accounted for 72.20% of the total variation ( Figure 6). The PCA confirmed that fractions of axial parenchyma, starch concentration, and NSC concentration displayed positive correlations with each other. The PCA showed that NSC was also negatively correlated with multiseriate rays. Fractions of total parenchyma were positively associated with ray parenchyma and its two components (multiseriate and thin ray). However, no correlation was identified between the fraction of ray parenchyma and the fraction of axial parenchyma.  Figure 5. Linear models test the effects of fraction of axial parenchyma, ray parenchyma, multiseriate ray parenchyma, thin ray parenchyma, and total parenchyma on the concentration of soluble sugar, starch, and total non-structural carbohydrates (NSC). The blue lines indicate that the relationships were significant (p < 0.05).

Figure 5.
Linear models test the effects of fraction of axial parenchyma, ray parenchyma, multiseriate ray parenchyma, thin ray parenchyma, and total parenchyma on the concentration of soluble sugar, starch, and total non-structural carbohydrates (NSC). The blue lines indicate that the relationships were significant (p < 0.05).
The results from the PCA show that the first and second axes accounted for 72.20% of the total variation ( Figure 6). The PCA confirmed that fractions of axial parenchyma, starch concentration, and NSC concentration displayed positive correlations with each other. The PCA showed that NSC was also negatively correlated with multiseriate rays. Fractions of total parenchyma were positively associated with ray parenchyma and its two components (multiseriate and thin ray). However, no correlation was identified between the fraction of ray parenchyma and the fraction of axial parenchyma.

Discussion
Currently, scientific literature investigating the variations in parenchyma fractions along environmental gradients at the intraspecific level is lacking, although several studies have conducted analyses across species [5,20,21]. In this study, a species of evergreen oak (Q. ciliaris) found in the northern subtropics of China was employed to investigate the variations in the parenchyma fractions along a temperature gradient. Our main findings show no changes in total parenchyma fractions with a colder climate, suggesting that the xylem storage capacity may not increase with the increasingly stressful conditions within the single species. A possible explanation for this may be that maintaining more living cells, as well as protecting them from frost damage, can be an energy-demanding process [5,22,23]. However, the trends of axial and ray parenchyma diverged with a colder climate, suggesting a functional differentiation of axial parenchyma and ray parenchyma. There is much indirect or direct evidence to support the suggestion in previous studies. Axial and ray parenchyma diverged, with climate having little influence over the amount of ray parenchyma, unlike axial parenchyma with regards to general trends over more than a thousand species [5]. Axial and ray parenchyma are produced by the fusiform and ray initials of the vascular cambium, respectively, and run perpendicular to each other [5,16]. The axially oriented axial parenchyma tends to perform the hydraulic function, since positive correlations between hydraulic efficiency and axial parenchyma fraction have been repeatedly found [17,18]. Ray parenchyma is oriented radially and thus can perform the functions of radial transfer [55] and mechanical support [4,17]. Unlike the axial parenchyma, the ray parenchyma fraction is negatively related to hydraulic efficiency [10,17,18].
Axial parenchyma fractions decreased with colder temperatures for Q. ciliaris adult trees (Table 1), which is consistent with the decreasing trend in axial parenchyma across angiosperm species on a global scale [5], and a tendency toward less abundant axial parenchyma at higher latitudes within the Ilex genus [56], as well as a tendency toward less abundant axial parenchyma in trees at higher latitudes in Brazil [57]. However, fractions of axial parenchyma for Buddleja cordata Kunth did not covary with latitude or elevation in Mexico [26], but increased with lower temperatures for Paubrasilia echinata (Lam.) Gagnon, H.C.Lima & G.P.Lewis in Brazil [24]. Therefore, while most studies confirmed that lower fractions of axial parenchyma would likely be observed at sites with low temperatures or high latitudes across species, there are not enough studies to establish the general trends within species. To our knowledge, the formation of air embolisms in vessels due to winter freezing is critical for the distribution of evergreen broad-leaved trees [58]. While all the axial parenchyma of Q. ciliaris were apotracheal parenchyma (Figure 2), each vessel of Quercus is sheathed by many layers of vasicentric tracheids, which serve as subsidiary conductive cells [16]. The abundant vasicentric tracheids were intermixed with the parenchyma (Figure 2), forming the connective tissue between vessels [59]. Quercus that evolved in climates with drier, hotter summers or colder winters had higher numbers of vasicentric tracheids than those from mesic climates, supporting the hypotheses that vasicentric tracheids are critical for water transport during drought [60]. Therefore, compared to a high total apotracheal parenchyma fraction, the abundant vasicentric tracheids intermixed with a small amount of parenchyma of Q. ciliaris ( Figure 2) may be more critical for stabilizing hydraulic capacity during freezing winters. Consequently, the less axial parenchyma fraction may not have much influence on hydraulic capacity during winter. Additionally, smaller vessels were linked with higher resistance to freezing-induced cavitation [6]. In our study, the decreasing fraction of axial parenchyma may also result from the decreased diameter of vessels with increasing elevation [61,62], since vessel diameter can be positively associated with axial parenchyma fraction, whether across species [24,63] or within species [24]. These trends in vessel properties of Q. ciliaris with temperature are beyond the scope of this paper and warrant further investigation.
Although ray parenchyma fractions in woody angiosperm species have been shown to decrease with colder temperatures or towards higher latitudes across species [20,21], our study showed an increasing trend in the fraction of ray parenchyma with colder temperatures within the single species (Table 1). A previous study with a focus on Juniperus thurifera (an evergreen conifer) indicated that the mean temperature in the previous October and in the current August was negatively associated with the fraction of ray parenchyma after controlling for ring-width [64]. There was also a similar result for ray parenchyma fraction in evergreen conifers, where the fraction tended to increase with colder temperatures, while the ray parenchyma fraction in an evergreen oak (Quercus chrysolepis) did not significantly vary with temperature within the species [21]. In our study, the increasing trend of ray parenchyma fraction was mainly driven by the multiseriate ray fraction (Table 1). It was hypothesized that wide and tall rays were able to adapt better to seasonal changes in water storage than narrow and short rays [16]. Consequently, we speculated that the increasing multiseriate ray parenchyma fraction towards more freezing winter climates may help repair vessel embolism by storing more water or radially transferring water from phloem to xylem.
In our study, total NSC concentration was positively correlated with fractions of axial parenchyma in trunks of Q. ciliaris, which was largely due to the positive correlation between starch concentration and axial parenchyma fraction (Figures 5 and 6). The axial parenchyma fraction may be a proxy for the concentration of NSC (especially for starch) for Q. ciliaris. Positive correlations between starch or total NSC concentration and axial parenchyma fraction were also observed in previous inter-species studies with a much higher R 2 (≥0.5) [10,65] than in our study (0.21 or 0.30 in our study). The low sample size in the study may partly explain the lower R 2 . While NSC concentration showed a weak positive correlation with the ray parenchyma fraction across species for three oaks [21], a significantly negative correlation was observed in our study for a single species. This negative trend was mainly caused by the opposing trends of axial parenchyma and ray parenchyma with temperature.
While most evidence to date has demonstrated a positive correlation between the parenchyma fraction and the concentration of starch or NSC, little evidence has been collected for the correlation between the parenchyma fraction and the soluble sugar concentration. While soluble sugar reaches its highest level in the dormant season [66], all previous studies have been sampled during the growing season. Therefore, a lack of observed correlation may be due to inappropriate sampling design, and tests occurring in the dormant season are therefore recommended. There is a known dynamic of NSC within the dormant season [67,68]. Starch is being mobilized long before bud burst [68,69]. The timing of sampling for starch/sugar analysis is crucial, and the choice of the right moment should be underpinned by detailed data on the periodicity of NSC dynamics. Therefore, for further study, it is better to carry out sampling monthly with large sample sizes in different phenological periods to verify the detailed dynamics of the movement and concentration of NSC during the whole year. Combined with visualizing starch grains in parenchyma at different phenological stages, the roles of different parenchyma in carbon storage can be more clearly distinguished.
In addition, previous studies have also shown that tree heights are associated with wood anatomy and hydraulic functions. For example, taller trees have predictably wider vessels [70]. Additionally, wood anatomy and hydraulic functions changed when measured at different points from the stem tip towards the base [71]. That means standardizing measurements by distance between the stem tip and tree height should improve comparability between different individuals [71]. Even though all the trees sampled in the current study were mature individuals, the effects of tree size were not taken into account. In further studies, standardizing measurements by distance between the stem tip and tree height is recommended to enhance comparability.

Conclusions
In this study, we employed a species of evergreen oak to investigate the variations in parenchyma fraction along a temperature gradient, along with the relationships between the parenchyma fraction and the NSC concentration at the end of a growing season in a northern subtropical region. The fraction of total parenchyma did not increase with colder temperatures. However, the axial parenchyma and ray parenchyma fractions changed in different trends with proximity to the elevational distribution limit of Q. ciliaris. Axial parenchyma fractions were lower for individuals in colder climates, while the opposite trend was true for ray parenchyma. The increasing trend of ray parenchyma fraction mainly resulted from the variations in the multiseriate ray fractions. As part of the structural basis of xylem storage, the axial parenchyma fraction was a proxy for the concentration of NSC, especially for starch concentration. While many studies investigating the links between the parenchyma fractions and their functions have been performed across species, such investigations have been rare within species. Therefore, to improve our understanding of the physiological significance of variations in parenchyma fraction, similar studies across and within species are needed.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.