Abstract

Background

   Study the leaf functional traits is highly important for understanding
   the survival strategies and climate adaptability of old trees. In this
   study, the old (over 100 years old) and mature trees (about 50 years
   old) of Pinus tabulaeformis in the Loess Plateau were studied, and the
   variation of 18 leaf functional traits (6 economic, 4 anatomical, 2
   photosynthetic and 6 physiological traits) was analyzed to understand
   the differences of survival strategies between old and mature trees.
   Combined with transcriptome and simple sequence repeats (SSR)
   techniques, the effects of soil property factors and genetic factors on
   leaf functional traits and the potential molecular mechanisms of traits
   differences were studied.

Results

   Compared with mature trees, old trees presented greater economic traits
   (except leaf phosphorus content), anatomical traits (except the
   stomatal density), and physiological traits (except superoxide
   dismutase activity) and lower photosynthetic traits, and their survival
   strategies were more conservative. The difference was mainly driven by
   soil property and genetic factors (common explanation rate was 67.89%),
   and the independent effect of genetic factors (10.09%) was slightly
   higher than that of soil property factors (2.88%). In addition, by
   constructing weighted gene co-expression networks analysis WGCNA), this
   research identified 24 candidate hub genes that regulate leaf
   functional traits, most of which are related to plant growth and
   development and the stress response, and can be used for further
   regulatory mechanism analysis.

Conclusions

   In conclusion, this study is helpful to understand the ecological
   adaptability of P. tabuliformis under the background of climate change
   in the Loess Plateau, and provides a theoretical basis related to leaf
   functional traits and molecular regulation for the protection of old
   trees.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12870-025-06130-8.

   Keywords: Arid and semi-arid regions, Environmental adaptation, Plant
   survival strategy, Trade-offs, WGCNA

Introduction

   The global climate has undergone significant transformation over the
   past century [[28]1], according to the Sixth Assessment Report of the
   United Nations Intergovernmental Panel on Climate Change (IPCC6), and
   human activities have caused a global warming of approximately 1.3 °C
   above preindustrial levels. Studies indicate that climate change is
   modifying the structure and function of forests globally [[29]2].
   However, because of the short growth time of common mature trees, the
   results do not adequately reflect the impact of long-term climate
   change on plants. Old trees (over 100 years old) are precious natural
   resources with high physiological, ecological and evolutionary values,
   and have attracted extensive attention from researchers in recent years
   [[30]3], especially in the study of plant climate adaptation [[31]4].
   Extremely old trees may contain information about climate changes over
   thousands of years [[32]5], and are extremely valuable materials for
   studying the climate sensitivity of forest growth [[33]6]. In the
   high-altitude forests of the Italian Alps, Larix decidua has
   age-related sensitivity to climate [[34]7]; when the climate warms, old
   trees have a stable growth rate and resistance to environmental changes
   compared to young and mature trees [[35]8]. In addition, due to the key
   role played by old trees in various ecosystem processes, such as
   biological carbon sequestration, the conservation of old trees is also
   receiving increasing attention [[36]9].

   Plant functional traits (PFTs) are relatively stable, measurable
   characteristics that influence plant growth, survival and reproduction,
   such as morphology, physiology and phenology, and can reflect plant
   responses to environmental changes [[37]10]. Leaves are the main organs
   of plants in direct contact with the environment, and their morphology,
   anatomy, physiology and other traits have strong plasticity and
   environmental sensitivity [[38]11, [39]12]. Generally, leaf functional
   traits are predictive of growth performance and ecological function
   [[40]13] and are important plant functional traits. For example, the
   leaf intracellular traits (physiological and biochemical parameters) of
   Populus euphratica under salt stress can be coordinated with the leaf
   economic and hydraulic traits to form a defense mechanism and reduce
   salt damage [[41]14]. Leaf functional traits analysis has become a
   basic tool for understanding the ecological model and function of
   plants, but in recent years, most research has focused on mature trees
   and has tended to analyze multiple traits of different organs [[42]15,
   [43]16]. Rawat et al. studied 13 leaf traits (3 morphological traits, 3
   chemical traits, 5 physiological traits and 2 stoichiometric traits) of
   9 tree species in Himalayan temperate forests, and accurately predicted
   the function of regional temperate forest ecosystem [[44]17]. And based
   on the measurement of 13 functional traits of leaves, petioles and
   branches of 101 broad-leaved trees in subtropical forests, it was found
   that specific leaf area and green leaf nutrient concentration were the
   main traits affecting leaf nitrogen and phosphorus uptake [[45]18].
   While these studies comprehensively illustrated the trade-offs for the
   whole plant, they lacked simultaneous assessments of the economic,
   anatomical, and biochemical properties of old trees and individual
   plant organs.

   For the first time, the leaf economic spectrum has quantified and
   summarized plant resource trade-off strategies through a series of
   closely related functional trait combinations (leaf nitrogen content,
   leaf phosphorus content, specific leaf mass, etc.), and quick
   investment-return-type plants generally have relatively high leaf
   nitrogen and phosphorus contents but relatively low leaf specific
   masses [[46]19]. Moreover, leaf dry matter content, leaf thickness, and
   leaf tissue density—traits linked to the capacity for resource
   acquisition, utilization, and conservation—have garnered significant
   attention from researchers [[47]20, [48]21]. The maximum leaf
   carboxylation rate is a pivotal parameter for assessing the
   photosynthetic capacity of plants, Chlorophyll content exhibits a
   significant positive correlation with the maximum leaf carboxylation
   rate [[49]22], and the relative leaf water content can better predict
   the photosynthetic rate of leaves [[50]23]. Leaf anatomical traits
   underpin leaf morphological variations and a range of physiological and
   biochemical processes [[51]24, [52]25]. Furthermore, the investigation
   of plant physiological traits enhances our comprehension of
   intracellular regulatory mechanisms, including stress-related enzymes
   (e.g., superoxide dismutase) and osmoregulatory substances (e.g.,
   proline) [[53]26, [54]27].

   The primary determinants of leaf functional traits are environmental
   and genetic factors. Leaf functional traits are highly responsive to
   environmental factors, including water, radiation, and temperature, and
   can be modulated to optimize the utilization of limited resources and
   enhance the survival capacity of plants [[55]28, [56]29]. The evolution
   of plant traits is influenced by phylogeny, and closely related
   individuals may be more similar in terms of leaf functional traits and
   ecological adaptability [[57]30]. However, evidence also suggests that
   variation in functional traits is weakly influenced by phylogeny
   [[58]31], but neglecting phylogenetic effects can result in
   overestimations of the number of independent observations and
   inaccurate assessments of relationships and trade-offs between traits
   [[59]32]. Therefore, combining environmental and genetic factors can
   help to fully reveal the regulatory mechanism of trait variation.

   Simple sequence repeat (SSR) markers are ideal tools for genetic
   diversity analysis, for example, with the help of SSR markers,
   functional SSR markers and provenances in natural Dalbergia odorifera
   populations were revealed [[60]33]. The development of SSR loci has
   also attracted extensive attention from researchers. Yang et al. used a
   large number of SSR markers selected from the transcriptome data of P.
   tabuliformis to select 24 optimal SSR markers and evaluate the genetic
   diversity of each generation in the seed orchard [[61]34]. Weighted
   gene co-expression network analysis (WGCNA) is a systems biology method
   used to describe patterns of gene association between different
   samples. It can be combined with transcriptome sequencing technology to
   identify highly covaried gene sets, and identify candidate biomarker
   genes based on the interconnectivity of gene sets and the association
   between gene sets and phenotypes, and is currently widely used to
   identify key genes related to plant phenotypic traits, responses to
   biological or abiotic stresses, and other aspects [[62]35]. For
   example, Ji et al. found 2 modules and 6 hub genes significantly
   associated with maize leaf phenotypes through this technique [[63]36].

   The Loess Plateau is located in the semihumid to semiarid transition
   zone in the temperate zone. It is famous for its fragile ecology and
   severe soil erosion [[64]37]. Pinus tabuliformis Carrière is a
   significant afforestation species on the Loess Plateau, known for its
   commendable ecological functions, including soil and water conservation
   [[65]38]. At present, in the Huanglong area, through afforestation
   projects with different periods, there are large areas of P.
   tabuliformis artificial and natural secondary forests (both belong to
   mature trees), which play important ecological roles in regional soil
   and water conservation [[66]39]. However, owing the large fluctuations
   in annual precipitation in the Huanglong, the stability of forest
   ecosystems is threatened; therefore, clarifying the adaptive mechanisms
   of P. tabuliformis is urgently needed. In addition, the study of the
   differences in functional traits between old and mature trees will help
   us better understand the survival strategies of plants in the context
   of climate change [[67]40, [68]41]. In this study, the leaf economic,
   photosynthetic, anatomical and physiological traits of old and mature
   trees (artificial and natural secondary forests) of P. tabuliformis in
   Huanglong area were studied. The aims were to: (1) investigate whether
   there are differences in survival strategies on the basis of leaf
   functional traits between old and mature trees; (2) explore the driving
   factors of the differences in leaf functional traits; and (3) screen
   key regulatory genes via transcriptome sequencing technology. The
   findings of this study contribute to the elucidation of the adaptive
   mechanisms of P. tabuliformis from the perspective of leaf functional
   traits, offering a scientific basis for the rational management of
   forests and the preservation of old trees.

