The experiment was conducted at the Experimental Forest Farm in Qingyuan County, Zhejiang Province, China (119°01’25″ E, 27°38’48″ N), with an altitude ranging from 480 m to 510 m. The general situation of the research area is shown in Figure 1 . The experimental site has a gentle slope. The region has an average annual temperature of 17.6 °C, an annual precipitation of 1721.3 mm, and a frost-free period of 245 days, characteristic of a typical subtropical monsoon climate. The experimental site comprises acidic red soil characterized by a pH of 4.67. Key soil properties include 7.8 g/kg organic matter, 0.45 g/kg total nitrogen, 0.27 g/kg total phosphorus, and available nutrient concentrations of 17.2 mg/kg nitrogen, 37.3 mg/kg potassium, and 1.02 mg/kg phosphorus.

A completely randomized block design with two factors (variety and fertilization) was used. Three P. massoniana hybrid varieties were selected: F32 (♀6627×♂1003), F35 (♀5163×♂1003), and F36 (♀5163×♂6627). These hybrid families were generated through artificial controlled pollination within a second-generation breeding population, with all three varieties sharing either a common paternal or maternal parent. The selected genotypes represented distinct geographic provenances: (1) genotype 6627, derived from Jiangxi provenance (central production region), exhibiting 2–3 annual shoot flushes; (2) genotype 1003, originating from Guangdong provenance (southernmost production region), characterized by 3–5 annual shoot flushes; and (3) genotype 5163, representing Zhejiang provenance (northern production region), typically showing 1–2 annual shoot flushes. Four phosphorus fertilizer treatments were established: 0 g·plant^-¹ (T0, as control), 100 g·plant^-¹ (T1), 250 g·plant^-¹ (T2), and 500 g·plant^-¹ (T3). The applied phosphorus fertilizer was diammonium phosphate ((NH₄)₂HPO₄: containing 46% phosphorus and 18% nitrogen). In March 2020, a measurement forest was established according to the scheme shown in Figure 2 . The 12 treatment combinations were replicated five times (blocks), with each block containing all treatments randomly arranged. Each treatment plot comprised 40 trees (4 rows × 10 plants, 2 m × 2 m spacing), totaling 2,400 trees (12 treatments × 5 blocks × 40 trees). Fertilization was applied in May 2020–2022 via a continuous circular ditch (15 cm from trunk, 20 cm deep, 10 cm wide) around each tree. Growth traits (fixed growth: predetermined bud flushes; free growth: post-flush shoot elongation) were assessed for all 40 trees per plot (no subsampling). Protective borders of open-pollinated seedlings surrounded each block. Annual weeding ensured minimal competition, with no additional nutrient/water inputs.

After the third growing season, tree height, diameter at breast height (DBH), and annual shoot characteristics were measured for each tree. The measurement of annual shoot characteristics was based on the whorled growth pattern of P. massoniana. After overwintering, the previous year’s winter buds began to sprout in the following spring. Shoot growth was determined by measuring the length between two spring branch whorls and the number of branch whorls, thereby calculating the annual height growth and shoot number. The growth of the initial shoot (fixed growth) was the distance from the base of the spring branch whorl to the base of the first shoot. For free whorl branches, the measurement was carried out upward in the same way to obtain the annual shoot number and growth. In addition, four individual trees were randomly selected from each treatment combination in each replicate, and the current year’s needles from the upper part of the tree crown were collected and brought back to the laboratory. The needles were deactivated at 105 °C for 30 minutes, then dried at 80°C to constant weight, and their mass was measured for the determination of nitrogen and phosphorus content. Phosphorus content was determined by the molybdenum-antimony colorimetric method, and nitrogen content was determined by the Kjeldahl method (Bremner, 1960; Zhang et al., 2025).

We conducted comprehensive analyses of growth traits, shoot morphology, and needle nutrient content in P. massoniana saplings across hybrid varieties and phosphorus (P) fertilization treatments. The descriptive statistics (mean, range, standard deviation) and coefficients of variation were calculated for all measured parameters. The calculation formulas are as follows:

