Main drivers of vertical and seasonal patterns of leaf photosynthetic characteristics of young planted Larix Olgensis trees

Photosynthetic characteristics are indicative of physiological parameters that correlate with forest primary production^[1] and drive carbon uptake^[2]. They directly participate in the natural carbon cycle but are sensitive to environmental conditions; they are often used to represent the resistance and resilience of vegetation to extreme climates^[3], pests and diseases^[4,5], fires^[6], and toxic metals^[7,8]. Convergent emergence or loss of photosynthetic phenotypes may facilitate adaptation to ecologically similar environments^[9].

In forestry, differences in photosynthetic characteristics among different tree species are prominent^[10]. The photosynthetic characteristics of leaves within the same species often have significant spatial heterogeneity owing to the complex spatial structure of the crown^[11], which is particularly evident in the vertical structure of the crown^[12,13]. Light^[14,15] and water potential^[16] are thought to be the key determining factors of vertical variation patterns in tree crown photosynthetic traits. However, the dominant roles of two species could shift under different forest densities and tree sizes. Generally, in closed canopies, light is the major factor that leads to the spatial heterogeneity of leaf photosynthetic traits^[17], causing leaves to alter their structure and physiological function to adapt to the lighting environment. Likewise differences appear between shade and sun leaves within the crown^[13]. Sun leaves have a high maximum net CO₂ assimilation (A_max), respiration rate (R_D), and light compensation point (LCP) but lower light utilization efficiency (LUE) and light saturation point (LSP)^[18]. Water potential tends to be the dominant factor for huge dominant trees^[19], as it represents the ability to transport the water from the root to the crown. Crown photosynthetic characteristics also show seasonal variation with individual development and environmental changes^[20,21]. For example, A_max usually shows a parabolic seasonal variation pattern due to its positive correlation with the temperature (T), solar radiation, and soil moisture^[22−24]. However, R_D exhibits an opposite 'U'-shaped trend^[18] due to the decline of cytochrome-mediated respiration^[21] and temperature sensitivity^[25].

Studies have shown that spatial and seasonal variations in canopy photosynthetic characteristics are closely associated with the comprehensive impact of light^[26], temperature^[24,27], humidity^[28], and seasonal patterns of leaf structural traits^[29,30]. Generally, adequate light and suitable temperature and humidity can promote photosynthesis, by improving light energy utilization^[31] and photosynthetic enzyme activity^[32], however, some environmental conditions will inhibit photosynthesis^[33]. For example, the common natural phenomena 'midday depressions' is a self protection mechanism by regulating blade osmotic pressure and maintaining mesophyll cell activity, in response to the stress of strong light, high temperature and low humidity^[34]. Additionally, the age^[22,35,36] and sex of dioecious tree species^[37] influence the photosynthetic characteristics of leaves to some extent due to the difference of mesophyll cell structure. When simulating the mechanism process model of canopy, forest productivity, and carbon absorption processes, the spatial and seasonal variations in canopy leaf photosynthetic characteristics need to be considered simultaneously, otherwise, incorrect results will be obtained^[16]. Therefore, the spatiotemporal heterogeneity and driving mechanisms of canopy photosynthetic characteristics need to be urgently clarified.

This study used the artificial forest of the main afforestation species (Larix olgensis) in the northern area of China as the research object that was tracked and monitored throughout its growing season. The factors of photosynthetic characteristics, leaf morphology, crown structure and environmental factor were collected. The overarching aims of this study are twofold: first, is there a significant difference in photosynthetic characteristics of different crown positions, growing periods and tree individuals? Secondly, what is the variation patterns of photosynthetic characteristics within the tree crown during the growing season? Finally, we provide a comprehensive assessment of the relationships between photosynthetic characteristics and crown structure, leaf morphology and environmental conditions.

The experiments were conducted in 2017 at the experimental forest farm of Northeast Forestry University in Maoershan, Haerbin, China (Northern latitude: 45°2′20″~45°18′16″, East longitude: 127°18′0″~127°41′6″; altitude 400 m above sea level). The climate in the Maer Mountain region belongs to the temperate continental monsoon climate, with an average annual temperature of 2.4 °C, the highest temperature being 34 °C and the lowest being −40 °C, approximately 125 d of frost-free period, an average annual precipitation of 700 mm, and dark brown soil as the main type. Total forest coverage is approximately 83.3%, including 14.7% plantation.