Materials and methods

Study area and sample

   The study area is located in Muchangqiao (E109°41′37″, N35°38′26″) and
   Caijiachuan (E109°56’19″, N35°54’1″) in Huanglong County (Yan’an city),
   Shaanxi Province (Fig. [69]S1). It is located within a temperate zone
   and is distinguished by a continental monsoon climate. The soil
   composition is predominantly mountain brown, the altitude is
   approximately 1842.4 m, the average annual precipitation (MAP) stands
   at 561.4 mm, with 66.7% of the total precipitation occurring in the
   summer months, and the average annual temperature (MAT) is 12.1 ℃
   [[70]42]. Mature trees of about 50 years old were selected. Since the
   Chinese government launched a nationwide restoration project in 1998,
   natural secondary forests of Pinus tabulaeformis Carr. have been
   effectively protected, and the artificial forests have a history of
   approximately 50 years [[71]43, [72]44]. The vegetation under the
   forests is mainly shrubs (such as Lonicera japonica Thunb., Celastrus
   orbiculatus Thunb., Prunus tomentosa Thunb., legume, etc.) and herbs
   (such as Carex spp., Artemisia argyi H. Lév. & Vaniot, etc.). The old
   tree group is located mainly in Muchangqiao village and is
   approximately 150–350 years old. It is close to the village and road,
   and is severely disturbed by human activities, with only a few herbs
   under the forest (Carex spp.). In this study, the age of all P.
   tabuliformis samples was determined using growth cones.

   At the end of June 2023 (after 20 days without rain), 6 fixed 20*20 m
   plots were established in the artificial forests (HR) and natural
   secondary forests (HT), 3 standard trees were selected for sampling in
   each plot, and each plot was treated with mixed samples. Considering
   that there were few old trees, 12 old trees (HG) in good health and
   away from the edge of the stand were randomly selected as research
   objects under the sampling principle. A total of 24 samples were
   collected, including 12 old trees, 6 artificial forests and 6 natural
   secondary forests, see Table [73]S1 for details. In order to better
   represent the average condition of the trees as a whole, annual healthy
   leaves were collected from four directions (east, west, north, and
   south) in the middle of the crown of each tree. With trees as the
   origin, soil samples of 0–20 cm were collected within a radius of 1 m
   in four directions [[74]45]. Soil and leaf samples collected from all
   four directions of the same standard tree need to be mixed. A total of
   24 soil and leaf samples were collected. The leaf samples were divided
   into 2 parts: one part stored at 4 ℃ for morphological and chemical
   determination, the other was placed on dry ice and stored at -80 ℃ for
   leaf physiological index determination and transcriptome sequencing.
   The soil samples were air-dried for soil physicochemical property
   analysis.