Where x i is the measured value for each individual plant, x ¯ is the mean value for a trait, n is the number of individual plants, δ G and δ E represent genetic and environmental variance components respectively, estimated using restricted maximum likelihood (REML) methods in a mixed-effects model framework.

The Shapiro-Wilk test was employed to assess the normal distribution of measured data across different hybrid varieties and phosphorus fertilizer application rates. One-way analysis of variance (ANOVA) and Tukey’s multiple comparison test were utilized to analyze the effects of varying phosphorus fertilizer intensities and varieties on needle nitrogen and phosphorus content. Partial least squares (PLS) regression was employed to examine multivariate relationships between functional traits and nutrient content across P gradients, addressing potential collinearity among predictors. All analyses were implemented in R 4.3.2. Effect sizes (η² for ANOVA, r for correlations) are reported alongside p-values.

After three years of cultivation, the P. massoniana exhibited an average tree height of 2.68 m, a diameter at breast height (DBH) of 24.22 mm, 1.12 flushes per year, a fixed growth of 85.02 cm, and a free growth of 42.39 cm. The free growth trait showed the highest genetic coefficient of variation (GCV=11.54%) and environmental coefficient of variation (ECV=52.99%), indicating that it has substantial plasticity. The growth and flushing characteristics of the hybrid varieties responded differently to phosphorus fertilizer application. Significant differences were observed among varieties and phosphorus treatments in tree height, DBH, and fixed growth ( Table 1 , Figure 3A ). The T3 treatment showed the best performance in tree height, DBH, and fixed growth, with increases of 6.14%, 5.49%, and 14.35%, respectively, compared to the T0 treatment ( Figure 4 ). This indicates that higher phosphorus fertilizer content significantly promoted the development of primary growth morphology in P. massoniana. No significant differences were found in free growth and annual flushing frequency among varieties and phosphorus treatments, suggesting that the multiple flushing behavior of hybrid varieties was not highly sensitive to phosphorus fertilizer.

The genetic coefficients of variation for needle nitrogen and phosphorus content and the nitrogen-to-phosphorus ratio were relatively small (1.03%-3.80%), lower than the environmental coefficients of variation (10.85%-19.75%) ( Table 1 , Figure 3A ). No significant differences were observed in needle nutrient content among varieties, but significant differences were found among fertilizer treatments (P<0.05). With increasing fertilizer intensity, needle nitrogen and phosphorus content and the nitrogen-to-phosphorus ratio showed an increasing trend ( Figure 4 ), indicates that increased external phosphorus supply promoted uptake of both nitrogen and phosphorus in P. massoniana needle.

The distribution of all indicators varied across three dimensions: variety (V), fertilization (P), and their interaction (V×P) ( Figure 3B ). The effect on DBH was in the order of variety effect > interaction effect (V×P) > fertilization (P) effect, with the variety effect contributing 51.22%. The effects on tree height, fixed growth, needle nitrogen content, phosphorus content, and nitrogen-to-phosphorus ratio were in the order of fertilization (P) effect > interaction effect (V×P) > variety effect. The contribution of the fertilization (P) effect ranged from 55.43% to 93.62%, and the contribution of the interaction effect (V×P) ranged from 3.13% to 22.42%. This indicates that fertilization is the dominant factor influencing initial growth and nutrient accumulation, while also being affected by the interaction effect. The interaction effect between variety and fertilization (V×P) had the highest contribution to free growth, reaching 74.1%, suggesting that the flushing behavior of different varieties exhibits heterogeneity in response to fertilization intensity.

The fixed growth of hybrid varieties 32 and 35 increased under fertilization treatments, peaking at T1 and T3 treatments, respectively; in contrast, clone 36 showed the opposite trend. The effects of different treatments on the free growth of each variety varied considerably ( Figure 5 ). This indicates that the increase in annual flushing length is attributed to the accumulation of fixed growth and free growth. Hybrid combinations with the same paternal or maternal parent exhibited different response patterns to phosphorus input. However, regardless of low or high phosphorus fertilizer input, the increase in annual flushing length of high-generation hybrid varieties was mainly attributed to the increase in fixed growth.

