The field experiment was conducted from March to July 2024 in Hexi Village, Guanshi Town, Hengnan County, Hengyang City, Hunan Province, China (26°43′12″N, 112°50′55″E; altitude 70 m), The region has a subtropical monsoon climate, with an average temperature of approximately 24 °C and precipitation of 800 mm during the tobacco growing season (March-July). The flue-cured tobacco variety used in this study was ‘Yunyan 87’. The site, previously cultivated with rice, had the following basic soil properties: pH 7.46, total nitrogen 2.47 g·kg^-¹, total potassium 1.69 g·kg^-¹, total phosphorus 0.95 g·kg^-¹, alkali-hydrolyzable nitrogen 90.44 mg·kg^-¹, available potassium 379.90 mg·kg^-¹, available phosphorus 49.63 mg·kg^-¹, and organic matter 43.16 g·kg^-¹.

The experiment was arranged in a randomized complete block design with three replicates. Under a constant total N application rate (165 kg·N·ha^-¹), five treatments with varying proportions of organic N were established: ON0 (0% organic N), ON10 (10%), ON20 (20%), ON30 (30%), and ON40 (40%), totaling 15 plots. Guard rows were established around the experimental area, and the planting density was maintained at 17,250 plants ha^-¹. The fertilizers used included a specialized basal fertilizer for tobacco (8% N, 10% P₂O₅, 11% K₂O), a starter fertilizer (20% N, 9% P₂O₅), a specialized topdressing fertilizer for tobacco (10% N, 32% K₂O), potassium sulfate (52% K₂O), potassium magnesium sulfate (25% K₂O), and superphosphate (12% P₂O₅). The organic N source was a bio-fermented rapeseed cake fertilizer (2.4% TN, 1.05% TP, 2.55% TK). All other agronomic practices followed the 2024 technical guidelines provided by the Hengyang City Company of the Hunan Provincial Tobacco Monopoly Bureau. Detailed fertilization data for each treatment are provided in Supplementary Table 1.

To ensure consistency, samples were collected at two distinct stages targeting specific physiological and quality indicators. First, 90 days after transplantation, fresh leaf samples were collected from the sixth leaf position (starting from the top), and 10 leaves were collected from each replicate. After removal of the veins, the leaves were stored at -80°C for determination of physiological and biochemical indexes, including enzyme activities (FAS, PEPC, PAP) and metabolite contents (PA, TP). Second, after the flue-curing process, cured leaves of the B2F grade (upper leaves) were selected to measure chemical components (PEE, TS, RS), physical properties (softness, tensile strength, thickness), and the sensory oil content score.

The contents of pyruvic acid (PA), total phenols (TP), and the activities of FAS, PEPC, and PAP were measured using reagent kits from Suzhou Grace Biotechnology Co., Ltd. (Suzhou, China). The mass of all measured samples was 0.1 g, and the volume of solvent was 1 mL.

The content of PA was determined based on its reaction with 2,4-dinitrophenylhydrazine to form 2,4-dinitrophenylhydrazone, a brownish-red hydrazone in an alkaline solution. Absorbance was measured at 520 nm.

The TP content was measured using the Folin-Ciocalteu method, based on phenol-induced reduction of tungstomolybdic acid, producing a blue complex (760 nm).

FAS activity was monitoring NADPH consumption (absorbance at 340 nm) during conversion of acetyl-CoA and malonyl-CoA to long-chain fatty acids.

PEPC activity was assessed through coupled reactions with malate dehydrogenase, tracking NADH oxidation at 340 nm. PAP activity was determined via β-glycerophosphate hydrolysis, detecting liberated inorganic phosphate using a phosphorus-specific reagent.

All assays were conducted in 96-well microtiter plates, with triplicate technical replicates.

The contents of total sugar (TS) and reducing sugar (RS) were determined following YC/T 159-2002, and petroleum ether extract (PEE) was measured according to YC/T 176-2003.