Analysis of soil physical and chemical properties

   The physical and chemical properties of the soil samples were
   determined via conventional methods [[75]45]. The soil pH was
   determined potentiometrically in a 1:2.5 soil-to-water solution via a
   glass electrode (PHSJ-4 A pH Acidometer, Shanghai, China). The soil
   moisture content (W, %) was calculated by weight loss when fresh soil
   was baked in a 105 ℃ oven to a constant weight (12–14 h), it was
   calculated as:
   graphic file with name M1.gif

   where, m[0] is the weight of soil before drying and m[1] is the weight
   of soil after drying. The content of the soil organic matter (SOC,
   g·kg^− 1) was determined via an external heating method that uses of
   potassium dichromate, according to the formula:
   graphic file with name M2.gif

   In the formula, c is the concentration of
   0.8000 mol·L^−1(1/6 K[2]Cr[2]O[7]) standard solution; 5 is the volume
   of potassium dichromate standard solution added (mL); V[0] is the
   volume (mL) of blank titration with FeSO[4] removed; V is the volume of
   FeSO[4] removed for sample titration (mL); 3.0 is the molar mass of 1/4
   carbon atom (g·mol·L^− 1); 1.1 is the oxidation correction factor; m is
   the quality of air-dried soil sample (g); k is the coefficient of
   converting air-dried soil into dried soil; 1.724 is the average
   conversion factor of soil organic carbon to soil organic matter. The
   soil total nitrogen content (TN, g·kg^−1) and the soil total phosphorus
   content (TP, g·kg^− 1) were determined by an AA3 continuous flow
   analyzer (AA3, SEAL Analytical, Germany) after digestion with
   H[2]SO[4]-H[2]O[2]. The TN and TP can be calculated as,
   graphic file with name M3.gif

   where, ρ is the mass concentration of the liquid to be measured
   (ug·mL^−1); V is the measured volume (mL) and m is the mass of dried
   soil sample (g). The soil ammonium nitrogen (NH[4]^+-N, mg·kg^−1) and
   nitrate nitrogen (NO[3]^−-N, mg·kg^− 1) were extracted with 1 M KCl and
   determined with an AA3 continuous flow analyzer (AA3, SEAL Analytical,
   Germany). The formula for calculating NH[4]^+-N is:
   graphic file with name M4.gif

   where, ρ is the mass concentration of the liquid to be measured
   (mg·L^− 1); V is the measured volume (mL) and m is the mass of dried
   soil sample (g), the calculation formula for NO[3]^−-N is the same as
   that for NH[4]^+-N. The soil available phosphorus (AP, mg·kg^− 1) was
   measured via the molybdenum blue colorimetric method after the samples
   were extracted with 0.5 M Na[2]CO[3]. It is calculated as:
   graphic file with name M5.gif

   where ρ is the mass concentration of P (mg·L^− 1) obtained from the
   standard curve; V is the constant volume (mL) at the time of color
   development; ts is the fraction multiple (that is, the ratio of the
   total volume of the extract to the volume of the extracted liquid at
   the time of color development); m is the mass of air-dried soil (g) and
   k is the coefficient that converts air-dried soil into dried soil mass.
   The soil total potassium content (TK, g·kg^− 1) was determined via a
   flame photometer (Flame photometer, Sherwood, Britain) after digestion
   with H[2]SO[4]-H[2]O[2]. The soil available potassium (AK, mg·kg^− 1)
   was shaken for 30 min with 1 M ammonium acetate solution (1:10 w/v) and
   then analyzed via flame photometry, AK and TK contents are calculated
   according to the standard curve.

Determination of leaf functional traits

   All indicators and abbreviations used in the study are shown in
   Table [76]1.

Table 1.

   Summary of leaf functional traits
   Type Trait Abbreviation Unit Function
   Leaf economic traits contents of nitrogen LNC g·kg^− 1 Comprehensive
   parameters to characterize plant growth and development
   contents of phosphorus LPC g·kg^− 1
   leaf thickness LT mm It is related to physiological activities such as
   nutrient absorption
   leaf mass per area LMA g·m^− 2 Reflects the trade-off between carbon
   gain and leaf life
   leaf dry matter content LDMC g·g^− 1 An important indicator of plant
   utilization of habitat resources
   leaf tissue density LTD g·cm^− 3 Reflects plant growth input strategy
   Leaf photosynthetic traits leaf relative water content LRWC % Reflects
   the photosynthetic status of trees
   chlorophyll content Chl mg·g^− 1
   Leaf anatomical traits stomatal number SN pcs Reflects the respiration
   and transpiration capacity of needles
   stomatal density SD pcs·mm^− 2
   resin canal number RCN pcs It reflects the transport capacity of
   secondary metabolites
   vascular tissue area VBA mm^2 It reflects the growth of the channeling
   tissue
   Leaf physiological traits superoxide dismutase activities SOD
   U·g^− 1FW·h^− 1 Two important protective enzymes used by plants to
   remove reactive oxygen species
   peroxidase activities POD u·g^− 1min^− 1
   contents of malondialdehyde MDA umol·g^− 1 Measuring the damage to the
   cell membrane of trees
   contents of proline Pro µg·g^− 1 Important osmoregulatory substances
   soluble protein SP mg·g^− 1
   soluble sugar SS %
   [77]Open in a new tab

   Leaf economic traits: Similar to those in the soil determination
   method, LNC and LPC were quantified using the method formulated by Bao
   [[78]36]. LT is the average leaf thickness at 1/4, 1/2, and 3/4. LMA
   and LDMC were determined based on the procedures delineated by
   Pérez-Harguindeguy et al. [[79]46]. LTD = leaf dry weight / leaf volume
   [[80]47].

   Leaf photosynthetic traits: LRWC was measured according to the method
   by Gadallah [[81]48]. Chl was determined via the 80% acetone extraction
   method [[82]49].

   Leaf anatomical traits: The methodology of Zeng et al. was employed to
   ascertain stomatal traits [[83]50]. Each needle was treated with
   transparent nail polish at the same position (2 cm away from the base).
   After the nail polish dried naturally. The oil film was carefully
   removed using tweezers and affixed to a slide, the SN on the 2 mm
   stomatal line (including concave and convex sides) was observed through
   a fluorescence microscope (Fluorescent biological microscope, COIC,
   China), and the statistics were performed via I[MAGE]J processing
   software (v.1.54i). The SD was the number of stomata per unit area
   divided by the corresponding field of view area. Needles 0.5 cm in
   length were cut 2 cm from the base and put into FAA fixing solution,
   and were manually sliced to make permanent sections. Images were
   obtained via fluorescence microscopy (Fluorescent biological
   microscope, COIC, China), and RCN and VBA were measured and analyzed
   via I[MAGE]J processing software (v.1.54i).

   Leaf physiological traits: SOD activity was determined via the nitro
   blue tetrazolium method [[84]51]. POD activity was determined via the
   guaiacol oxidation method [[85]52]. MDA was determined via the
   thiobarbituric acid reaction [[86]53]. Pro was analyzed via the
   ninhydrin colorimetric method [[87]54]. The SP content was determined
   via the Coomassie Brilliant Blue G-250 method [[88]55]. The SS content
   was analyzed via the anthrone colorimetric method [[89]56].