As the amount of phosphorus fertilizer application increased, more correlations were observed between tree height, DBH, and nutrient factors ( Figure 6 ). Flushing frequency showed mostly significant negative correlations with other traits (P< 0.05), especially under T1 and T2 treatments. The N:P ratio was positively correlated with tree height and DBH, but negatively correlated with other traits. Furthermore, under the T4 treatment, with increasing fertilization, nitrogen, phosphorus, and the N:P ratio showed positive correlations with multiple growth indicators, particularly with tree height (P< 0.01). This indicates that the nutritional environment influences the genetic expression of traits, thereby affecting the interrelationships between traits.

Phosphorus fertilizer application significantly promoted the growth of P. massoniana young forests, notably improving traits such as tree height, DBH, and fixed growth ( Figure 3 ). As a key nutrient, phosphorus input also increased nitrogen and phosphorus contents in needles, especially phosphorus uptake ( Figure 4 ), directly enhancing biomass accumulation and morphogenesis during early growth (Na et al., 2014; Yi et al., 2023). Additionally, phosphorus input indirectly improved nitrogen use efficiency by enhancing soil nutrients or root activity (Mulligan and Sands, 1988). Growth of P. massoniana is limited by nitrogen and phosphorus deficiency in southern acidic soils, particularly phosphorus (Xu et al., 2018). This study confirms that exogenous phosphorus alleviates this limitation, improves soil nutrients and plant uptake, and significantly boosts P. massoniana growth, especially under low soil phosphorus conditions (Xu et al., 2018). Research has shown that fertilization’s effect on needle nutrients is mainly driven by environmental factors, with fertilization treatments accounting for 52.99%–64.08% of the variation in nitrogen and phosphorus content ( Table 1 ), while genotypic influence is relatively limited. This supports Guo W. et al. (2023), who suggested needle nutrient content reflects environmental adaptation more than genotype. Under significant environmental changes, plants exhibit high plasticity in nutrient uptake and allocation, highlighting the importance of fertilization management in forestry (Rubio and Lavado, 1999; Andersen et al., 2017). Fertilization supports optimal growth in variable environments (Freiberger et al., 2013; Niu et al., 2020; Rennenberg et al., 2009).

Under low phosphorus conditions, a negative correlation was observed between DBH and needle phosphorus content ( Figure 5 ); as phosphorus supply increased, this shifted to a significant positive correlation, reflecting P. massoniana nutrient allocation strategy under phosphorus stress. When phosphorus is limited, the species may prioritize allocation to growth organs (e.g., roots, stems), reducing leaf phosphorus to sustain growth (Turner, 2008). This adaptive mechanism allows plants to maintain growth under stress (DeLoose et al., 2024). By analyzing needle nitrogen and phosphorus content under different phosphorus levels and their correlations with traits like flushing growth and frequency, the study found that nutrient environment strongly influences trait expression and interactions. Without phosphorus fertilizer, traits such as tree height and DBH showed strong positive correlations with other growth traits, which weakened with increased phosphorus input. Thus, appropriate phosphorus fertilization not only alters needle nutrient levels but also indirectly affects flushing traits via the nutritional environment, deepening our understanding of P. massoniana growth mechanisms under phosphorus variation.

This study revealed significant differences in the growth responses of P. massoniana hybrid varieties to phosphorus fertilizer intensity, highlighting the critical role of genotype–environment interaction (V×P) in shaping young forest traits. Hybrid combination 32 performed best under moderate phosphorus treatment (T2), with significantly greater tree height and DBH, while combinations 35 and 36 showed better growth under high (T3) and low (T1) phosphorus levels, respectively ( Figure 4 ). This demonstrates genotypic variation in phosphorus response and the influence of fertilization intensity on growth potential (Debell et al., 1986). Combination 32, derived from southernmost and central provenances, exhibited high growth potential, while combinations 35 and 36, with northern provenance genes, showed different responses—indicating that growth performance is linked to parent origin. These findings align with Stovall et al. (2011), who reported that genotype × fertilization interaction significantly influenced height and stem growth in slash pine. Other studies confirm that environmental factors, especially fertilization intensity, not only directly affect growth but also interact with genotype responses, emphasizing the importance of G×E effects in forest development (Keller et al., 2018; Zhao et al., 2016).