During the germination period, the upper surface of the first mature leaf (from top to bottom) was taken, stained with 0.2% (w/v) Rhodamine B for 30 minutes, rinsed three times with distilled water, and blotted dry. Glandular trichome density was quantified using a digital microscope (VHX-500F, KEYENCE Corporation, Japan) by analyzing three random fields from the central upper epidermis.

The first mature leaf (from top to bottom) was taken during the germination period. Leaf discs with a diameter of 5 cm were cut from the leaves. Fifty leaf discs constituted one biological replicate, and three biological replicates were prepared for each material. The leaf discs were extracted by dipping them eight times (2 s each) into a dichloromethane solution. A 1 mL internal standard solution (a mixture of 2.020 mg·ml^-¹ sucrose octaacetate and 2.542 mg·ml^-¹ n-heptadecanol) was added to the extract. After thorough mixing, the solution was filtered and then concentrated using a rotary evaporator (LCA-RN-1300BE, Shanghai Lichen Instrument Technology Co., Ltd.). The concentrate was dried under a nitrogen stream using a nitrogen evaporator (LC-DCY-12GP, Shanghai Lichen Instrument Technology Co., Ltd.). The dried residue was subsequently silylated and analyzed by GC-MS (7890B-5977A, Agilent Technologies Inc., USA).

The leaf softness value (mN) was measured longitudinally and transversely using a softness tester (YT-RRY1000, Hangzhou Yante Technology Co., Ltd.). The tensile strength was determined with a force gauge (SF-10, Wenzhou Weidu Electronics Co., Ltd.), and the thickness was determined using a high-precision thickness gauge (547-401, Mitutoyo Corporation, Japan).

Oil content was scored based on GB 2635-92 (“Grading of Flue-cured Tobacco”), supplemented by expert input: abundant (9–10 points), present (6–8 points), slightly present (3–5 points), and sparse (0–2 points). Evaluators were blinded to the specific treatments to minimize bias. The scores were treated as a continuous variable for regression modeling, a common practice in sensory analysis to capture subtle quality gradients.

The model was developed using the RandomForestRegressor class from the scikit-learn library in Python. The standardized key indicators were used as input features, and the original oil content scores served as the target variable. To overcome the limitations of a small sample size (N = 15 biological replicates) and to capture robust non-linear relationships, we implemented a comprehensive modeling framework involving data augmentation, recursive feature elimination (RFE), and hyperparameter optimization.

Data were processed using Excel 2021, SPSS 19.0, and Pycharm 2024. Line graphs and bar charts were produced with OriginPro 2024. One-way analysis of variance (ANOVA) followed by Duncan’s new multiple range test was used for variance analysis. The Spearman method was used to calculate correlation coefficients.

During the budding stage, the morphology and density of glandular trichomes on the upper leaves were observed and quantified for each treatment. Figure 1 shows the density of segmented long-stalked glandular hairs gradually increases from ON0 to ON30, but significantly decreases at ON40. The results showed that glandular trichome density followed the order of ON30 > ON20 > ON10 > ON40 > ON0. Specifically, ON30 exhibited a 50.98% increase in total trichome density compared to ON0 (Figure 2A). The density pattern of long-stalked glandular trichomes, key sites for secretion, mirrored that of total trichomes, with ON30 showing a 53.24% increase compared to ON0. Short-stalked trichome density was highest in ON20; while differences among organic nitrogen treatments were not significant, all exceeded ON0 significantly. Non-secretory epidermal trichomes showed no treatment differences.

GC–MS analysis of trichome secretions revealed that the contents of all measured compounds in ON20, ON30, and ON40 were significantly higher than in ON0 (Figure 2B). The largest increase occurred in ON30, where cembratriene-diol (CBT-diol) content rose by 68.13% and sucrose esters by 47.39% compared to ON0, paralleling the elevated density of the long-stalked trichomes.