Genotyping and genetic diversity analysis

   DNA was extracted from each leaf sample via a DNA extraction kit
   (BioTeke, Beijing, China). The SSR primers directly referenced 8 pairs
   of primers with stable amplification and high polymorphism rates
   obtained in previous studies [[90]57] (Table [91]S2). An ABI-Veriti
   96-well gradient PCR was used for DNA amplification. The PCR system was
   25 µl, and the reaction procedure for PCR amplification as follows:
   predenaturation at 98 ℃ for 2 min; denatured for 10 s at 98 ℃;
   annealing at 50 ~ 60 ℃ for 15 s, delay at 72 ℃ for 20 s; 33 cycle; and
   finally, extension for 5 min at 72 ℃ and storage at 4 ℃. The PCR
   products were then subjected to capillary electrophoresis. The original
   fluorescence detection data were analyzed with GeneMarker v.1.65
   (Applied Biosystems) software, the original data of each sample were
   read, and the corresponding results were recorded in Excel. The number
   of alleles (N[a]), expected heterozygosity (H[e]), Shannon index (I),
   interpopulation Nei’s genetic distance (D), and inbreeding coefficient
   (F[is]) were calculated using GenAlex v.6.502 software. The
   polymorphism information index (PIC) was calculated via PowerMarker
   v.3.25 software, and cluster analysis was performed on the basis of
   Nei’s genetic distance. The evolutionary tree was constructed via MEGA
   (v.6.0) software.

   The 8 primers utilized in the study amplified an average of 2.5 alleles
   per primer pair, with the exception of J20 and J29, for which the PIC
   values of the other primers exceeded 0.25, indicating moderate
   polymorphism, which could effectively reveal the genetic diversity of
   P. tabuliformis germplasm resources. More parameters related to the
   genetic diversity of the primers are shown in Table [92]S3.

RNA sequencing, weighted correlation network analysis and gene network
visualization

   Total RNA was extracted from frozen conifer samples, and analyzed by
   Personalbio Biotechnology Co., Ltd. (Shanghai, China). Results from Niu
   et al. were used as a reference genome [[93]58]. Following the
   exclusion of genes with low expression levels (FPKM < 1), weighted gene
   co-expression network analysis (WGCNA) was conducted using the “WGCNA”
   R package to analyze expression data from 24 different samples. The
   soft threshold was set to 6, mergeCutHeight to 0.2, minModuleSize to
   100, and deepSplit to 2. Cytoscape (v.3.10.1) was then used to
   visualize the gene networks within the module and present the
   biological interactions of the hub genes.

Validation of hub genes

   Both RNA and cDNA were extracted using TIANGEN (Beijing, China) kits.
   PtACTIN12 of P. tabuliformis was used as an internal reference gene
   [[94]59]. All primer sequences are listed in Table [95]S4. Real-time
   fluorescence quantitative PCR assay system was used to determine the
   results.

Statistical analysis

   A one-way analysis of variance (ANOVA) and a least significant
   difference (LSD) multiple comparison test were conducted using IBM SPSS
   Statistics (v.26.0.0.0) to ascertain the differences in leaf functional
   traits between old and mature trees. The Pearson correlation
   coefficient was employed to calculate the correlation between leaf
   functional traits and the correlation between leaf functional traits
   and soil property factors. The correlation heatmap was generated using
   the “Corrplot” and “ggplot2” R packages. Principal component analysis
   (PCA) was conducted using the “FactoMineR” R package to elucidate the
   synergistic and trade-off relationships among leaf functional traits.
   Redundancy analysis (RDA) was employed to assess the correlation
   between leaf functional traits and soil property variables across the
   three populations, with the significance of each soil property factor
   determined through a permutation test. Variance decomposition analysis
   (VPA) was used to quantify the proportions of soil property and genetic
   factors that explained the variation in leaf functional traits. Except
   for ANOVA and LSD, the remaining calculations were obtained using R v.
   4.0.2.

Results

Changes in leaf functional traits of old and mature trees

   Differences in leaf functional traits were observed between old and
   mature trees, with the interspecific CV values for 18 leaf functional
   traits ranging from 2.02 to 55.82% (Fig. [96]1; Table [97]S5). The
   economic traits of old trees were significantly higher than those of
   mature trees except LPC (P < 0.05) (Fig. [98]1a).The photosynthetic
   traits of old trees were notably lower than those of mature trees,
   while certain anatomical traits (RCN and VBA) were significantly
   greater than those of mature trees (P < 0.05) (Fig. [99]1b, c).
   Regarding physiological traits, with the exception of SOD, the
   physiological traits of old trees were higher to those of mature trees,
   and significant differences were observed in the majority of these
   traits (P < 0.05) (Fig. [100]1d).

Fig. 1.

   [101]Fig. 1
   [102]Open in a new tab

   Histogram of (a) economic traits, (b) photosynthetic traits, (c)
   anatomical traits, and (d) physiological traits of 3 populations.
   Different letters indicated significant differences in leaf traits
   among different populations (LSD test, P < 0.05). LNC, leaf nitrogen
   content; LPC, leaf phosphorus content; LT, leaf thickness; LMA, leaf
   mass per area; LTD, leaf tissue density; LDMC, leaf dry matter content;
   LRWC, leaf relative water content; Chl, chlorophyll content; SD,
   stomatal density; SN, stomatal number; RCN, resin canal number; VBA,
   vascular tissue area; SOD, superoxide dismutase activity; POD,
   peroxidase activity; MDA, malondialdehyde content; Pro, proline
   content; SP, soluble protein content; SS, soluble sugar content