Sample selection

Five fixed plots of 20 m × 30 m with same site quality in the young L. olgensis plantation are set up, and the trees with a diameter at breast height (DBH) larger than 5 cm in each plot were measured. The specific measurement factors include tree height (H), DBH, crown width (CW) and the relative coordinates (x_i, y_i) of each tree were investigated. Then the average DBH of five plots were calculated according to the per tree measurement data. A scaffold was built around the sample tree, ensuring that the sample tree was completely surrounded by it, ensuring that all branches of the sample wood and each position of each branch can be measured on the scaffold, each layer is connected by a walkway. After each measurement, the upper tread was removed to avoid the influence on the measurement result caused by the blocking of light. The crown length of each sample tree was divided into several vertical sections based on the whorls from tree top to bottom and numbering began from V1st. Three healthy and fully expanded needles located within each section in the middle of the foliated branches in sunny, semisunny and shaded crowns were selected, according to the sample selection principles based on Liu et al.^[38]. The photosynthetic characteristics in each vertical section was the average values of measurements taken from different directions (sunny, semisunny and shaded crown).

Photosynthetic gas exchange measurements

The photosynthetic light response (PLR) curves were measured twice per month (the beginning of the month and mid-month) during the growing season (from approximately May 15^th to September 10^th). All photosynthetic properties were measured with a portable steady-state photosynthesis system (LI-6400XT, LI-COR, Inc., Lincoln, NE, USA) equipped with a standard LED light source (6400-02B, LI-COR, Inc., Lincoln, NE, USA). Sample chamber was acclimated for 20 min at a CO₂ concentration of 390 ppm with a CO₂ mixer (6400-01, LI-COR, Inc., Logan, NE, USA) to maintain a stable CO₂ supply. All sample cluster needles were acclimated under a PAR of 1,400 μmol·m⁻²·s⁻¹ for 10 to 20 min by the LED light source (6400-02B, LI-COR, Inc., Lincoln, NE, USA). PLR curve was measured at 10 PAR gradients: 2,000, 1,500, 1,200, 1,000, 500, 200, 150, 100, 50, 0 μmol·m⁻²·s⁻¹. Sample cluster needles were allowed to equilibrate for a minimum of 2 min at each measurement before data was logged, and a calibration (match) was performed after each count. At the same time, the temperature of the leaf (Air temperature, Tair), the relative humidity (Relative humidity, RH) and the vapor pressure deficit (Vapor pressure deficit, VPD) were recorded. After the measurements, the depth into the crown (Depth into crown, DINC) of each measured sample was recorded in the crown, and the relative depth into the crown (Relative depth into crown, RDINC) was calculated according to the crown length (Crown length, CL): RDINC = DINC/CL.

Needle morphology measurements

Once the photosynthetic gas exchange measurements were completed, the sample cluster needles were immediately taken back to the laboratory for measuring the needle mass per area (NMA, g·m⁻²). Each cluster sample was scanned immediately after collection and then surveyed with an image analysis software (Image-Pro Plus 6.0, Media Cybernetics, Inc., Bethesda, USA), resulting in the projected needle area (NA, m²), needle length (l) and needle width (w), and consequently obtained the ratio of needle length to needle width (r_lw). Then, the scanned samples were dried to a constant weight at 85 °C and weighed to dry weight (WD). The NMA was calculated: NMA = WD/LA.

Photosynthetic parameters

The light-saturated CO₂ assimilation (A_max, μmol·m⁻²·s⁻¹) and dark respiration (R_D, μmol·m⁻²·s⁻¹) were estimated from the PLR curves using the modified Mitscherlich model^[39]:

where A_n is the net CO₂ assimilation (μmol·m⁻²·s⁻¹), A_max is the light-saturated net CO₂ assimilation (μmol·m⁻²·s⁻¹), α is the apparent quantum yield, PAR is the photosynthetically active radiation (μmol·m⁻²·s⁻¹), and R_D is the dark respiration rate (μmol·m⁻²·s⁻¹).

The light saturated stomatal conductance (g_s-sat, mol·m⁻²·s⁻¹) was determined as the corresponding gs value of A_max. Water-use efficiency (WUE_sat, mmol CO₂ mol H₂O⁻¹) was calculated as the ratio of A_max to g_s-sat. As the environment conditions were not maintained under a certain value during the measurement of PLR curves except CO₂ concentration (stabilized at 390 ppm).