Organic nitrogen ratio significantly influenced lipid-related enzyme activities (Table 1). FAS activity peaked in ON30, showing a 27.8% increase over ON0 and exceeded ON10 and ON20 significantly. Although slightly lower in ON40, FAS activity remained significantly above ON0. PAP activity significantly increased only in ON30 (by 44.2% compared to NO0), with other treatments showing non-significant increases. Similarly, PEPC activity also peaked in ON30, rising 82.3% over ON0, and was significantly higher than in all other treatments. Taken together, ON30 (a moderately high level of organic nitrogen) maximized FAS, PAP, and PEPC activities.

Different lowercase letters following the data in the same column indicate significant differences among treatments at the 0.05 level. The same notation applies to subsequent tables. FAS: Fatty acid synthase, PAP: Phosphatidic acid phosphatase, PEPC: Phosphoenolpyruvate carboxylase.

The contents of PEE, PA, TP, TS, and RS were analyzed in cured leaves (Table 2). PEE content was highest in ON30, showing an 18.44% increase compared to Ono, while other treatments did not differ significantly from ON0. PA content decreased significantly in ON30 and ON40 (-17.63% and -5.93%, respectively), whereas ON10 showed a slight increase over ON0. TP content decreased significantly across all treatments compared to ON0, likely reflecting enhanced pyruvate metabolism and downstream diterpenoid biosynthesis, suppressing the phenylpropanoid pathway. Both TS and RS increased with organic nitrogen up to ON30, where they peaked (+43.19% and +64.31% vs. ON0, respectively), before declining in ON40.

Leaf softness, tensile strength, and thickness responded non-linearly to organic nitrogen (Table 3). Softness value (inversely related to softness) decreased then increased with rising organic nitrogen, reaching its lowest (softest) point in ON30, significantly below all other treatments. Tensile strength and leaf thickness followed a rising–falling trend, both peaking in ON30. All treatments except ON0 produced leaves within the optimal thickness range (120–140 μm) for upper leaves, with ON30 yielding the thickest leaves.

Oil content scores of cured upper leaves displayed a non-linear response to organic nitrogen (Figure 3). ON30 and ON40 achieved the highest scores, both significantly exceeding ON0.

Spearman correlation analysis (Figure 4) was conducted to explore the relationships among oil content and physicochemical indicators. As shown in the heatmap, oil content exhibited significant positive correlations (p< 0.05) with PEE, FAS, TS, RS, Glandular trichomes density (GT), CBT-diol, PEPC, PAP, and tensile strength (TsS). Conversely, it showed significant negative correlations with PA content and softness value (Sof). Notably, strong multicollinearity was observed between certain features, such as TS and RS, necessitating feature selection for subsequent modeling.

To construct a robust predictive model despite the limited sample size, we employed a data augmentation strategy combined with rigorous feature selection, and encoded the treatments as 0, 10, 20, 30, 40.

Recursive Feature Elimination (RFE) with cross-validation was first applied to the augmented dataset. As shown in Figures 5A, the model’s performance stabilized when the number of features was reduced to five. Consequently, PEE, CBT-diol, Sof, RS, and GT were selected as the optimal input features for the final model, effectively reducing redundancy while retaining critical biological information.

The Random Forest model trained on this optimized feature set and augmented data demonstrated excellent performance. Grid search optimization (Figure 5E) identified the optimal hyperparameters (e.g., n_estimators=300, max_depth=10, min_samples_leaf=1), which were used for the final training. The model achieved an OOB R² of 0.781 and a 3-fold CV R² of 0.819. In stark contrast, the model trained on the original small dataset (N = 15) failed to generalize, yielding a negative CV R² (-0.039) (Figures 5C, D). This comparison underscores the necessity and effectiveness of our data augmentation strategy. The OOB error curve (Figure 5B) further confirmed that the error rate stabilized and minimized around 150–200 trees, validating the sufficiency of setting n_estimators to 300.

Model diagnostics indicated that the assumptions for regression were met. The distribution of residuals was approximately normal (Figure 5F), and the Q-Q plot (Figure 5G) showed that the standardized residuals aligned well with the theoretical quantiles, with no significant deviations. Sensitivity analysis (Figure 5H) demonstrated that the model’s performance remained relatively stable across different noise levels (1% - 3%) used in augmentation, although performance degraded at very high noise levels (10%), confirming the appropriateness of the chosen 3% noise parameter.