Correlations between leaf functional traits

   Differences were observed in the correlation and trade-off models of
   functional traits between old and mature trees (Figs. [103]2 and
   [104]3). In old trees, LT, LMA and LTD were strongly correlated with
   leaf anatomical traits (except SD) (P < 0.05), whereas LDMC was
   negatively correlated with photosynthetic traits (LWRC) (P < 0.001)
   (Fig. [105]2a). PCA of the leaf functional traits on the first two axes
   revealed that the interpretative variances of the two axes were 28.59%
   and 17.43%, respectively. LT, LTD, LMA, SN, VBA and RCN contributed
   more to the first principal component axis, whereas LTD had the
   opposite synergistic direction with the other traits. LDMC and LRWC
   contributed more to the second principal component axis, and the
   direction of coordination was the opposite (Fig. [106]3a). Artificial
   and natural secondary forests showed different synergies and trade-off
   modes. Compared with old trees, leaf functional traits of artificial
   forests were more correlated (Fig. [107]2b), some economic traits (LPC,
   LMA and LDMC) and photosynthetic traits contributed more to the first
   principal component axis (Fig. [108]3b). The correlation between traits
   in natural secondary forests was weaker than that in artificial
   forests, but stronger than old trees (Fig. [109]2c).

Fig. 2.

   [110]Fig. 2
   [111]Open in a new tab

   Heat maps of leaf functional traits of (a) old trees, (b) artificial
   forests, and (c) natural secondary forests. ***P < 0.001; **P < 0.01;
   *P < 0.05. LNC, leaf nitrogen content; LPC, leaf phosphorus content;
   LT, leaf thickness; LMA, leaf mass per area; LTD, leaf tissue density;
   LDMC, leaf dry matter content; LRWC, leaf relative water content; Chl,
   chlorophyll content; SD, stomatal density; SN, stomatal number; RCN,
   resin canal number; VBA, vascular tissue area; SOD, superoxide
   dismutase activity; POD, peroxidase activity; MDA, malondialdehyde
   content; Pro, proline content; SP, soluble protein content; SS, soluble
   sugar content

Fig. 3.

   [112]Fig. 3
   [113]Open in a new tab

   Principal components analyses (PCA) of community-level (a) old trees,
   (b) artificial forests, and (c) natural secondary forests. LNC, leaf
   nitrogen content; LPC, leaf phosphorus content; LT, leaf thickness;
   LMA, leaf mass per area; LTD, leaf tissue density; LDMC, leaf dry
   matter content; LRWC, leaf relative water content; Chl, chlorophyll
   content; SD, stomatal density; SN, stomatal number; RCN, resin canal
   number; VBA, vascular tissue area; SOD, superoxide dismutase activity;
   POD, peroxidase activity; MDA, malondialdehyde content; Pro, proline
   content; SP, soluble protein content; SS, soluble sugar content

Driving factors of leaf functional traits variation

   Traits are usually determined by a combination of soil property and
   genetic factors. Through RDA, the explanatory values of the first two
   axes were 69.42% and 6.66%, respectively. A significant relationship
   was observed between leaf functional traits and W, SOC, TP and C/N
   (P < 0.05), which are the primary soil property factors contributing to
   the variation in leaf functional traits (Fig. [114]4a).

Fig. 4.

   [115]Fig. 4
   [116]Open in a new tab

   Analysis of driving factors of difference of leaf functional traits.
   (a) RDA ranking of leaf functional traits by environmental factors.
   ***P < 0.001; **P < 0.01; *P < 0.05. (b) Principal coordinate analysis
   for 3 populations of P. tabuliformis. (c) UPGMA dendrogram of P.
   tabuliformis germplasm resources. (d) VPA analysis of soil property and
   genetic factors. W, the soil moisture content; SOC, the content of the
   soil organic matter; TP, the soil total phosphorus content; TK, the
   soil total potassium content; AP, the soil available phosphorus; AK,
   the soil available potassium; NO[3]^−-N, the soil nitrate nitrogen;
   NH[4]^+-N, the soil ammonium nitrogen; pH, soil pH; C/N, ratio of soil
   carbon to nitrogen content

   SSRs were employed to further explore the genetic variation of P.
   tabuliformis across different populations. Genetic diversity varied
   among different populations (Table [117]2). Molecular analysis of
   variance (AMOVA) revealed that 81.41% of the variation was within the
   population (Table [118]3). Cluster analysis was performed according to
   Nei’s genetic distance (Fig. [119]4c), when the genetic distance was
   0.028, the different populations were divided. Old and mature trees are
   independent of each other and have a long genetic distance
   (Fig. [120]4b).

Table 2.

   Parameters of genetic diversity for 3 populations of P. tabuliformis
   Populations N [a] N [e] I    H [o] H [e] F [is]
   HG           2.00  1.63 0.53  0.39  0.35  -0.12
   HR           1.88  1.36 0.38  0.31  0.23  -0.33
   HT           2.00  1.42 0.43  0.38  0.27  -0.39
   Mean         1.96  1.47 0.45  0.36  0.28  -0.28
   [121]Open in a new tab

   N[a], observed number of alleles; N[e], effective number of alleles; I,
   Shannon’s information index; H[o], observed heterozygosity; H[e],
   expected heterozygosity; F[is], inbreeding coefficient at the
   population level

Table 3.

   Genetic variation among and within populations of P. tabuliformis based
   on AMOVA analysis
   Source of variation Degree of freedom Sum of squares Mean squares
   Percent of total variation
   Among populations           2             11.54          5.77     0.19
   Within populations         21             44.67          2.13     0.81
   Total                      23             56.21                   1.00
   [122]Open in a new tab

   The contribution of soil property factors and genetic factors to leaf
   functional traits was evaluated by VPA method (Fig. [123]4d). Genetic
   factors accounted for 10.09% of the variation in leaf functional
   traits, approximately three times greater than that of soil property
   factors (2.88%), and the common explanation of both factors was 67.89%.