Statistical analysis

Table 1 shows the data summary. Statistical analyses were performed using R software 4.2.2^[40]. A three-way repeated-measures analysis of variance (ANOVA) was performed on all experimental variables to evaluate the effects of individual tree (T), period (P), and crown layer (L) on light-saturated CO₂ assimilation (A_max), dark respiration (R_D), light-saturated stomatal conductance (g_s-sat), and water-use efficiency (WUE_sat). Pearson's correlation analysis was used to test the relationships among all the measured variables. The significance of all the statistical analyses was at α = 0.05 level. All figures were drawn using the ggplot2 package in R software 4.2.2.

Photosynthetic parameters (A_max, R_D, g_s-sat and WUE_sat) and morphological parameters (LMA and r_lw) differed significantly among the different measurement phases, individual trees, and vertical locations of the crown (Table 2). Considering the average pattern across the five sample trees, nearly all physiological and morphological parameters of the needles exhibited a similar vertical profile, which decreased noticeably with increasing RDINC (Fig. 1). However, the mean WUE_sat of the five sampled trees followed the opposite trend (Fig. 1f). A_max significantly decreased with RDINC (Fig. 1a), and the mean A_max in the top crown (12.84 μmol·m⁻²·s⁻¹) was almost 2.7 times higher than that in the bottom crown (4.81 μmol·m⁻²·s⁻¹). Although there were significant tree-specific differences in g_s-sat, NMA, and R_D (Table 2), their tendencies demonstrated a pronounced decrease with increasing RDIINC (Fig. 1b, c, & d). The mean values of R_D, g_s-sat, and NMA varied by 2.5-fold, 3.3-fold, and 2.3-fold, respectively, from the top to the bottom of the crown. Mean r_lw exhibited a slight decrease with RDINC when RDINC was lower than 0.3 (upper crown) (Fig. 1e) but then sharply decreased when RDINC was greater than 0.4 (middle and lower crown). In contrast, mean WUE_sat showed a weak upward trend with increasing RDINC (Fig. 1f), varying by only 0.015 mmol CO₂ mol H₂O⁻¹ from top to bottom.

Figure 1. 

Vertical profiles of (a) light-saturated net photosynthetic rate (A_max); (b) light-saturated stomatal conductance (g_s-sat); (c) needle mass per area (NMA); (d) dark respiration (R_D); (e) ratio of length to width of needles (r_lw) and (f) light-saturated water use efficiency (WUE_sat). Data points represented seasonal mean values (solid bars represented stand error). Black solid line represented mean values of five sample trees (a)−(f).

As mentioned above, photosynthetic and morphological parameters were significantly affected by the vertical location of the crown. However, it is unknown whether the same pattern remains during the entire growth period. Analysis of variance was performed on photosynthetic and morphological parameters based on the vertical locations (upper, middle, and lower crown) in each measurement phase, and the results are summarized in Table 3. r_lw was the only parameter that showed a significant vertical difference across the entire growth season (upper > middle > lower crown). A_max, g_s-sat, and NMA demonstrated a similar vertical pattern in June, but the mean values of A_max and g_s-sat were not significantly different between the upper and middle crown during the early period of needle expansion (PI, May). Mean NMA showed no significant vertical difference within the crown. R_D showed no significant difference with respect to vertical location in the crown during the early period of needle expansion but showed significantly greater values in the upper crown than in the middle and lower crowns after June. WUE_sat only showed a slight vertical difference at PVI, PVII, and PIV.

Seasonal variation of photosynthetic and morphological parameters of needles

All photosynthetic and morphological parameters differed significantly among the individual sample trees and fluctuated during the growing seasons (Fig. 2). Mean daily A_max increased with time until late summer (at early August) to a maximum of nearly 9.42 μmol·m⁻²·s⁻¹ (Fig. 2a). For the remaining season, mean A_max ranged from 8.25 μmol·m⁻²·s⁻¹ to 8.34 μmol·m⁻²·s⁻¹. The mean daily R_D exhibited a significant decrease over time in early summer to a minimum of near 0.77 μmol·m⁻²·s⁻¹ and was then restored to 1.04 μmol·m⁻²·s⁻¹ (Fig. 2d). g_s-sat and NMA exhibited a similar time course with an increase during the growing season (Fig. 2b and 2c), but an abnormal peak appeared in early August. Mean WUE_sat demonstrated a time course that is opposite to that of g_s-sat and NMA (Fig. 2f). Mean NMA increased abruptly at the early period of needle expansion (PI, May), then remained stable until the second half of August (PVII), but finally increased at the end of growth.