Feature importance analysis based on Mean Decrease in Impurity (MDI) (Figure 6) revealed the relative contribution of each indicator to the oil content prediction. PEE was identified as the most influential predictor, followed closely by CBT-diol and Sof. The feature importance rankings between the trained models based on the raw data and enhanced data are almost consistent, further verifying the robustness of identifying key features.

To gain a deeper understanding of how key indicators influence oil content, we interpreted the RF model using Partial Dependence Plots (PDP) and SHAP value analysis.

The PDP (Figure 7) revealed the marginal effects of the top five features. PEE and CBT-diol demonstrated strong positive associations with oil content: as their values increased, the predicted oil content score rose sharply before leveling off, suggesting a saturation effect. GT also showed a positive trend. In contrast, Sof exhibited a clear negative relationship, where lower softness values (indicating softer leaves) were associated with higher oil scores. PA showed a threshold effect, where oil scores dropped significantly when PA content exceeded a certain level.

SHAP analysis provided further insights into individual feature contributions and interactions. The SHAP summary plot (Figure 8A) confirmed the global importance trends, showing that high values of PEE, CBT-diol, and GT had a positive impact on the model output (pushing it to the right), while high values of Sof and PA had a negative impact. We further explored the interactions among the top four features using SHAP dependence plots (Figures 8B–E), which revealed that the influence of individual indicators is highly context-dependent. For instance, the positive contribution of PEE to oil content was notably modulated by CBT-diol, specifically, samples characterized by both high PEE and high CBT-diol levels exhibited elevated SHAP values, suggesting a synergistic effect between these lipid-related metabolites. A similar but negative, more complex, interaction was observed between CBT-diol and leaf softness. Differently, the negative effects of Sof and PA on oil content were found to vary depending on the levels of PEE and Sof, respectively. These results underscore that tobacco oil content is not determined by single factors in isolation, but rather by the complex interplay among lipid metabolites, physical properties, and precursor availability.

t-SNE visualization (Figure 9) was employed to explore the underlying data structure. Figure 9A shows that samples from different organic nitrogen treatments (ON0-ON40) formed distinct clusters or gradients in the 2D space, confirming that the selected five features contain sufficient information to distinguish between treatments. Figure 9B displays the features in the t-SNE space, colored by their importance. The proximity of certain features (e.g., GT and CBT-diol) suggests they may share similar information or biological roles in determining leaf quality.

Finally, we compared the treatment levels with the RF model’s cross-validation predictions (Figure 10). The average predicted oil content for each treatment level accurately reproduced the non-linear trend observed in the actual field experiment: rising from ON0 to a peak at ON30, and then slightly declining at ON40. This confirms that our machine learning model successfully captured the true biological response of tobacco oil content to varying organic nitrogen application ratios.

The secretory trichomes (glandular trichomes) on the leaf surface of flue-cured tobacco at specific developmental stages are a key factor influencing oil content and overall quality. Glandular trichomes are the sites for the synthesis and secretion of numerous very-long-chain fatty acids, their derivatives, and alicyclic compounds (Philippe et al., 2022; Wang, 2014). Among these, cembratriene-diol and sucrose esters are important secretions of the glandular trichomes. Both are crucial aroma precursors that can impart desirable aroma and sweetness to flue-cured tobacco (Qu et al., 2018), and they enhance the “oily” sensation of the leaf surface, contributing to the improvement of the tobacco’s appearance quality. As the primary secretory organs for terpenoids, an increased density of long-stalked glandular trichomes directly thickens the waxy layer on the leaf surface, thereby enhancing the tactile perception of oiliness (e.g., smoothness) and the aroma (Wang et al., 2022). In this study, all treatments with increased organic nitrogen ratios showed significantly greater densities of all types of glandular trichomes and higher secretion contents compared to the control. Notably, in ON30, the density of long-stalked glandular trichomes increased by 53.24% compared to ON0, with a concurrent significant rise in the contents of CBT-diol and sucrose ester secretions. This is likely the reason why the final oil content score for ON30 was significantly higher than that for ON0. Our Random Forest model further validated this biological observation: Recursive Feature Elimination (RFE) identified both Glandular Trichome density (GT) and CBT-diol content as top-tier predictors (Figure 5A), confirming their indispensable role in determining the final oil score. Some researchers suggest that the content of glandular trichome secretions may be related to the regulation of cell wall structure by the combined application of organic nitrogen (Zhou et al., 2024). For instance, the increased leaf thickness observed in ON30 might enhance the mechanical strength of the leaf, reducing physical damage to the trichomes during growth and thus allowing for the retention of more secretions. However, this mechanism requires further verification.