Construction and analysis of the gene co-expression network

   In this study, PCA and cluster analysis were performed on the
   transcriptome data, and the results revealed that the old and mature
   trees were divided into two groups (Fig. [124]S2). To further
   understand the regulation of leaf functional trait changes in old and
   mature trees, WGCNA analysis was performed. Transcriptome sequencing
   was performed on 24 samples, filtered out genes whose average FPKM
   expression was lower than 1, obtained 23,308 genes. A similarity
   threshold of Fold > 0.2 for module fusion to construct a co-expression
   network module. A total of 28 expression modules were identified based
   on similar expression patterns (Fig. [125]5a). The heatmap of
   module-trait correlations (Fig. [126]5b) indicated that most of the
   economic traits, anatomical traits and physiological traits presented
   significant positive correlations with the MEblue module and MEsalmon
   module genes, whereas the photosynthetic traits presented the greatest
   correlations with the MEdarkred module genes. In this study, KEGG
   enrichment analysis was performed for genes significantly associated
   with these three modules (Fig. [127]6) and discovered that the genes
   within the MEblue module are predominantly involved in vesicular
   transport; N-glycan biosynthesis; D-amino acid metabolism; proteasome,
   pantothenate and CoA biosynthesis; and 10 other kinds of life
   activities (Fig. [128]6a). The genes within the MEdarkred module are
   involved mainly in Sphingolipid metabolism; fatty acid elongation;
   steroid biosynthesis; lipoic acid metabolism; glycosaminoglycan
   degradation and 13 other kinds of life activities (Fig. [129]6b). The
   genes within the MEsalmon module are implicated in aminoacyl-tRNA
   biosynthesis; carotenoid biosynthesis; nicotinate and nicotinamide
   metabolism (Fig. [130]6c). These findings suggest that the genes within
   these three modules are primarily associated with the variations in the
   functional traits of P. tabuliformis needles and are implicated in
   crucial life processes such as plant growth and development, as well as
   stress resistance.

Fig. 5.

   [131]Fig. 5
   [132]Open in a new tab

   Correlation analysis of leaf functional traits and transcriptomics. (a)
   Dendrogram showing co-expression modules (clusters) identified by
   WGCNA. The major tree branches constitute 28 modules labeled with
   different colors. (b) Heat maps showing module-trait correlations. Each
   row corresponds to a module in a different color. Each column
   corresponds to a functional trait. Red color indicates a positive
   correlation between the cluster and the tissue. Blue color indicates a
   negative correlation

Fig. 6.

   [133]Fig. 6
   [134]Open in a new tab

   Pathway enrichment analysis of three important modules. (a) MEblue
   module, (b) MEdarkred module, and (c) MEsalmon module

Analysis of the interaction network of hub genes in the module

   In this study, gene network visualization and gene connectivity
   analysis for the MEblue, MEsalmon, and MEdarkred module genes. The
   network visualization map of genes with the top 200 weight values in
   the module was drawn via Cytoscape software, and genes with a
   degree > 50 were selected as the hub genes (Fig. [135]7; Table [136]4).
   Among the 10 hub genes in the MEblue module, Pt3G28290 is a gene with
   an unknown function, while other genes are involved in oxidoreductase
   (Pt2G33120), material transport (Pt2G50170, Pt5G28900, PtXG14450, and
   Pt6G34200), secondary metabolite synthesis (Pt7G63210, Pt7G22770, and
   Pt8G05170), and endoplasmic reticulum formation (Pt6G53630)
   (Fig. [137]7a). Among the 10 hub genes in the MEdarkred module,
   PtJG20980 is a gene with unknown function, whereas the other genes are
   involved in energy synthesis (Pt1G44120, Pt2G56940, and Pt7G27100),
   ubiquitination (Pt7G44390, Pt5G35270, and Pt5G44890), lipid metabolism
   (Pt7G00280), and stress response (Pt3G40810 and Pt3G08250)
   (Fig. [138]7b). In the MEsalmon module, Pt2G64380 and Pt4G57300 are
   genes with an unknown function, and Pt3G16430 is an immune protein.
   Pt7G16010 is related to the metabolism of niacin and niacinamide, and
   participates in various metabolic activities in plants (Fig. [139]7c).

Fig. 7.

   [140]Fig. 7
   [141]Open in a new tab

   Analyzing the interactions between hub gene networks within the
   co-expression module. (a) Interaction analysis of core genes in the
   MEblue module, (b) Interaction analysis of core genes in the MEdarkred
   module, and (c) Interaction analysis of core genes in the MEsalmon
   module

Table 4.

   Hub genes related to leaf functional traits
   Module Gene ID Description
   MEdarkred Pt2G33120 Glucose-methanol-choline (GMC) oxidoreductase
   family protein
   Pt3G28290 -
   Pt2G50170 N-MYC downregulated-like 3
   Pt5G28900 lipid-transfer protein
   Pt8G05170 galactan synthase 2
   PtXG14450 ATP-binding cassette B2
   Pt7G22770 3-ketoacyl-CoA synthase
   Pt6G34200 solute carrier family 35, member F5
   Pt6G53630 Reticulon family protein
   Pt7G63210 S-adenosyl-L-methionine-dependent methyltransferases
   superfamily protein
   Meblue Pt1G44120 proton pump interactor 1
   Pt7G44390 F-box and leucine-rich repeat protein 2/20
   Pt2G56940 Phosphoribulokinase / Uridine kinase family
   Pt7G00280 fatty acid omega-hydroxylase
   PtJG20980 -
   Pt7G27100 serine/threonine-protein kinase 24/25/MST4
   Pt3G40810 zinc finger (C2H2 type) family protein
   Pt5G35270 Ubiquitin carboxyl-terminal hydrolase family protein
   Pt3G08250 Endosomal targeting BRO1-like domain-containing protein
   Pt5G44890 E3 ubiquitin-protein ligase Arkadia
   Mesalmon Pt2G64380 -
   Pt3G16430 Mannose-binding lectin superfamily protein
   Pt4G57300 -
   Pt7G16010 pyrimidine and pyridine-specific 5’-nucleotidase
   [142]Open in a new tab

Validation of the hub genes via RT‒qPCR

   To validate the reproducibility and authenticity of the RNA-seq data, 8
   genes with large fold change (|Fold Change|>1.5) were randomly selected
   from the 24 hub genes obtained by WGCNA analysis for qRT-PCR analysis
   (Table [143]S7). The qRT-PCR results for these 8 genes aligned with the
   expression patterns observed in the RNA-seq data (Fig. [144]S4).

Discussion