Figure 2. 

Seasonal evolution of (a) light-saturated net photosynthetic rate (A_max); (b) light-saturated stomatal conductance (g_s-sat); (c) needle mass per area (NMA); (d) dark respiration (R_D); (e) ratio of length to width of needles (r_lw) and (f) light-saturated water use efficiency (WUE_sat) for five sample trees. Data points represent seasonal mean values (solid bars represented stand error). Black solid line represented corrective mean values of five sample trees (a)−(f).

Correlations between photosynthetic and meteorological parameters

The relationships between photosynthetic and main meteorological parameters are shown in Fig. 3. A_max significantly correlated to T_leaf, RH, and VPD in the entire treatment (Fig. 3a−c), in which A_max positively correlated to T_leaf and RH but negatively correlated to VPD. The correlation was stronger between A_max and RH (r = 0.46) than between A_max versus T_leaf (r = 0.27) and VPD (−0.28). R_D positively and linearly correlated with VPD (r = 0.31, Fig. 3i), but a stronger nonlinear relationship was observed between R_D and T_leaf (r = 0.61, Fig. 3g). RH poorly correlated with R_D (r = 0.07, Fig. 3h). g_s-sat significantly and nonlinearly correlated with RH (positive) and VPD (negative) (Fig. 3e & f). In contrast, WUE_sat negatively correlated with RH (r = −0.65, Fig. 3k), positively correlated with VPD (r = 0.64, Fig. 3l), and weakly correlated with T_leaf (r = 0.3, Fig. 3j).

Figure 3. 

Relationship between light-saturated net photosynthetic rate (A_max), light saturated stomatal conductance (g_s-sat), dark respiration (R_D), and light-saturated water use efficiency (WUE_sat) and (a), (d), (g), (j) leaf temperature (T_leaf); (b), (e), (h), (k) relative humidity (RH); and (c), (f), (i), (l) vapor pressure deficit (VPD) for five trees.

Relationships between physiological parameters and LMA

A_max and R_D exhibited a positive and significant correlation with LMA (Fig. 4a & c), but there were slight differences in the correlation coefficients among the individual sample trees. Although g_s-sat showed a similar correlation with LMA as A_max and R_D, the correlation was weaker (r = 0.35, Fig. 4b). WUE_sat only significantly correlated with LMA for two sample trees, even although a significantly negative relationship was observed for all the sample trees (Fig. 4d). WUE_sat negatively correlated with LMA but was more significant seasonally than spatially.

Figure 4. 

Relationships between (a) light-saturated net photosynthetic rate (A_max) and needle mass per area (NMA); (b) light-saturated stomatal conductance (g_s-sat) and NMA; (c) dark respiration (R_D) and NMA; (d) light-saturated water use efficiency (WUE_sat) and NMA. R values are the Pearson correlation coefficients. Solid lines represent the fitting result and are based on linear equations.