Carbohydrate-derived pyruvic acid (PA) serves as a centralxnode linking glycolysis to lipid biosynthesis, acting as a precursor for acetyl-CoA, fatty acid chains, sterols, and terpenes via the tricarboxylic acid (TCA) cycle and methylerythritol phosphate (MEP) pathway (Li-Beisson et al., 2010). The synthesis of TP, secondary metabolites derived from PA metabolism, competes with lipid synthesis for common precursors like phosphoenolpyruvate (PEP). In ON30, PA levels declined sharply (−17.6%), suggesting enhanced carbon flux toward acetyl-CoA and downstream lipid synthesis (Table 2). This reduction is consistent with the hypothesis that carbon flux may be reallocated from the PA pool towards acetyl-CoA via the TCA cycle for lipid synthesis or towards other downstream products (Chen and Zhu, 2023), although this metabolic flux shift is inferred from static concentration changes. The enhancement of this metabolic pathway might indirectly lead to the downregulation of the phenylpropanoid pathway, thus explaining the decrease in TP content. PEPC is a key enzyme in carbon fixation and the supply of precursors for amino acid and organic acid synthesis, influencing overall metabolic flux and ultimately affecting the production of secondary metabolites (Baena et al., 2021). An increase in its activity not only accelerates the conversion of PA but also promotes the accumulation of TS and RS by enhancing carbon-nitrogen synergy. FAS plays a central role in the de novo synthesis of fatty acids, and its increased activity directly promotes fatty acid accumulation (Kitajima-Koga et al., 2020). Meanwhile, PAP is involved in lipid metabolism and the regulation of phospholipid levels, and the reactions it participates in are crucial steps in the synthesis and storage of triacylglycerols (Torres-Rodríguez et al., 2021). This study found that ON30 significantly enhanced the activities of FAS, PAP, and PEPC, which is likely a key factor contributing to the significant increase in the oil content score for this treatment. The significant increase in PEE content in ON30 and its strong correlation with oil content further corroborate this conclusion, as PEE comprises nearly all secondary metabolites in tobacco leaves (Jiang et al., 2025) and is often used as a proxy indicator for the “oil content” of flue-cured tobacco. However, a notable decline in enzyme activities (FAS, PEPC) and lipid metabolites (PEE, CBT-diol) was observed at the highest organic nitrogen level (ON40). This “diminishing return” or inhibitory effect likely stems from the carbon-nitrogen imbalance hypothesis. Excessive organic nitrogen input in ON40 may lead to an overabundance of nitrogen relative to carbon skeletons. According to the C/N balance theory, plants prioritize carbon resources for amino acid and protein synthesis to assimilate excess nitrogen, thereby diverting acetyl-CoA away from the lipid biosynthesis pathway (Fritz et al., 2006). Additionally, the significantly lower soil pH observed under high organic loading (ON40) might create a sub-optimal rhizosphere environment, suppressing the activity of soil microorganisms and root enzymes essential for nutrient uptake and metabolism (Feng et al., 2023). Meanwhile, it is crucial to acknowledge a potential spatial mismatch: while our enzyme activity assays represent whole-leaf protein levels, lipid biosynthesis relevant to “oiliness” is highly concentrated in the glandular trichomes. Future studies utilizing trichome-specific metabolite profiling are needed to directly link enzyme expression with localized lipid accumulation.