The survival strategies of old and mature P. tabuliformis differ

   To adapt to different habitat conditions, plant functional traits and
   their combinations differ, resulting in different survival strategies
   [[145]60]. Our findings indicated that variations existed in leaf
   functional traits between old and mature trees. (Fig. [146]1).
   According to the leaf economic spectrum, old trees resource utilization
   strategies are more conservative [[147]19]. Generally, when the
   survival strategy tends to be conservative, the leaf photosynthetic
   trait value decreases due to the decrease of plant transpiration and
   photosynthesis. However, the value of traits associated with the
   construction of vascular bundle tissue will increase, because the
   developed vascular bundle tissue can ensure efficient transport and
   adequate water and mineral supplies, thereby enhancing the resilience
   to environmental stress [[148]24]. Our results are consistent with this
   conclusion (Fig. [149]1b, c). Moreover, previous studies have suggested
   that leaf photosynthetic and anatomical traits are usually coupled with
   leaf economic traits [[150]61–[151]63]. However, in our study, only the
   leaf photosynthetic traits of artificial forests exhibited a strong
   correlated with leaf economic traits (Fig. [152]2). This may be due to
   the fact that the living conditions of old trees and natural secondary
   forests are worse than those of artificial forests, and the overall
   living environment of artificial forests is better under thinning and
   other tending measures, while natural secondary forests have high stand
   density and strong inter-specific competition, and old trees are
   seriously affected by human interference and age factors, so
   multi-dimensional trait relationships are needed to maintain their
   survival. Multiple trait dimensions may enhance species’ adaptation to
   multiple ecological niche dimensions, thereby fostering species
   coexistence [[153]64]. Leaf physiological traits are frequently
   employed as potential indicators of plant survival strategies [[154]65,
   [155]66]. When plants are stressed, a large number of reactive oxygen
   species are produced in the cells, which damages the integrity of the
   cell membrane [[156]67]. In this case, cells possess a timely defense
   system in response to stimuli, such as increased peroxidase activity
   and soluble substance content, which can effectively reduce ROS levels
   and maintain a stable growth state in plants (Fig. [157]1d). Notably,
   SOD constitutes the primary line of defense in the enzymatic clearance
   system, with its activity progressively diminishing under severe plant
   stress [[158]68]. When old trees are subjected to the same
   environmental pressure as mature trees, they are also subjected to many
   self-generated stresses, such as aging, which may be one of the reasons
   for the decrease in the SOD value of old trees. In conclusion, our
   results fully illustrate the first goal of our study, namely, that old
   trees have more conservative survival strategies, and our study further
   confirms that leaf functional traits have integrity, that different
   traits can be related to each other, and that plants should be regarded
   as interconnected organisms. This method of exploring trait variation
   on the basis of comprehensive traits is highly important.

Leaf functional traits are driven by both soil property and genetic factors

   Our study systematically examined the impacts of soil property and
   genetic factors on the leaf functional traits of P. tabuliformis, and
   the findings indicated that these traits were modulated by both soil
   property and genetic influences (Fig. [159]4). In terms of soil
   property factors, the RDA results demonstrated that W and TP were the
   main driving factors of leaf functional traits variation
   (Fig. [160]4a), with an increase in the TP content promoting increases
   in leaf economic, anatomical, and physiological traits, whereas the
   effect of W was the opposite (Fig. [161]S3). Changes in the soil water
   supply can alter rooting patterns, which directly affect plant anatomy
   and function, and thus affect the construction of other functional
   traits of plants [[162]69]. TP is closely associated with plant growth
   and ecosystem productivity, and plant species exhibiting relatively
   elevated leaf nutrient concentrations and metabolic rates under
   conditions of unrestricted TP availability [[163]70].The action mode of
   W was similar to that of previous studies [[164]71], but there were
   differences in the TP. Cui et al. reported that species with relatively
   high SLA, LNC and LPC values are relatively common in soils with high
   phosphorus contents, whereas species with relatively high LDMC values
   are relatively common in soils with low phosphorus contents [[165]72].
   This difference may be related to the severe phosphorus deficiency
   caused by highly weathered soil in the Loess Plateau [[166]73]. In
   terms of genetic factors, this study explored the genetic diversity of
   different habitats of P. tabuliformis and carried out genetic distance
   analysis. The results showed that although the 3 populations are
   closely related to each other, there is still a certain genetic
   distance between the old trees and mature trees (Fig. [167]4b, c),
   which may also be an important reason for the variation in leaf
   functional traits. Compared with mature trees, old trees have greater
   genetic diversity (Table [168]2), possibly because the gene pool of old
   trees is a product of environmental heterogeneity over a long period of
   evolution, possibly leading to the survival of individuals with high
   genetic diversity [[169]74]. Moreover, the degree of genetic diversity
   influences a species’ adaptability to environmental changes, with
   greater genetic variation correlating with enhanced environmental
   adaptability [[170]75]. This correlation also suggests that old trees
   possess robust environmental adaptability. The VPA results revealed
   that the common explanation rate of soil property factors and genetic
   factors was 67.89%, and the single explanation rate of genetic factors
   was higher than that of soil property factors (Fig. [171]4d). These
   findings suggest that the synergistic effect on the leaf functional
   traits of P. tabuliformis cannot be ignored. Similar findings were also
   found in other related studies, for example, the stomatal traits of
   Betula ermanii were regulated by both environmental and genetic
   factors, and stomatal size was more regulated by genetics [[172]76];
   the leaf morphological traits of old Ginkgo biloba L. are also more
   influenced by genetic factors, but environmental factors can induce
   morphological variation by affecting the accumulation of important
   secondary metabolites [[173]77]. However, the specific synergistic
   mechanisms need to be investigated in future studies. Furthermore,
   human-mediated interference is also an important factor affecting plant
   traits, but it is difficult to predict such interference because it is
   usually the co-occurrence and compound effect of multiple interference
   factors [[174]78], which also leads to the limited interpretation of
   leaf functional trait variation in our study.