Tree height showed higher heritability (h² = 0.37) than free growth (h² = 0.30) ( Table 1 ), suggesting stronger genetic control over apical growth, while free growth was more influenced by environmental factors, particularly fertilization intensity. Its high plasticity (CV = 54.08%) reflects sensitivity to phosphorus availability, with trees under low phosphorus prioritizing structural growth over flushing—evidenced by the negative DBH–needle P correlation. This resource allocation is hormonally regulated via gibberellin–auxin interactions and elevated ABA under phosphorus stress (Liu et al., 2024). These patterns align with the “source–sink redistribution” model (Suriyagoda et al., 2023), highlighting the need to consider fertilization effects on plastic traits like free growth when managing P. massoniana. Genetic variation in flushing traits (h² = 0.30) suggests differential phosphorus sensitivity among clones, shaped by provenance adaptation. Free growth was most affected by variety × fertilization interaction (74.1%), while annual flushing frequency showed low responsiveness, indicating it is less sensitive to phosphorus input. Overall, fertilization interacts with genotype to regulate growth plasticity in P. massoniana.

This study found that hybrid combination 32 performed best under moderate phosphorus levels, while combinations 35 and 36 showed better growth under high and low phosphorus, respectively. This suggests that phosphorus concentration has genotype-specific effects on P. massoniana growth, with each hybrid maximizing its potential under optimal conditions (Conroy et al., 1990). Such differences likely arise from genotypic variation in nutrient absorption, utilization, and allocation (Liu et al., 2023; McKeand et al., 2023). Therefore, selecting genotypes suited to specific phosphorus levels can enhance fertilization efficiency and forest productivity (Tie et al., 2024; Yang et al., 2025).

Shoot-sprouting behavior can reflect the growth rate of P. massoniana, with fixed growth and free growth being closely related to genetic and environmental factors, respectively. For hybrid varieties with high-generation genetic improvement, the variation in flushing frequency may not be obvious, but mainly manifested as changes in fixed growth. This means that during the genetic improvement process, the genetic stability of varieties is enhanced, and the buffering capacity against environmental changes is improved, making the flushing frequency relatively stable. Free growth, on the other hand, is more influenced by environmental stimuli. Under different phosphorus fertilizer intensities, the variation in free growth among different clones varied greatly and showed no obvious pattern, reflecting the complex role of environmental factors on it. In forestry practice, this study provides a theoretical basis for “precision genotype-environment matching”. For example, hybrid combination 32 is preferentially planted on sites with moderate phosphorus availability to maximize its root efficiency and phosphorus absorption potential; while hybrid combination 35 is promoted in high-phosphorus soils to avoid phosphorus toxicity risks and optimize resource utilization. This is consistent with the “suitable genotype for suitable site” strategy proposed by Chen et al. (2016), which emphasizes improving plantation productivity and sustainability through precise matching of genotype and environment.

In addition, in actual forestry production, rational fertilization is one of the effective ways to improve the nutrient status and promote the growth of P. massoniana (Mao et al., 2021; Guo Q. et al., 2023). More attention should be paid to the regulation of DBH growth, and the timber yield of P. massoniana should be increased by optimizing fertilization and selecting suitable genotypes (Qiao et al., 2022). In the future, genomic selection technology should be combined to screen families with high phosphorus response potential, which can accelerate the genetic improvement process of P. massoniana. Future research can combine genomics approaches (such as transcriptomics and metabolomics) to analyze the molecular mechanisms of nutrient response in P. massoniana, and monitor the feedback effects of fertilization on soil fertility in the long term, to achieve sustainable management of plantations.