Previous studies have suggested that light^[15] and water potential^[19,41] are the most important factors affecting the vertical pattern of leaf physiology and morphology, but the primary driver between these two factors is still being debated^[42]. Recently, it has become increasingly accepted that the effects of these two factors vary with tree height^[43]. Light reportedly affects leaf functions and structures in short trees^[15,19], but for tall trees, a decrease in water potential considerably limits their leaf expansion and photosynthetic rate^[44]. Our results showed that A_max and g_s-sat decreased significantly from crown top to bottom (Fig. 1a & c), which is consistent with that of other studies^[45]. Martin et al. reported that shade leaves growing under less irradiance had lower leaf stomatal conductance than sun leaves^[18]. R_D negatively correlated with RDINC (Fig. 1b), as previously documented for different species^[46], because leaves generally adapted to dark environments by reducing R_D and Non-photochemical quenching. High leaf tissue density^[15] decreased mesophyll conductivity to gas, and restrained R_D^[45], which was also proved by the higher positive correlation between R_D and NMA (Fig. 4c). The vertical pattern of R_D partly decreased A_max with increasing tree height but was not significant. NMA is one of the main morphological traits that changes in response to light variations^[47]; thus, in our study, NMA followed the same pattern as A_max, g_s-sat, and R_D (Fig. 1e), suggesting that needles synthesize more photosynthetic tissue with increasing height to maximally use sufficient illumination. Studies revealed that the universal NMA gradient within the tree crowns or forest canopies is likely driven by solute content, leaf thickness^[48], leaf turgor pressure^[45], and leaf tissue density^[15], which reflect the plasticity and adaptability of foliage to the environment. Physiological variations in needles are usually accompanied by a corresponding change in their external form^[19]. Our results showed that r_lw significantly decreased with increasing RDINC (Fig. 1f), which further implied that trees maximized their photosynthetic efficiency by adjusting their foliage morphology to adapt to different environments in the vertical direction. Variations in WUE_sat originate from variations in photosynthetic rate, stomatal conductance, or both^[49]. Studies showed that WUE_sat negatively correlated with g_s-sat^[50]. In this study, WUE_sat showed an opposite vertical tendency to g_s-sat (Fig. 1d), which increased slightly and positively with RDINC, further supporting the opinion that foliage in the lower crown or canopy usually compensates for low resources by improving the utilization efficiency of site resources^[51].

Seasonal variation of needle physiology and morphology

Understanding the effect of seasonal variations on physiological and morphological parameters is critical for accurate modeling of carbon dioxide uptake by ecosystems, which can then be used to determine the magnitude of ecosystem carbon fluxes^[52]. Neglect of this variation may result in incorrect simulations of carbon uptake^[53]. Previous studies revealed that A_max and g_s-sat generally show a trend similar to a typical parabolic curve during the growing season^[22,54] although contrary results have been reported^[55]. However, A_max strongly correlated with g_s-sat in the above investigations, indicating that stomatal behavior has a pronounced impact on A_max. Our results show that all the physiological and morphological parameters fluctuated during the growth season. A_max had a parabolic seasonal variation, similar to that reported in other studies^[22] though accompanied by slight fluctuations in different trees (Fig. 2a), which was probably caused by the high correlation between A_max and seasonal variation of environment conditions (Fig. 3)^[38]. Kunert et al. confirmed that short-term exposure to high temperatures poses a considerable threat to conifer species in Central European forest production systems^[56]. During spring, an increase in A_max resulted from a gradual increase in photosynthetic capacity^[37]. A decrease in A_max during needle senescence is associated with a decrease in mesophyll conductivity to carbon dioxide owing to the increasing size of chloroplasts, starch grains, plastoglobuli, and the resorption of nitrogen^[57]. Seasonal variations in leaf photosynthetic traits, including maximum photosynthesis rate, maximum carboxylation rate, and mesophyll and stomatal conductance, can be well explained based on photoperiod variations^[2]. In addition, under both winter and drought stress, the main challenge for plants is that electron acceptor regeneration processes markedly slow down compared to primary photosynthetic processes, and this creates an imbalance between absorption and utilization of light energy^[28].

Previous studies have suggested that seasonal variations in R_D are mainly driven by seasonal patterns in temperature and NMA^[58]. It is well known that starch and soluble sugars are the main reactants in the respiratory process; thus, their content directly affects respiration. Temperature also indirectly limits respiration by affecting the activity of enzymes that participate in respiration^[59]. A previous study showed that a reduction in leaf expansion phase was due to a decrease in cytochrome-mediated respiration^[21]. Our results implied that R_D was significantly correlated with T_leaf (Fig. 3b) and NMA (Fig. 4b). The seasonal pattern of R_D showed an obvious reduction in the leaf expansion phase and then slightly recovered with a little fluctuation (Fig. 2b), following a similar trend in other studies^[60]. NMA showed a progressive increase throughout the growing season (Fig. 2e) owing to the accumulation of structural proteins and calcium. Seasonal variation in WUE_sat was different from that in g_s-sat (Fig. 2d), probably because WUE_sat and g_s-sat are negatively correlated^[49]. Previous studies on seasonal patterns of leaf length, width, and thickness have shown that they universally follow a saturation or parabolic curve^[58] throughout the growth season, reflecting the dynamic nature of photosynthetic acclimation^[61]. We also observed variations in leaf shape and found that the ratio of length to width (r_lw) showed a saturated tendency (Fig. 2f). The increase at the end of the growing season indicated that the needles started to senescence.