Molecular mechanism of leaf functional trait variation

   The transcriptome furnishes the molecular foundation for the analysis
   of gene expression across biological processes. By conducting
   transcriptome and WGCNA analyses of P. tabuliformis leaves from 3
   populations, 3 important modules and 24 candidate hub genes were
   selected (Figs. [175]5 and [176]7). The results of the KEGG enrichment
   module analysis indicated that the significantly enriched pathways
   affected mainly material and energy metabolism and distribution, signal
   transduction, stress response and other life activities (Fig. [177]6),
   which were closely related to the construction of leaf functional
   traits of P. tabuliformis. For example, the carotenoid biosynthesis
   (MEsalmon) pathway plays a variety of roles related to photosynthesis,
   photoprotection, development and stress hormones in plants [[178]79].
   Park et al. found that overexpression of ZEP in Arabidopsis thaliana
   enhanced plant tolerance to osmotic stress [[179]80]. In this study,
   the expression of ZEP (Pt8G30020) gene in old trees was higher than
   that in mature trees, which was similar to the results of functional
   traits. In the arid area of the Loess Plateau, the high expression of
   ZEP gene may regulate the leaf physiological traits of old trees to
   better adapt to environmental changes. The proteasome (MEblue) pathway
   also plays an important role in plant growth and development, although
   there is no direct evidence that the proteasome is related to leaf
   functional traits, studies in Arabidopsis and soybean have found that
   it may be involved in regulating the level of cytokinin in plants
   [[180]81], there is a significant correlation between cytokinin and
   leaf traits formation and senescence [[181]82, [182]83]. Both Pt9G38890
   and Pt7G08750 genes belong to RPN12 (26 S proteasome regulatory subunit
   N12). Although the function of RPN12 in P. tabuliformis has not yet
   been verified, in Arabidopsis, the leaf growth rate and morphology of
   Rpn12a-1 mutant changed [[183]84], so the significant enrichment of
   RPN12 in P. tabuliformis may be one of the factors inducing the change
   of leaf traits. Subsequently, the annotation analysis of hub genes was
   carried out in this study, and the functions of some hub genes were
   been verified and determined to be related to plant growth and
   development and abiotic stress. Auxin plays an important role in plant
   leaf formation [[184]85]. N-MYC downregulation-like 3 (NDL3)
   (Pt2G50170) and Ubiquitin carboxyl-terminal hydrolase family protein
   (UCH) (Pt5G35270) have been found to be associated with the auxin level
   of Arabidopsis [[185]86, [186]87]. In this study, Pt2G50170 and
   Pt5G35270 had opposite effects on leaf functional traits
   (Fig. [187]S5). However, the high expression of Pt5G35270 in old trees
   is positively correlated with the leaf economic, anatomical and
   physiological traits, which may be due to the synergistic effect of UCH
   and various auxin [[188]88], which indirectly affects the variation of
   functional traits. KCS (Pt7G22770) is a key rate-limiting enzyme in
   plant wax synthesis and actively participates in physiological and
   biochemical reactions at various stages of plant growth and
   development. For example, Arabidopsis HIC gene is involved in the
   regulation of stomatal number [[189]89], and the expression pattern of
   KCS gene in barley under drought stress is ultimately manifested as
   changes in the wax content and composition of Hordeum vulgare L. skin
   [[190]90]. In this study, it was found that the expression of KCS gene
   in old trees was significantly lower than that in mature trees, and
   there was a significant positive correlation with the variation of
   photosynthetic traits, indicating that KCS gene may affect
   photosynthetic traits through chloroplast development [[191]91].
   CYP86B1 (Pt7G00280) belongs to fatty acid ω-hydroxylase, it localized
   in chloroplast outer envelope [[192]92]. In this study, CYP86B1 gene
   expression and photosynthetic trait value showed a uniform trend (old
   trees are significantly smaller than mature trees). These results
   suggest that cork lipoprotein biosynthesis may lead to changes in leaf
   cell wall function and thus affect leaf photosynthetic traits
   [[193]93]. In addition, some of the secondary metabolites associated
   with growth within the cell cannot be transported through vesicles, and
   the presence of lipid-transfer protein (Pt5G28900), ATP-binding
   cassette B2 (PtXG14450) and proton pump interactor 1 (Pt1G44120) allows
   these substances to be widely distributed [[194]94–[195]96]. Other Hub
   genes whose functions are still unclear may also play significant roles
   in the variation in leaf functional traits. Further experimental
   verification is needed. In conclusion, leaf functional traits need to
   be regulated through complex and highly coordinated molecular
   processes, and the screening of hub genes helps us to obtain a
   preliminary understanding of these processes. Further studies need to
   be conducted on these candidate hub genes to elucidate the molecular
   mechanisms underlying their regulation and to harness their full
   potential.

   Our study shows that there are differences in leaf functional traits
   between old and mature trees in the Loess Plateau, W and TP are the
   main driving factors of leaf functional traits variation, and the
   possible regulatory genes of functional traits are discussed. Based on
   the above research results, in the protection of old trees, we can
   regulate the variation of leaf functional traits of old trees by
   controlling the soil microenvironment. Reasonable nutrient addition and
   water control can effectively promote the change of survival strategies
   of old trees and better enhance their adaptability to the environment.
   The study of genes related to functional traits is helpful to enhance
   our understanding of the mechanism of trait variation, and provide
   genome resources for developing adaptive strategies to cope with the
   environment of P. tabuliformis while formulating scientific and
   rational conservation strategies.

Conclusion

   In this study, the leaf functional traits (economic, anatomical,
   photosynthetic and physiological traits) of P. tabuliformis were used
   as indicators to evaluate the differences in survival strategies
   between old and mature trees. In contrast, the resource acquisition
   strategy of old trees was more conservative, and the all leaf traits
   changed with changes in survival strategy. This transformation is
   influenced by both soil property and genetic factors, with genetic
   factors exerting a more pronounced effect; however, the nature of their
   interaction remains ambiguous. In addition, constructed a complex
   regulatory network by synthesizing transcriptome data, and identified
   candidate hub genes related to leaf functional trait variation. The
   findings of this study are helpful to further understand the survival
   mechanism of old P. tabuliformis in the Loess Plateau, and provide
   scientific basis for formulating effective conservation strategies.
   Subsequent research should delve deeper into the mode of action of the
   driver and function of the candidate hub genes to better understand the
   regulatory mechanism of functional variation in the leaves of P.
   tabuliformis.

Electronic supplementary material

   Below is the link to the electronic supplementary material.
   [196]Supplementary Material 1^ (1.1MB, docx)

Acknowledgements