The sensory evaluation of oil content depends not only on chemical composition but is also closely related to the physical properties of the tobacco leaf (Cai et al., 2023). In our study, ON30 significantly reduced the softness value (Table 3), indicating that the leaves from this treatment were softer, which is consistent with the typical tactile characteristics of high-oil-content leaves (Li et al., 2008). The increase in tensile strength reflects an optimization of the leaf’s fibrous structure, which may be attributed to the regulation of cell wall components (such as cellulose and pectin) by the combined application of organic nitrogen (Ogden et al., 2018). Furthermore, the significant increase in leaf thickness for ON30, while remaining within the standard range (120-140 μm), ensured both good combustion performance and avoided difficulties in rolling due to excessive thickness. These physical traits contribute directly to the perception of oiliness, highlighting the multifactorial nature of oil content, where chemical and structural factors act synergistically. The deterioration of these physical traits (increased hardness, decreased thickness) in ON40 further supports the notion of a physiological stress response to nutrient over-supply, where rapid vegetative growth might occur at the expense of structural refinement and secondary metabolite accumulation. Importantly, our RF model identified Softness (Sof) as a top-three predictor. While softness itself is a physical state rather than a metabolic driver like FAS, it serves as a critical “accompanying trait” that integrates the outcomes of cell wall structure and moisture retention, making it a reliable indicator for sensory quality.

Despite the inherent challenges of small sample sizes in agricultural field trials, the Random Forest (RF) model, augmented with a rigorous data generation strategy (noise injection and linear interpolation), demonstrated robust explanatory power. The stark contrast between the poor performance on original data (CV R² = -0.039) and the high accuracy on augmented data (CV R² = 0.819) highlights the necessity of data augmentation in capturing latent biological patterns that are otherwise obscured by limited sampling (Shorten and Khoshgoftaar, 2019). By integrating RFE for feature selection, we distilled the model down to five core indicators—PEE, CBT-diol, Sof, PA, and GT—which collectively explain the majority of the variance in oil content.

The feature importance analysis (Figure 6) and PDP (Figure 7) quantified the relative contributions of these markers. PEE and CBT-diol emerged as the most dominant positive drivers, confirming that the abundance of lipid-related metabolites is the material basis of oiliness. Interestingly, PA exhibited a threshold-dependent negative effect, consistent with its role as a transient metabolic precursor that is depleted during active lipid synthesis (Zeng et al., 2025).

Crucially, the SHAP interaction analysis (Figure 8) moved beyond single-factor effects to reveal context-dependent regulation. We observed that the positive contribution of PEE to oil content was synergistic with CBT-diol levels. This reflects the physiological reality where PEE (a solvent mixture) and CBT-diol (a specific solute) co-vary to determine the sensory “oiliness.” Similarly, the negative impact of Softness appeared to be modulated by PEE levels, suggesting that high lipid content might partially compensate for less-than-ideal leaf texture. These non-linear and interactive insights provide a more nuanced understanding of quality formation than traditional linear correlation analyses.

While this study offers valuable insights into the regulation of tobacco oil content, several limitations must be acknowledged. First, the field experiment was conducted at a single location over one growing season. Given the known influence of environmental factors (temperature, rainfall, soil type) on secondary metabolism, multi-year and multi-location trials are necessary to validate the generalizability of the optimal 30% organic nitrogen ratio. Second, although the data augmentation strategy successfully mitigated overfitting and improved model convergence (as evidenced by the OOB error curve), the original sample size (N = 15 biological replicates) remains small for machine learning applications. The findings should be interpreted as exploratory hypothesis generation rather than definitive predictive rules. Future research should prioritize larger-scale datasets and incorporate multi-omics approaches (transcriptomics and metabolomics) to construct a comprehensive regulatory network spanning genes, metabolites, and sensory traits.