A further analysis of seasonal difference in photosynthetic rates among different canopy positions was conducted (Table 4), and the result indicated that A_max, WUE_sat and g_s-sat were significantly different in the upper crown but not significant in the lower crown. Conversely, R_D showed significant seasonal difference in the lower crown but not significant in the upper crown. Rare studys mentioned relevant results, but some research have proved that photosynthesis was more sensitive to light intensity and respiration was mainly affected by temperature^[49−51]. In our study, the closed-canopy caused an obvious seasonal change of light intensity in upper crown, but weak in lower crown. However, the seasonal variation of temperature was evident in the whole crown. Thus, A_max, WUE_sat and g_s-sat showed different seasonal difference compared to R_D.

Relationship between needle physiology and environment

T_leaf showed a significant parabolic correlation with A_max (Fig. 3a), corroborating the results of many other studies that focused on different tree species such as, Quercus crispul^[57], Picea mariana^[62], Pinus cembra^[63]. Both RH and VPD showed a significant relationship with A_max (Fig. 3e & i), particularly in the upper crown, presumably because needles in the upper crown are exposed to environmental stresses more frequently, and the variations in RH and VPD in the upper crown are more sensitive and greater^[64] than that in other crowns. R_D exhibited a typical exponential relationship with T_leaf (Fig. 3b), corroborating the results of other studies^[46,58,62]. Both g_s-sat and WUE_sat significantly correlated with RH (Fig. 3g & h) and VPD (Fig. 3k & l), but the tendencies were diametrically opposite. Similar processes have been observed in other studies^[37,50]. Some studies have suggested that variations in WUE_sat across dates are primarily driven by g_s-sat^[49] and our investigation showed that WUE_sat significantly correlated to g_s-sat (r = −0.71). Thus, we suggest that the relationship among WUE_sat, versus RH and VPD is likely caused by the influence of RH and VPD on g_s-sat.

Relationship between needle physiology and NMA

NMA plays an important role in predicting foliar physiological function, serves as a parameter in ecosystem modeling, and is used as an indicator for potential growth rate^[65]. Our results showed that A_max and g_s-sat both had a significant positive correlation with NMA (Fig. 4a & c), similar to the results of previous studies on other species^[66]. Other studies have shown a negative relationship between mass-based A_max and g_s-sat versus NMA, probably because of the vertical pattern of NMA^[44]. Han found a negative relationship between A_max and g_s-sat versus NMA for Pinus densiflora, probably because of a higher NMA value (> 200 g·m⁻²) than that in our study (< 130 g·m⁻²)^[67]. Moreover, a NMA value of up to 500 g·m⁻² for Pinus monticola ^[68], 800 g·m⁻² for Sequoia sempervirens ^[19], and 1,000 g·m⁻² for Pseudotsuga menziesii and Pinus ponderosa^[68]; almost all these species are categorized as tall tree species with a negative relationship between A_max and NMA. Therefore, we speculate that the relatively low NMA in our study may not be sufficient to limit mesophyll conductivity to carbon dioxide, and consequently, A_max and R_D is significantly and positively correlated with NMA, probably because of starch and soluble sugar contents^[58]. However, WUE_sat had a stronger correlation with RH and VPD than with NMA (Figs 3h, l & 4d), indicating that variations in WUE_sat across the growing season were primarily driven by the environment rather than by the needle morphology in this study.

Our study found that the spatial and seasonal variations of crown photosynthetic parameters for Larix olgensis were directly influenced by NMA, RH and VPD, in which NMA generally reflected the adaptability of leaves to the environmental factors. Thus, clarifying the response relationships between micro-environment and thinning intensity will contribute to the determination of optimal stand density. In addition, our results make it feasible to estimate the crown photosynthetic production and is helpful to further determine the contribution of branches to the trunk, which is the basis when making a pruning plan to produce no-knots wood and enhance the carbon sink capacity of young forests.

The authors confirm contribution to the paper as follows: study conception and design: Liu Q, Li F, Xie L; data collection: Liu Q, Xie L; analysis and interpretation of results: Liu Q, Zhang Z, Xie L; draft manuscript preparation: Liu Q, Wang D, Zhang Z, Xie L. All authors reviewed the results and approved the final version of the manuscript.

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.