Metabolomic and transcriptomic analyses reveal differences in fatty acids in tobacco leaves across cultivars and developmental stages
Tobacco is a key economic crop with a long cultivation history, grown widely around the world. As a model species within the Solanaceae family, it is of significant scientific value for studies on gene function, molecular breeding, and other areas of research . Tobacco products undergo various curing processes, including flue curing, sun curing, and air curing, each of which significantly affects the chemical composition and sensory properties of the tobacco. The levels of starch, pigments, polyphenols, fatty acids, and other compounds in tobacco vary markedly under different curing conditions, with the curing process also driving the transformation of several substances . Fatty acids are an essential class of compounds in plants, playing a crucial role in their growth and development . They are key components of cell membrane lipids, important energy sources, and precursors of signaling molecules. Additionally, fatty acids are involved in cell recognition, species specificity, and tissue immunity. They also contribute to indirect insect resistance in plants , by inducing defense responses and utilizing pre-formed physical barriers, such as the cuticle and cell wall, to combat pathogens and other threats . Fatty acids are classified based on chain length into short-chain (C1-C6), medium-chain (C6-C12), and long-chain (> C12) fatty acids. Short- and medium-chain fatty acids and their derivatives exhibit extensive structural diversity and are incorporated into various biomolecules, serving as components of antibiotics, insect pheromones, and plant storage lipids . Long-chain fatty acids play a crucial role in maintaining normal plant development and preventing organ fusion .In tobacco, the composition, content, and distribution of fatty acids play a key regulatory role in its adaptation to various stress conditions . In addition to being essential for plant growth, development, and cell structure, fatty acids also significantly influence the quality and flavor of tobacco leaves . As tobacco leaves mature, both the content and composition of fatty acids change. Different tobacco varieties exhibit distinct patterns of organic acid variation, as confirmed by previous studies . The study found that spice tobacco contains higher levels of 3-methylvaleric acid, valeric acid, and isovaleric acid compared to flue-cured tobacco, while formic acid and acetic acid are relatively lower. In contrast, all organic acid levels in burley tobacco and Maryland tobacco are lower than those in spice and flue-cured tobaccos. Current research on tobacco fatty acids primarily focuses on leaf fatty acid composition , metabolic changes in primary, secondary, and lipid metabolism before and after topping , and the effects of curing methods (such as air-curing, sun-curing, and flue-curing) on the chemical composition of tobacco leaves . However, there is a lack of comprehensive studies analyzing fatty acid content variations across different growth stages and curing methods. We selected four tobacco varieties: K326, Basma, Samsun, and Cuba1. K326 is a flue-cured tobacco ( Nicotiana tabacum L.) variety known for its superior quality . Basma, an aromatic type, and Samsun, a taste type, are representative varieties of aromatic tobacco, originating from the Mediterranean region in the 1560 s . Cuba1 is one of the primary varieties used in cigar production. This study provides a comprehensive analysis of the fatty acid composition during the growth, development, and processing of four tobacco varieties. It investigates the differential expression of fatty acid related genes at specific stages, offering molecular insights into the accumulation patterns of fatty acids in tobacco. Moreover, the findings provide new perspectives for selecting roasting methods tailored to different tobacco varieties. 

Plant materials Four varieties samples from the seedling and transplant stages in the seedbed, as well as two stages during field growth, and applied two curing methods: flue-curing and air-curing. Using K326 as the reference variety, we compared the four tobacco varieties based on transcriptomics, targeted fatty acid metabolomics, and targeted short-chain fatty acid metabolomics. For all omics analyses, three independent biological replicates were used, with each replicate consisting of 3–6 tobacco plants of similar growth vigor. Samples were immediately flash-frozen in liquid nitrogen and stored at -80℃ until further use. Detailed information on the time, stage, usage, and omics category of material selection is shown in Table . Schematic diagram of the materials is shown in Fig. . The seeds of K326, Samsun, Basma, and Cuba1 are stored at the Technology Center of Hunan Tobacco Industry Co., Ltd. The seedbed stage cultivation conditions were as follows: soil cultivation in an artificial growth chamber with a temperature of 26–28 °C, humidity around 75%, light intensity of 4.3 K Lux, a 16-h photoperiod with 8 h of darkness, and a carbon dioxide concentration of approximately 672 PPM. The field growth stage was conducted in the experimental field at Liuyang City, Hunan Province. Curing status: Flue-cured tobacco was processed using a three-stage curing method, while air-cured tobacco was cured until the main veins were dry. The curing procedures for both flue-cured and air-cured tobacco followed the specific method outlined by Jie Chen et al. . Construction of the differential fatty acid metabolic profile for different tissue types of the K326 variety, samples were collected from mature plants at the squaring stage (55 days after transplanting). A total of 15 tissue samples were analyzed (Fig. a), including: flower buds (FB), upper stem (US), upper leaf 1 (UL), upper leaf 2 (U2L), upper leaf 1 vein (ULV), upper leaf 2 vein (U2LV), middle stem (MS), middle leaf (ML), middle leaf vein (MLV), lower leaf 1 (LL), lower leaf 2 (L2L), lower leaf 1 vein (LLV), lower leaf 2 vein (L2LV), lower stem (LS), and root (R). Based on the average total number of leaves at the squaring stage, the stems were divided into three sections according to plant height, and the leaves were classified into five parts according to their position on the plant. Tissue samples were collected from three plants of uniform growth, with three biological replicates. 

Four varieties samples from the seedling and transplant stages in the seedbed, as well as two stages during field growth, and applied two curing methods: flue-curing and air-curing. Using K326 as the reference variety, we compared the four tobacco varieties based on transcriptomics, targeted fatty acid metabolomics, and targeted short-chain fatty acid metabolomics. For all omics analyses, three independent biological replicates were used, with each replicate consisting of 3–6 tobacco plants of similar growth vigor. Samples were immediately flash-frozen in liquid nitrogen and stored at -80℃ until further use. Detailed information on the time, stage, usage, and omics category of material selection is shown in Table . Schematic diagram of the materials is shown in Fig. . The seeds of K326, Samsun, Basma, and Cuba1 are stored at the Technology Center of Hunan Tobacco Industry Co., Ltd. The seedbed stage cultivation conditions were as follows: soil cultivation in an artificial growth chamber with a temperature of 26–28 °C, humidity around 75%, light intensity of 4.3 K Lux, a 16-h photoperiod with 8 h of darkness, and a carbon dioxide concentration of approximately 672 PPM. The field growth stage was conducted in the experimental field at Liuyang City, Hunan Province. Curing status: Flue-cured tobacco was processed using a three-stage curing method, while air-cured tobacco was cured until the main veins were dry. The curing procedures for both flue-cured and air-cured tobacco followed the specific method outlined by Jie Chen et al. . Construction of the differential fatty acid metabolic profile for different tissue types of the K326 variety, samples were collected from mature plants at the squaring stage (55 days after transplanting). A total of 15 tissue samples were analyzed (Fig. a), including: flower buds (FB), upper stem (US), upper leaf 1 (UL), upper leaf 2 (U2L), upper leaf 1 vein (ULV), upper leaf 2 vein (U2LV), middle stem (MS), middle leaf (ML), middle leaf vein (MLV), lower leaf 1 (LL), lower leaf 2 (L2L), lower leaf 1 vein (LLV), lower leaf 2 vein (L2LV), lower stem (LS), and root (R). Based on the average total number of leaves at the squaring stage, the stems were divided into three sections according to plant height, and the leaves were classified into five parts according to their position on the plant. Tissue samples were collected from three plants of uniform growth, with three biological replicates. 

Targeted short-chain fatty acid metabolomics The method for short-chain fatty acid determination is based on references . After sample extraction, short-chain fatty acids (SCFAs) were analyzed by LC–MS/MS in multiple reaction monitoring (MRM) mode. The stock solution of individual SCFAs was mixed and prepared in an SCFA-free matrix to generate a series of SCFA calibrators. Isotope-labeled standards were prepared and mixed as the internal standard (IS). For sample preparation, 0.02 g of the ground sample was weighed and suspended in 80 μL of 80% methanol–water solution, followed by mixing. The sample was then centrifuged at 12,000 rpm for 10 min at 4 °C. The supernatant (50 μL) was transferred and mixed with 150 μL of derivatization reagent. The derivatized samples were incubated at 40 °C for 40 min. After derivatization, the samples were diluted with 80% methanol–water solution at 1 × and 100 × concentrations. A total of 95 μL of the supernatant was mixed with 5 μL of a mixed internal standard solution in 80% methanol–water, and the mixture was thoroughly mixed before LC–MS analysis. Short-chain fatty acids (SCFAs) were quantified using an ultra-high performance liquid chromatography coupled to tandem mass spectrometry (UHPLC-MS/MS) system (Vanquish™ Flex UHPLC-TSQ Altis™, Thermo Scientific Corp., Germany), operated by Novogene Co., Ltd. (Beijing, China). Separation was performed on a Waters ACQUITY UPLC BEH C18 column (2.1 × 100 mm, 1.7 μm), maintained at 40 °C. The mobile phase consisted of 10 mM ammonium acetate in water (solvent A) and a 1:1 mixture of acetonitrile and isopropanol (solvent B), delivered at a flow rate of 0.30 mL/min. The solvent gradient was programmed as follows: initial 25% B for 2.5 min; 25–30% B for 3 min; 30–35% B for 3.5 min; 35–38% B for 4 min; 38–40% B for 4.5 min; 40–45% B for 5 min; 45–50% B for 5.5 min; 50–55% B for 6.5 min; 55–58% B for 7 min; 58–70% B for 7.5 min; 70–100% B for 7.8 min; 100–25% B for 10.1 min; and 25% B for 12 min. The mass spectrometer was operated in negative multiple reaction monitoring (MRM) mode with the following parameters: IonSpray Voltage (-4500 V), Sheath Gas (35 psi), Ion Source Temperature (550℃), Auxiliary Gas (50 psi), and Collision Gas (55 psi). The selection of differential short-chain fatty acids (SCFAs) was based on fold change (FC) and P -value, with the following thresholds: FC > 1.2 or FC < 0.833, and P -value < 0.05 . Targeted fatty acid metabolomics analysis The method for fatty acid determination is referenced from . The stock solution of individual fatty acids was prepared by mixing them in a fatty acid-free matrix to generate a series of fatty acid calibrators. Isotope-labeled standards were prepared and mixed as the internal standard (IS). Samples (100 mg) were frozen in liquid nitrogen, homogenized with 300 μL of isopropanol/acetonitrile (1:1) solution, and subjected to ultrasound treatment for 10 min. After 60 min at -20 °C, the samples were centrifuged at 12,000 rpm for 10 min. The supernatant (50 μL) was mixed with 150 μL of derivatization reagent and incubated at 40 °C for 40 min. Following derivatization, 47.5 μL of the supernatant was combined with 2.5 μL of a mixed internal standard solution. The final preparation was injected into the LC–MS/MS system for analysis. Fatty acids were quantified using an ultra-high performance liquid chromatography coupled to tandem mass spectrometry (UHPLC-MS/MS) system (ExionLC™ AD UHPLC-QTRAP 6500 + , AB SCIEX Corp., Boston, MA, USA) at Novogene Co., Ltd. (Beijing, China). Separation was carried out on a Waters ACQUITY UPLC BEH C18 column (2.1 × 100 mm, 1.7 μm), maintained at 40 °C. The mobile phase consisted of 0.1% formic acid in acetonitrile/water (1:1) (solvent A) and isopropanol/acetonitrile (1:1) (solvent B), delivered at a flow rate of 0.30 mL/min. The solvent gradient was programmed as follows: initial 45% B for 1 min; 45–70% B for 4.5 min; 70–75% B for 9 min; 75–80% B for 12.5 min; 80–100% B for 14 min; 100–45% B for 15.1 min; and 45% B for 17 min.The mass spectrometer was operated in negative multiple reaction monitoring (MRM) mode with the following parameters: IonSpray Voltage (-4500 V), Curtain Gas (35 psi), Ion Source Temperature (550℃), Ion Source Gas 1 and 2 (60 psi). Total RNA extraction and cDNA synthesis Four varieties of Middle leaf blade in the Squaring stage were used as materials (Table , 55 days after transplanting, -SS). Total RNA was extracted using the RNAprep Pure Plant Plus Kit (Tiangen, China), and cDNA synthesis was carried out using the All-in-One First-Strand cDNA Synthesis SuperMix Kit (NovoScript Plus, China). Detailed protocols are provided in the kit manuals and are based on previously published methods . The RNA concentration was measured using a NanoDrop 2000 and the RNA integrity was assayed by agarose gel electrophoresis. Only RNA samples with an OD260/280 ratio between 1.8 and 2.2 and an OD260/230 ratio greater than 2.0 that showed three discrete bands of 28S, 18S and 5S were used for cDNA synthesis. The reverse transcription system of each sample was as follows: RNA template (1 μg), gDNA Purge (1 μL), Supermix (10 μL) and RNase Free Water (up to 20 μL) at 50 °C for 30 min and 75 °C for 5 min. Quantitative RT–PCR validation of candidate genes Quantitative RT–PCR (qRT–PCR) was performed on a Roche LightCycler 480. Gene-specific primers were designed using Primer Premier 5.0 software and the NCBI primer BLAST tool (Table ). Actin7 was used as the reference gene (LOC107795948) . Detailed procedures are outlined in the LightCycler® Multiplex Master qPCR Reaction Mix kit manual. Relative gene expression levels were calculated using the 2 −ΔΔCt method. Detailed protocols are provided in the kit manuals and are based on previously published methods . Transcriptome sequencing and analysis RNA sequencing was conducted by Novogene (Beijing, China, https://en.novogene.com/ ). The reference genome and gene model annotation files were obtained directly from the genome website ( http://lifenglab.hzau.edu.cn/Nicomics/ ) . Gene function annotation was performed using the following databases: KEGG ( https://www.kegg.jp/kegg/ko.html ) and GO ( http://geneontology.org/ ). Differential expression analysis was carried out using DESeq2 software (version 1.20.0) with the following criteria: padj ≤ 0.05 and |log2FoldChange|≥ 1.0 for each comparison . For detailed methods, analysis procedures, and software, showed in Supplementary Material 5. The original data of differentially expressed genes and pathway enrichment information are shown in Table . Data analysis Data analysis was performed using GraphPad Prism 8 and Excel 2016, results are presented as the mean ± SD. One-way analysis of variance (ANOVA) was conducted, followed by Tukey’s post-hoc test, with statistical significance set at P < 0.05. All data were derived from three biological replicates. For the metabolomics differential metabolite clustering heatmap, the data were standardized using the following formula: (content value—metabolite mean) / standard deviation. For the transcriptomics differential gene clustering heatmap, FPKM-normalized values were used (Log2(FPKM + 1)) . 

The method for short-chain fatty acid determination is based on references . After sample extraction, short-chain fatty acids (SCFAs) were analyzed by LC–MS/MS in multiple reaction monitoring (MRM) mode. The stock solution of individual SCFAs was mixed and prepared in an SCFA-free matrix to generate a series of SCFA calibrators. Isotope-labeled standards were prepared and mixed as the internal standard (IS). For sample preparation, 0.02 g of the ground sample was weighed and suspended in 80 μL of 80% methanol–water solution, followed by mixing. The sample was then centrifuged at 12,000 rpm for 10 min at 4 °C. The supernatant (50 μL) was transferred and mixed with 150 μL of derivatization reagent. The derivatized samples were incubated at 40 °C for 40 min. After derivatization, the samples were diluted with 80% methanol–water solution at 1 × and 100 × concentrations. A total of 95 μL of the supernatant was mixed with 5 μL of a mixed internal standard solution in 80% methanol–water, and the mixture was thoroughly mixed before LC–MS analysis. Short-chain fatty acids (SCFAs) were quantified using an ultra-high performance liquid chromatography coupled to tandem mass spectrometry (UHPLC-MS/MS) system (Vanquish™ Flex UHPLC-TSQ Altis™, Thermo Scientific Corp., Germany), operated by Novogene Co., Ltd. (Beijing, China). Separation was performed on a Waters ACQUITY UPLC BEH C18 column (2.1 × 100 mm, 1.7 μm), maintained at 40 °C. The mobile phase consisted of 10 mM ammonium acetate in water (solvent A) and a 1:1 mixture of acetonitrile and isopropanol (solvent B), delivered at a flow rate of 0.30 mL/min. The solvent gradient was programmed as follows: initial 25% B for 2.5 min; 25–30% B for 3 min; 30–35% B for 3.5 min; 35–38% B for 4 min; 38–40% B for 4.5 min; 40–45% B for 5 min; 45–50% B for 5.5 min; 50–55% B for 6.5 min; 55–58% B for 7 min; 58–70% B for 7.5 min; 70–100% B for 7.8 min; 100–25% B for 10.1 min; and 25% B for 12 min. The mass spectrometer was operated in negative multiple reaction monitoring (MRM) mode with the following parameters: IonSpray Voltage (-4500 V), Sheath Gas (35 psi), Ion Source Temperature (550℃), Auxiliary Gas (50 psi), and Collision Gas (55 psi). The selection of differential short-chain fatty acids (SCFAs) was based on fold change (FC) and P -value, with the following thresholds: FC > 1.2 or FC < 0.833, and P -value < 0.05 . 

The method for fatty acid determination is referenced from . The stock solution of individual fatty acids was prepared by mixing them in a fatty acid-free matrix to generate a series of fatty acid calibrators. Isotope-labeled standards were prepared and mixed as the internal standard (IS). Samples (100 mg) were frozen in liquid nitrogen, homogenized with 300 μL of isopropanol/acetonitrile (1:1) solution, and subjected to ultrasound treatment for 10 min. After 60 min at -20 °C, the samples were centrifuged at 12,000 rpm for 10 min. The supernatant (50 μL) was mixed with 150 μL of derivatization reagent and incubated at 40 °C for 40 min. Following derivatization, 47.5 μL of the supernatant was combined with 2.5 μL of a mixed internal standard solution. The final preparation was injected into the LC–MS/MS system for analysis. Fatty acids were quantified using an ultra-high performance liquid chromatography coupled to tandem mass spectrometry (UHPLC-MS/MS) system (ExionLC™ AD UHPLC-QTRAP 6500 + , AB SCIEX Corp., Boston, MA, USA) at Novogene Co., Ltd. (Beijing, China). Separation was carried out on a Waters ACQUITY UPLC BEH C18 column (2.1 × 100 mm, 1.7 μm), maintained at 40 °C. The mobile phase consisted of 0.1% formic acid in acetonitrile/water (1:1) (solvent A) and isopropanol/acetonitrile (1:1) (solvent B), delivered at a flow rate of 0.30 mL/min. The solvent gradient was programmed as follows: initial 45% B for 1 min; 45–70% B for 4.5 min; 70–75% B for 9 min; 75–80% B for 12.5 min; 80–100% B for 14 min; 100–45% B for 15.1 min; and 45% B for 17 min.The mass spectrometer was operated in negative multiple reaction monitoring (MRM) mode with the following parameters: IonSpray Voltage (-4500 V), Curtain Gas (35 psi), Ion Source Temperature (550℃), Ion Source Gas 1 and 2 (60 psi). 

Four varieties of Middle leaf blade in the Squaring stage were used as materials (Table , 55 days after transplanting, -SS). Total RNA was extracted using the RNAprep Pure Plant Plus Kit (Tiangen, China), and cDNA synthesis was carried out using the All-in-One First-Strand cDNA Synthesis SuperMix Kit (NovoScript Plus, China). Detailed protocols are provided in the kit manuals and are based on previously published methods . The RNA concentration was measured using a NanoDrop 2000 and the RNA integrity was assayed by agarose gel electrophoresis. Only RNA samples with an OD260/280 ratio between 1.8 and 2.2 and an OD260/230 ratio greater than 2.0 that showed three discrete bands of 28S, 18S and 5S were used for cDNA synthesis. The reverse transcription system of each sample was as follows: RNA template (1 μg), gDNA Purge (1 μL), Supermix (10 μL) and RNase Free Water (up to 20 μL) at 50 °C for 30 min and 75 °C for 5 min. 

Quantitative RT–PCR (qRT–PCR) was performed on a Roche LightCycler 480. Gene-specific primers were designed using Primer Premier 5.0 software and the NCBI primer BLAST tool (Table ). Actin7 was used as the reference gene (LOC107795948) . Detailed procedures are outlined in the LightCycler® Multiplex Master qPCR Reaction Mix kit manual. Relative gene expression levels were calculated using the 2 −ΔΔCt method. Detailed protocols are provided in the kit manuals and are based on previously published methods . 

RNA sequencing was conducted by Novogene (Beijing, China, https://en.novogene.com/ ). The reference genome and gene model annotation files were obtained directly from the genome website ( http://lifenglab.hzau.edu.cn/Nicomics/ ) . Gene function annotation was performed using the following databases: KEGG ( https://www.kegg.jp/kegg/ko.html ) and GO ( http://geneontology.org/ ). Differential expression analysis was carried out using DESeq2 software (version 1.20.0) with the following criteria: padj ≤ 0.05 and |log2FoldChange|≥ 1.0 for each comparison . For detailed methods, analysis procedures, and software, showed in Supplementary Material 5. The original data of differentially expressed genes and pathway enrichment information are shown in Table . 

Data analysis was performed using GraphPad Prism 8 and Excel 2016, results are presented as the mean ± SD. One-way analysis of variance (ANOVA) was conducted, followed by Tukey’s post-hoc test, with statistical significance set at P < 0.05. All data were derived from three biological replicates. For the metabolomics differential metabolite clustering heatmap, the data were standardized using the following formula: (content value—metabolite mean) / standard deviation. For the transcriptomics differential gene clustering heatmap, FPKM-normalized values were used (Log2(FPKM + 1)) . 

Fatty acid profiling in K326 at the squaring stage In this study, targeted fatty acid metabolomics analysis was conducted on K326 plants at the flowering stage (55 days after transplanting), with 15 tissue samples analyzed. A total of 11 short-chain fatty acids (SCFAs) were measured. Among the 15 tissue samples, the contents of SCFAs in flower buds (FB) and upper leaves (UL) were higher than in other tissues, except for 3-Methylvalerate and 2-Methylvalerate. Cluster analysis (Fig. b) revealed that the samples were grouped into five clusters, with FB and UL forming one group, exhibiting significantly higher relative content compared to other tissues. Flower buds and upper leaves (UL) clustered together, and the 11 measured SCFAs were grouped into four categories: 2-Methylvalerate and 3-Methylvalerate, which showed lower levels in flower and upper tissues compared to other parts; 4-Methylvaleric acid, 2-Methylbutyrate, and Isobutyric acid, which were higher in flower tissues compared to other parts. The total content of short-chain fatty acids in the leaves followed the pattern: UL > U2L > LL > ML > L2L (Fig. c), and in the stems: US > MS > LS. The total short-chain fatty acid content in the upper leaves, stems, and flower buds was higher than in the middle and lower parts of the plant. This suggests that the short-chain fatty acid levels in K326 at the flowering stage were generally higher in the upper parts of the plant compared to the middle and lower sections, although the root did not show the lowest content. In the targeted fatty acid metabolomics analysis, 50 fatty acids were measured, including all but gamma-Linolenic acid, across the various tissue samples (C8:0-C22:6). The fatty acid clustering heatmap (Fig. d) indicated that the tissue samples could be grouped into three categories: roots as one group, flower buds, all leaves, and the main veins of upper leaves as another group, and stems and other main leaf veins as the third group. Significant differences in fatty acid content were observed between the upper, middle, and lower leaf, vein, and stem tissues, with distinct separation based on fatty acid content. The root had higher levels of nine fatty acids, including Octadecanoic acid, Hexadecanoic acid, and Heptadecanoic acid, compared to other tissues. Flower buds exhibited elevated levels of Linoelaidic acid, Linoleic acid, and cis-13,16-Docosadienoic acid, compared to other tissues. Upper leaves (UL) contained higher levels of Caprylic acid, cis-10-Heptadecenoic acid, and trans-10-Pentadecenoic acid compared to other tissues. The specific content data for short-chain fatty acids and fatty acid metabolomics measurements are provided in Table . Analysis of short-chain fatty acid content variations in leaves during the seedling stage The short-chain fatty acid content in the leaves at two stages of the seedling bed period was measured. Among the 11 short-chain fatty acids detected, Acetic acid had the highest content in all tobacco varieties, accounting for 78–95% of the total short-chain fatty acids, with a higher proportion in the (transplant stag) ML stage compared to the (Seedling stage) S stage (Table ). In terms of total short-chain fatty acid content, K326 had the highest content in the S stage, followed by Basma, Cuba1, and Samsun. In the ML stage, Basma showed the highest content, closely followed by K326, while Cuba1 and Samsun had lower levels. For 3-Methylvaleric acid, Cuba1 had the highest proportion (0.18%) at the S stage, compared to the other three varieties. At the ML stage, Basma and Samsun had significantly higher proportions (0.18% and 0.16%, respectively) compared to K326 and Cuba1 (0.04%). In the case of Isobutyric acid, Cuba1 exhibited a lower proportion at the S stage but had a significantly higher proportion at the ML stage than the other varieties. For the C5 isomer/normal ratio, Cuba1 had the highest value at the S stage, while the ratio in the ML stage was lower than that of the other varieties. The C6 isomer/normal ratio in K326 was notably lower in the ML stage compared to the other varieties. At the S stage, K326 showed higher levels of differential short-chain fatty acids compared to the other three varieties. At the ML stage, Basma and Cuba1 exhibited downregulation of 2-Methylbutyrate compared to K326. Specific fold changes (FC > 1.2 or FC < 0.833 and P -value < 0.05) for the differential short-chain fatty acids are provided in Table . As shown in Fig. a, principal component analysis (PCA) based on PC1 (54.03%) and PC2 (21.09%) clearly differentiated the short-chain fatty acid profiles of Cuba1 and Basma. In the MS stage, PCA (Fig. d) with PC1 (36.46%) and PC2 (21.08%) showed that Samsun and Basma had similar profiles, with reduced differences in the first principal component compared to the S stage. Figure b, e shows the specific content of 11 short-chain fatty acids at two stages. According to the heatmap of normalized short-chain fatty acid content (Fig. c, f), most short-chain fatty acids in K326 at the S stage were higher than in the other varieties. Comparison between stages revealed that Acetic acid, Isovaleric acid, and 2-Methylbutyrate were higher in the ML stage, while Propionic acid, Isobutyric acid, Butyric acid, and Valeric acid showed the opposite trend (Fig. g). Gene expression analysis of tobacco leaves at the seedling bed stage Differential gene expression was analyzed in tobacco leaves from four varieties at two stages of the seedling bed period: seedling stage (S) and mature leaf stage (ML). A total of 158.25 GB of clean data was obtained from 28 samples, with an average of 6.59 GB per sample. The average total mapping rate to the reference genome was 93.52%. At the S stage, compared to K326, Samsun showed 1973 differentially expressed genes (DEGs), with 906 upregulated and 1067 downregulated (DESeq2 padj ≤ 0.05, |log2FoldChange|≥ 1.0). Basma exhibited 3433 DEGs, with 1212 upregulated and 2221 downregulated. Cuba1 showed 2705 DEGs, with 961 upregulated and 1744 downregulated (Fig. a). At the ML stage, Samsun compared to K326 had 3383 DEGs, with 1969 upregulated and 1414 downregulated. Basma had 5020 DEGs, with 2161 upregulated and 2859 downregulated, while Cuba1 showed 2574 DEGs, with 1308 upregulated and 1266 downregulated (Fig. b). Basma showed the highest number of DEGs at both stages among the three varieties compared to K326. PCA indicated distinct separation of S and ML stages along PC1, with a more scattered distribution in the ML stage compared to S. Intra-variety comparisons in the ML stage showed greater variability (Fig. c). Venn diagram analysis of the common DEGs between K326 and the other three varieties revealed 710 DEGs at the S stage and 879 at the ML stage, with 457 shared DEGs between the two stages. This DEG set represents the major differences between K326 and the other varieties during the seedling bed period (Fig. d). GO enrichment analysis showed significant enrichment in “signal transduction” and “signaling” pathways (padj < 0.05, Fig. f). KEGG pathway analysis indicated significant enrichment in pathways such as “Various types of N-glycan biosynthesis,” “N-Glycan biosynthesis,” and “RNA degradation” ( p -value < 0.05, padj = 0.52, Fig. e). In the fatty acid-related biosynthesis pathways, genes involved in “Fatty acid degradation,” “Fatty acid biosynthesis,” and “Fatty acid metabolism” ( p -value < 0.24, padj = 0.52–0.61) were enriched. Three genes ( Nta01g31980 , Nta08g22780 , and Nta23g11140 ) were identified in these pathways, annotated as acyl-activating enzyme, alcohol dehydrogenase, and enoyl-[acyl-carrier-protein] reductase. These genes were not expressed in K326 at either stage, but were expressed in the other three varieties. Based on gene expression data, a fatty acid pathway heatmap for the four varieties during the seedling bed period was generated, showing pathways such as “Fatty acid biosynthesis” (nta00061), “Fatty acid elongation” (nta00062), and “Fatty acid degradation” (nta00071) (Fig. ). In the “Fatty acid biosynthesis” pathway (nta00061), the relative expression levels of 13 ACC genes, 2 MT genes, and 4 enoyl-[acyl-carrier-protein] reductase I genes were significantly higher in the S stage compared to the ML stage. Four KAS III genes and six TE genes showed half of the genes with higher relative expression in either the S or ML stage. Varieties showed varying expression patterns across different genes. For example, Nta23g00340 (KAS II) was not expressed in the ML stage, and Nta12g11030 (acyl-[acyl-carrier-protein] desaturase) was not expressed in the S stage. In the “Fatty acid elongation” pathway (nta00062), two palmitoyl-protein thioesterases and two palmitoyl-CoA hydrolases had higher relative expression in the ML stage compared to the S stage, while two long-chain 3-oxoacyl-CoA reductases showed higher relative expression in the S stage. K326 showed higher relative expression than the other varieties for these genes. Nta06g11470 and Nta05g12800 were not expressed in the ML stage. In the “Fatty acid degradation” pathway (nta00071), the Nta12g16270 and Nta11g17320 were not expressed in the ML stage for all varieties. Two acetyl-CoA C-acyltransferase genes showed significantly higher expression in K326 at the ML stage compared to the other varieties and stages. Two long-chain fatty acid omega-monooxygenase genes had higher expression in K326 during the S stage compared to the other varieties and stages. Analysis of the difference in fatty acid content in leaves before and after topping The short-chain fatty acid (SCFA) content in the middle leaves of four varieties was measured at the budding stage (55 days post-transplantation) and 10 days after topping (70 days post-transplantation). Among the 11 measured SCFAs, the total content at the budding stage was ranked as Basma > Samsun > Cuba1 > K326, while the ranking reversed after topping. The proportion of 3-methylvaleric acid was significantly higher in Basma and Samsun than in Cuba1 and K326 both before and after topping. The trend of isovaleric acid proportion was consistent at the budding stage, but after topping, all four varieties showed similar trends. The C5 isomer/normal ratio was higher in Basma and Samsun than in Cuba1 and K326 at the budding stage, while after topping, K326 and Basma had significantly higher ratios than Samsun and Cuba1 (Table ). Differential metabolite analysis (FC > 1.2 or FC < 0.833, P -value < 0.05) and within-group significance analysis revealed that at the budding stage, the primary differentially expressed SCFAs were valeric acid and hexanoic acid, which showed similar trends across varieties. After topping, Basma and Cuba1 had significantly higher levels of 3-methylvalerate than K326. Detailed information on the differentially expressed SCFAs across varieties is provided in Table . As shown in Fig. a, the PCA plot based on PC1 (83.7%) and PC2 (7.15%) clearly clustered the samples into two groups: the K326 and Cuba1 groups were closer to each other, while the Samsun and Basma groups were closer. After topping, K326 and Cuba1 were more distinctly separated, while Basma and Samsun clustered together (Fig. d). Figure c, f shows the specific content of 11 short-chain fatty acids at two stages. The contents of acetic acid, isobutyric acid, butyric acid, valeric acid, isovaleric acid, 2-methylbutyrate, hexanoic acid, and 2-methylvalerate were significantly lower after topping (Fig. g). After topping, K326 showed higher levels of all fatty acids compared to other varieties, except for 2-methylvalerate, propionic acid, and 3-methylvalerate. Samsun had significantly higher levels of 2-methylvalerate, propionic acid, and 3-methylvalerate compared to the other varieties. Cluster heatmap analysis (Fig. b, e) revealed that after topping, the short-chain fatty acid content patterns in Basma and Samsun became more similar compared to those before topping. Additionally, the content of most short-chain fatty acids in K326 after topping was higher than in the other varieties. The PCA plot for the budding stage (Fig. h), based on PC1 (54.43%) and PC2 (15.68%), showed that K326 was clearly separated from the other three varieties along PC1. After topping, Basma and Cuba1 clustered together, while K326 and Samsun grouped similarly. The two groups were distinctly separated along PC2. Among the 49 measured fatty acids, Basma, Samsun, and Cuba1 had 13, 16, and 5 differentially expressed fatty acids compared to K326 at the budding stage, respectively, and 6, 13, and 22 differentially expressed fatty acids after topping (Fig. i). At the budding stage, Samsun showed a downregulation of fatty acids compared to K326, with pentadecanoic acid showing the highest downregulation (Log2FC = -2.09). Basma also showed a downregulation of fatty acids compared to K326, except for linoelaidic acid, with trans-10-heptadecenoic acid being the most downregulated (Log2FC = -2.07). Cuba1 showed the highest downregulation of cis-11,14,17-eicosatrienoic acid (Log2FC = -2.02) compared to K326. Common differentially expressed fatty acids between Basma and Samsun compared to K326 included pentadecanoic acid, palmitoleic acid, and heptadecanoic acid. After topping, Samsun showed downregulation of 5 fatty acids compared to K326, with caprylic acid being the most downregulated (Log2FC = -0.62). Basma downregulated 12 fatty acids compared to K326, with cis-5,8,11,14,17-eicosapentaenoic acid showing the highest downregulation (Log2FC = -1.64). Cuba1 downregulated 21 fatty acids, with pentadecanoic acid showing the highest downregulation (Log2FC = -1.30). The common differentially expressed fatty acids between Basma, Samsun, and Cuba1 compared to K326 included decanoic acid and hendecanoic acid. Based on the clustering heatmap of the 49 measured fatty acids (Fig. j), the fatty acid content differences were classified into three groups: Basma and Cuba1 after topping, Basma and Samsun at the budding stage, and the remaining four groups. Gene expression analysis of leaves in the Squaring stage Differential gene analysis was performed on the leaves of four tobacco varieties at the squaring stage. A total of 78.04 G of clean bases were obtained from 12 samples, with an average of 6.50 G per sample. The average total alignment to the reference genome was 94.45%. Principal component analysis (PCA) revealed that Basma and Samsun were closely related, while K326 and Cuba1 were clearly distinct (Fig. a). In the differential gene analysis, 18,600 DEGs were identified between Samsun and K326, with 9,067 upregulated and 9,533 downregulated (DESeq2 padj ≤ 0.05 |log2FoldChange|≥ 1.0). For Basma vs. K326, 17,770 DEGs were detected, including 8,520 upregulated and 9,250 downregulated. For Cuba1 vs. K326, 15,537 DEGs were found, with 6,678 upregulated and 8,859 downregulated (Fig. b). Using K326 as the control, differential gene analysis was conducted for the other three varieties, and a Venn diagram was used to identify a common set of 9,456 DEGs across the two stages, which represents the major differential gene set for K326 vs. the other varieties at the squaring stage (Fig. c). Enrichment analysis of this gene set showed that 37 GO terms were significantly enriched (padj < 0.05). Notably, the terms “sequence-specific DNA binding,” “unfolded protein binding,” and “calcium ion binding” were most enriched (Fig. d). KEGG analysis identified 10 pathways with significant enrichment (padj < 0.05), with the most significant being “MAPK signaling pathway – plant,” “Photosynthesis—antenna proteins,” and “Endocytosis” (Fig. e). In the fatty acid-related synthesis and metabolism pathways, 35 genes, including KAT1 , KCS1 , AIM1 , and FAD2 were enriched in pathways such as “Fatty acid degradation,” “Fatty acid biosynthesis,” and “Fatty acid metabolism.” Clustering analysis of these 35 genes revealed similar expression patterns in Basma and Samsun, distinct from K326 and Cuba1. Notably, genes such as Nta02g06150 ( KAT1 ), Nta12g29320 ( KAT1 ), Nta12g07350 ( KCS11 ), and Nta11g08380 ( KCS11 ) showed significantly higher expression in K326 compared to the other varieties (Fig. f). Six differentially DEGs were selected for qRT-PCR validation of their expression levels, including ACOT8 ( Nta02g19100 ), AL7B4 ( Nta05g24430 ), KCS11 ( Nta16g12070 ), KCS20 ( Nta20g17200 ), KCS4 ( Nta08g20960 ) and LACS8 ( Nta17g22980 ). These genes are DEGs from four fatty acid-related metabolic pathways (Fig. f).The expression data obtained by RNA-Seq and qRT–PCR were compared. The genes related to the fatty acid pathway ACOT8 , AL7B4 , KCS4 , LACS8 had exactly the same trends as the transcriptome results (Fig. ). Analysis of fatty acid content in leaves subjected to different curing methods Short-chain fatty acids and fatty acids were analyzed in the middle leaves of four tobacco varieties following two different curing methods. The PCA plot based on 11 short-chain fatty acids and principal components PC1 (43.51%) and PC2 (30.95%) clearly differentiated between the curing methods (Fig. a). Except for the air-cured Cuba1 (Cuba1_AC), the second principal component effectively separated the air-cured and flue-cured samples, with air-cured samples clustering at the top and flue-cured samples at the bottom. The Cuba1_AC sample exhibited a more scattered distribution in the principal components. Figure b shows the number of differential short-chain fatty acids both between varieties and within varieties under different curing methods. Generally, air-cured (AC) samples of the same variety exhibited fewer downregulated or unchanged short-chain fatty acids compared to flue-cured (FC) samples, indicating higher short-chain fatty acid content in air-cured leaves. The differential short-chain fatty acid clustering heatmap (Fig. c) revealed that Samsun_AC and Basma_AC clustered together, as did Samsun_FC and K326_FC, while Basma_FC and Cuba1_FC formed another cluster. Overall, the short-chain fatty acids were grouped into two categories: Butyric acid, Hexanoic acid, Isobutyric acid, Valeric acid, and 2-Methylbutyrate clustered together, showing upregulation in Samsun_AC, Basma_AC, and K326_AC, and downregulation in the remaining samples. Figure d presents the trends in short-chain fatty acid content before and after curing. With the exception of Propionic acid, 2-Methylvalerate, and 3-Methylvalerate, all other short-chain fatty acids showed varying degrees of increase after curing. For K326, the content of short-chain fatty acids (except for Propionic acid and 2-Methylvalerate) followed the order: air-cured > flue-cured > pre-curing. Similar trends were observed for Isobutyric acid, Butyric acid, the three C5 acids, and Hexanoic acid in Basma and Samsun. Additionally, the changes in the three C5 acids before and after curing were consistent across all four varieties. A targeted fatty acid metabolomics analysis was performed on the processed tobacco leaves, identifying a total of 47 fatty acids. A PCA plot (Fig. a) based on PC1 (51.45%) and PC2 (17.51%) showed a significant distinction in fatty acid composition between air-cured and flue-cured tobacco leaves along the first principal component, with the flue-cured samples being more tightly clustered. Differential analysis of the fatty acids revealed that, under air-curing, Samsun, Basma, and Cuba1 had 21, 8, and 13 differential fatty acids compared to K326, respectively. Under flue-curing, Samsun, Basma, and Cuba1 showed 7, 13, and 5 differential fatty acids compared to K326. Notably, air-cured Samsun and Basma had more downregulated differential fatty acids compared to K326, while flue-cured samples showed the opposite trend (Fig. b). When comparing different curing methods for the same variety, K326 showed 2 upregulated and 15 downregulated differential fatty acids under flue-curing compared to air-curing. The most downregulated fatty acid was Tridecanoic acid (log2FC: -2.66). Four fatty acids, including Tridecanoic acid, cis-10-Heptadecenoic acid, Docosanoic acid, and Heneicosanoic acid, showed a decrease (log2FC < -1) across all varieties. In contrast, gamma-Linolenic acid and Arachidonic acid showed upregulation (log2FC > 0.5). For air-cured varieties, Samsun, Basma, and Cuba1 showed three shared downregulated differential fatty acids compared to K326, namely alpha-Linolenic acid, cis-11,14,17-Eicosatrienoic acid, and homo-gamma-Linolenic acid. Fatty acid clustering analysis (Fig. c) based on differences in fatty acid content showed that Samsun_AC and Basma_AC clustered together, while K326_AC and Cuba1_AC formed another group. The remaining samples each formed separate clusters. Overall, the measured fatty acids were divided into three categories: one category showed downregulation in flue-cured compared to air-cured, another showed upregulation, and the third showed upregulation in air-cured K326 and Cuba1 and flue-cured Basma, with downregulation in the other five sample groups. 

In this study, targeted fatty acid metabolomics analysis was conducted on K326 plants at the flowering stage (55 days after transplanting), with 15 tissue samples analyzed. A total of 11 short-chain fatty acids (SCFAs) were measured. Among the 15 tissue samples, the contents of SCFAs in flower buds (FB) and upper leaves (UL) were higher than in other tissues, except for 3-Methylvalerate and 2-Methylvalerate. Cluster analysis (Fig. b) revealed that the samples were grouped into five clusters, with FB and UL forming one group, exhibiting significantly higher relative content compared to other tissues. Flower buds and upper leaves (UL) clustered together, and the 11 measured SCFAs were grouped into four categories: 2-Methylvalerate and 3-Methylvalerate, which showed lower levels in flower and upper tissues compared to other parts; 4-Methylvaleric acid, 2-Methylbutyrate, and Isobutyric acid, which were higher in flower tissues compared to other parts. The total content of short-chain fatty acids in the leaves followed the pattern: UL > U2L > LL > ML > L2L (Fig. c), and in the stems: US > MS > LS. The total short-chain fatty acid content in the upper leaves, stems, and flower buds was higher than in the middle and lower parts of the plant. This suggests that the short-chain fatty acid levels in K326 at the flowering stage were generally higher in the upper parts of the plant compared to the middle and lower sections, although the root did not show the lowest content. In the targeted fatty acid metabolomics analysis, 50 fatty acids were measured, including all but gamma-Linolenic acid, across the various tissue samples (C8:0-C22:6). The fatty acid clustering heatmap (Fig. d) indicated that the tissue samples could be grouped into three categories: roots as one group, flower buds, all leaves, and the main veins of upper leaves as another group, and stems and other main leaf veins as the third group. Significant differences in fatty acid content were observed between the upper, middle, and lower leaf, vein, and stem tissues, with distinct separation based on fatty acid content. The root had higher levels of nine fatty acids, including Octadecanoic acid, Hexadecanoic acid, and Heptadecanoic acid, compared to other tissues. Flower buds exhibited elevated levels of Linoelaidic acid, Linoleic acid, and cis-13,16-Docosadienoic acid, compared to other tissues. Upper leaves (UL) contained higher levels of Caprylic acid, cis-10-Heptadecenoic acid, and trans-10-Pentadecenoic acid compared to other tissues. The specific content data for short-chain fatty acids and fatty acid metabolomics measurements are provided in Table . 

The short-chain fatty acid content in the leaves at two stages of the seedling bed period was measured. Among the 11 short-chain fatty acids detected, Acetic acid had the highest content in all tobacco varieties, accounting for 78–95% of the total short-chain fatty acids, with a higher proportion in the (transplant stag) ML stage compared to the (Seedling stage) S stage (Table ). In terms of total short-chain fatty acid content, K326 had the highest content in the S stage, followed by Basma, Cuba1, and Samsun. In the ML stage, Basma showed the highest content, closely followed by K326, while Cuba1 and Samsun had lower levels. For 3-Methylvaleric acid, Cuba1 had the highest proportion (0.18%) at the S stage, compared to the other three varieties. At the ML stage, Basma and Samsun had significantly higher proportions (0.18% and 0.16%, respectively) compared to K326 and Cuba1 (0.04%). In the case of Isobutyric acid, Cuba1 exhibited a lower proportion at the S stage but had a significantly higher proportion at the ML stage than the other varieties. For the C5 isomer/normal ratio, Cuba1 had the highest value at the S stage, while the ratio in the ML stage was lower than that of the other varieties. The C6 isomer/normal ratio in K326 was notably lower in the ML stage compared to the other varieties. At the S stage, K326 showed higher levels of differential short-chain fatty acids compared to the other three varieties. At the ML stage, Basma and Cuba1 exhibited downregulation of 2-Methylbutyrate compared to K326. Specific fold changes (FC > 1.2 or FC < 0.833 and P -value < 0.05) for the differential short-chain fatty acids are provided in Table . As shown in Fig. a, principal component analysis (PCA) based on PC1 (54.03%) and PC2 (21.09%) clearly differentiated the short-chain fatty acid profiles of Cuba1 and Basma. In the MS stage, PCA (Fig. d) with PC1 (36.46%) and PC2 (21.08%) showed that Samsun and Basma had similar profiles, with reduced differences in the first principal component compared to the S stage. Figure b, e shows the specific content of 11 short-chain fatty acids at two stages. According to the heatmap of normalized short-chain fatty acid content (Fig. c, f), most short-chain fatty acids in K326 at the S stage were higher than in the other varieties. Comparison between stages revealed that Acetic acid, Isovaleric acid, and 2-Methylbutyrate were higher in the ML stage, while Propionic acid, Isobutyric acid, Butyric acid, and Valeric acid showed the opposite trend (Fig. g). 

Differential gene expression was analyzed in tobacco leaves from four varieties at two stages of the seedling bed period: seedling stage (S) and mature leaf stage (ML). A total of 158.25 GB of clean data was obtained from 28 samples, with an average of 6.59 GB per sample. The average total mapping rate to the reference genome was 93.52%. At the S stage, compared to K326, Samsun showed 1973 differentially expressed genes (DEGs), with 906 upregulated and 1067 downregulated (DESeq2 padj ≤ 0.05, |log2FoldChange|≥ 1.0). Basma exhibited 3433 DEGs, with 1212 upregulated and 2221 downregulated. Cuba1 showed 2705 DEGs, with 961 upregulated and 1744 downregulated (Fig. a). At the ML stage, Samsun compared to K326 had 3383 DEGs, with 1969 upregulated and 1414 downregulated. Basma had 5020 DEGs, with 2161 upregulated and 2859 downregulated, while Cuba1 showed 2574 DEGs, with 1308 upregulated and 1266 downregulated (Fig. b). Basma showed the highest number of DEGs at both stages among the three varieties compared to K326. PCA indicated distinct separation of S and ML stages along PC1, with a more scattered distribution in the ML stage compared to S. Intra-variety comparisons in the ML stage showed greater variability (Fig. c). Venn diagram analysis of the common DEGs between K326 and the other three varieties revealed 710 DEGs at the S stage and 879 at the ML stage, with 457 shared DEGs between the two stages. This DEG set represents the major differences between K326 and the other varieties during the seedling bed period (Fig. d). GO enrichment analysis showed significant enrichment in “signal transduction” and “signaling” pathways (padj < 0.05, Fig. f). KEGG pathway analysis indicated significant enrichment in pathways such as “Various types of N-glycan biosynthesis,” “N-Glycan biosynthesis,” and “RNA degradation” ( p -value < 0.05, padj = 0.52, Fig. e). In the fatty acid-related biosynthesis pathways, genes involved in “Fatty acid degradation,” “Fatty acid biosynthesis,” and “Fatty acid metabolism” ( p -value < 0.24, padj = 0.52–0.61) were enriched. Three genes ( Nta01g31980 , Nta08g22780 , and Nta23g11140 ) were identified in these pathways, annotated as acyl-activating enzyme, alcohol dehydrogenase, and enoyl-[acyl-carrier-protein] reductase. These genes were not expressed in K326 at either stage, but were expressed in the other three varieties. Based on gene expression data, a fatty acid pathway heatmap for the four varieties during the seedling bed period was generated, showing pathways such as “Fatty acid biosynthesis” (nta00061), “Fatty acid elongation” (nta00062), and “Fatty acid degradation” (nta00071) (Fig. ). In the “Fatty acid biosynthesis” pathway (nta00061), the relative expression levels of 13 ACC genes, 2 MT genes, and 4 enoyl-[acyl-carrier-protein] reductase I genes were significantly higher in the S stage compared to the ML stage. Four KAS III genes and six TE genes showed half of the genes with higher relative expression in either the S or ML stage. Varieties showed varying expression patterns across different genes. For example, Nta23g00340 (KAS II) was not expressed in the ML stage, and Nta12g11030 (acyl-[acyl-carrier-protein] desaturase) was not expressed in the S stage. In the “Fatty acid elongation” pathway (nta00062), two palmitoyl-protein thioesterases and two palmitoyl-CoA hydrolases had higher relative expression in the ML stage compared to the S stage, while two long-chain 3-oxoacyl-CoA reductases showed higher relative expression in the S stage. K326 showed higher relative expression than the other varieties for these genes. Nta06g11470 and Nta05g12800 were not expressed in the ML stage. In the “Fatty acid degradation” pathway (nta00071), the Nta12g16270 and Nta11g17320 were not expressed in the ML stage for all varieties. Two acetyl-CoA C-acyltransferase genes showed significantly higher expression in K326 at the ML stage compared to the other varieties and stages. Two long-chain fatty acid omega-monooxygenase genes had higher expression in K326 during the S stage compared to the other varieties and stages. 

The short-chain fatty acid (SCFA) content in the middle leaves of four varieties was measured at the budding stage (55 days post-transplantation) and 10 days after topping (70 days post-transplantation). Among the 11 measured SCFAs, the total content at the budding stage was ranked as Basma > Samsun > Cuba1 > K326, while the ranking reversed after topping. The proportion of 3-methylvaleric acid was significantly higher in Basma and Samsun than in Cuba1 and K326 both before and after topping. The trend of isovaleric acid proportion was consistent at the budding stage, but after topping, all four varieties showed similar trends. The C5 isomer/normal ratio was higher in Basma and Samsun than in Cuba1 and K326 at the budding stage, while after topping, K326 and Basma had significantly higher ratios than Samsun and Cuba1 (Table ). Differential metabolite analysis (FC > 1.2 or FC < 0.833, P -value < 0.05) and within-group significance analysis revealed that at the budding stage, the primary differentially expressed SCFAs were valeric acid and hexanoic acid, which showed similar trends across varieties. After topping, Basma and Cuba1 had significantly higher levels of 3-methylvalerate than K326. Detailed information on the differentially expressed SCFAs across varieties is provided in Table . As shown in Fig. a, the PCA plot based on PC1 (83.7%) and PC2 (7.15%) clearly clustered the samples into two groups: the K326 and Cuba1 groups were closer to each other, while the Samsun and Basma groups were closer. After topping, K326 and Cuba1 were more distinctly separated, while Basma and Samsun clustered together (Fig. d). Figure c, f shows the specific content of 11 short-chain fatty acids at two stages. The contents of acetic acid, isobutyric acid, butyric acid, valeric acid, isovaleric acid, 2-methylbutyrate, hexanoic acid, and 2-methylvalerate were significantly lower after topping (Fig. g). After topping, K326 showed higher levels of all fatty acids compared to other varieties, except for 2-methylvalerate, propionic acid, and 3-methylvalerate. Samsun had significantly higher levels of 2-methylvalerate, propionic acid, and 3-methylvalerate compared to the other varieties. Cluster heatmap analysis (Fig. b, e) revealed that after topping, the short-chain fatty acid content patterns in Basma and Samsun became more similar compared to those before topping. Additionally, the content of most short-chain fatty acids in K326 after topping was higher than in the other varieties. The PCA plot for the budding stage (Fig. h), based on PC1 (54.43%) and PC2 (15.68%), showed that K326 was clearly separated from the other three varieties along PC1. After topping, Basma and Cuba1 clustered together, while K326 and Samsun grouped similarly. The two groups were distinctly separated along PC2. Among the 49 measured fatty acids, Basma, Samsun, and Cuba1 had 13, 16, and 5 differentially expressed fatty acids compared to K326 at the budding stage, respectively, and 6, 13, and 22 differentially expressed fatty acids after topping (Fig. i). At the budding stage, Samsun showed a downregulation of fatty acids compared to K326, with pentadecanoic acid showing the highest downregulation (Log2FC = -2.09). Basma also showed a downregulation of fatty acids compared to K326, except for linoelaidic acid, with trans-10-heptadecenoic acid being the most downregulated (Log2FC = -2.07). Cuba1 showed the highest downregulation of cis-11,14,17-eicosatrienoic acid (Log2FC = -2.02) compared to K326. Common differentially expressed fatty acids between Basma and Samsun compared to K326 included pentadecanoic acid, palmitoleic acid, and heptadecanoic acid. After topping, Samsun showed downregulation of 5 fatty acids compared to K326, with caprylic acid being the most downregulated (Log2FC = -0.62). Basma downregulated 12 fatty acids compared to K326, with cis-5,8,11,14,17-eicosapentaenoic acid showing the highest downregulation (Log2FC = -1.64). Cuba1 downregulated 21 fatty acids, with pentadecanoic acid showing the highest downregulation (Log2FC = -1.30). The common differentially expressed fatty acids between Basma, Samsun, and Cuba1 compared to K326 included decanoic acid and hendecanoic acid. Based on the clustering heatmap of the 49 measured fatty acids (Fig. j), the fatty acid content differences were classified into three groups: Basma and Cuba1 after topping, Basma and Samsun at the budding stage, and the remaining four groups. 

Differential gene analysis was performed on the leaves of four tobacco varieties at the squaring stage. A total of 78.04 G of clean bases were obtained from 12 samples, with an average of 6.50 G per sample. The average total alignment to the reference genome was 94.45%. Principal component analysis (PCA) revealed that Basma and Samsun were closely related, while K326 and Cuba1 were clearly distinct (Fig. a). In the differential gene analysis, 18,600 DEGs were identified between Samsun and K326, with 9,067 upregulated and 9,533 downregulated (DESeq2 padj ≤ 0.05 |log2FoldChange|≥ 1.0). For Basma vs. K326, 17,770 DEGs were detected, including 8,520 upregulated and 9,250 downregulated. For Cuba1 vs. K326, 15,537 DEGs were found, with 6,678 upregulated and 8,859 downregulated (Fig. b). Using K326 as the control, differential gene analysis was conducted for the other three varieties, and a Venn diagram was used to identify a common set of 9,456 DEGs across the two stages, which represents the major differential gene set for K326 vs. the other varieties at the squaring stage (Fig. c). Enrichment analysis of this gene set showed that 37 GO terms were significantly enriched (padj < 0.05). Notably, the terms “sequence-specific DNA binding,” “unfolded protein binding,” and “calcium ion binding” were most enriched (Fig. d). KEGG analysis identified 10 pathways with significant enrichment (padj < 0.05), with the most significant being “MAPK signaling pathway – plant,” “Photosynthesis—antenna proteins,” and “Endocytosis” (Fig. e). In the fatty acid-related synthesis and metabolism pathways, 35 genes, including KAT1 , KCS1 , AIM1 , and FAD2 were enriched in pathways such as “Fatty acid degradation,” “Fatty acid biosynthesis,” and “Fatty acid metabolism.” Clustering analysis of these 35 genes revealed similar expression patterns in Basma and Samsun, distinct from K326 and Cuba1. Notably, genes such as Nta02g06150 ( KAT1 ), Nta12g29320 ( KAT1 ), Nta12g07350 ( KCS11 ), and Nta11g08380 ( KCS11 ) showed significantly higher expression in K326 compared to the other varieties (Fig. f). Six differentially DEGs were selected for qRT-PCR validation of their expression levels, including ACOT8 ( Nta02g19100 ), AL7B4 ( Nta05g24430 ), KCS11 ( Nta16g12070 ), KCS20 ( Nta20g17200 ), KCS4 ( Nta08g20960 ) and LACS8 ( Nta17g22980 ). These genes are DEGs from four fatty acid-related metabolic pathways (Fig. f).The expression data obtained by RNA-Seq and qRT–PCR were compared. The genes related to the fatty acid pathway ACOT8 , AL7B4 , KCS4 , LACS8 had exactly the same trends as the transcriptome results (Fig. ). 

Short-chain fatty acids and fatty acids were analyzed in the middle leaves of four tobacco varieties following two different curing methods. The PCA plot based on 11 short-chain fatty acids and principal components PC1 (43.51%) and PC2 (30.95%) clearly differentiated between the curing methods (Fig. a). Except for the air-cured Cuba1 (Cuba1_AC), the second principal component effectively separated the air-cured and flue-cured samples, with air-cured samples clustering at the top and flue-cured samples at the bottom. The Cuba1_AC sample exhibited a more scattered distribution in the principal components. Figure b shows the number of differential short-chain fatty acids both between varieties and within varieties under different curing methods. Generally, air-cured (AC) samples of the same variety exhibited fewer downregulated or unchanged short-chain fatty acids compared to flue-cured (FC) samples, indicating higher short-chain fatty acid content in air-cured leaves. The differential short-chain fatty acid clustering heatmap (Fig. c) revealed that Samsun_AC and Basma_AC clustered together, as did Samsun_FC and K326_FC, while Basma_FC and Cuba1_FC formed another cluster. Overall, the short-chain fatty acids were grouped into two categories: Butyric acid, Hexanoic acid, Isobutyric acid, Valeric acid, and 2-Methylbutyrate clustered together, showing upregulation in Samsun_AC, Basma_AC, and K326_AC, and downregulation in the remaining samples. Figure d presents the trends in short-chain fatty acid content before and after curing. With the exception of Propionic acid, 2-Methylvalerate, and 3-Methylvalerate, all other short-chain fatty acids showed varying degrees of increase after curing. For K326, the content of short-chain fatty acids (except for Propionic acid and 2-Methylvalerate) followed the order: air-cured > flue-cured > pre-curing. Similar trends were observed for Isobutyric acid, Butyric acid, the three C5 acids, and Hexanoic acid in Basma and Samsun. Additionally, the changes in the three C5 acids before and after curing were consistent across all four varieties. A targeted fatty acid metabolomics analysis was performed on the processed tobacco leaves, identifying a total of 47 fatty acids. A PCA plot (Fig. a) based on PC1 (51.45%) and PC2 (17.51%) showed a significant distinction in fatty acid composition between air-cured and flue-cured tobacco leaves along the first principal component, with the flue-cured samples being more tightly clustered. Differential analysis of the fatty acids revealed that, under air-curing, Samsun, Basma, and Cuba1 had 21, 8, and 13 differential fatty acids compared to K326, respectively. Under flue-curing, Samsun, Basma, and Cuba1 showed 7, 13, and 5 differential fatty acids compared to K326. Notably, air-cured Samsun and Basma had more downregulated differential fatty acids compared to K326, while flue-cured samples showed the opposite trend (Fig. b). When comparing different curing methods for the same variety, K326 showed 2 upregulated and 15 downregulated differential fatty acids under flue-curing compared to air-curing. The most downregulated fatty acid was Tridecanoic acid (log2FC: -2.66). Four fatty acids, including Tridecanoic acid, cis-10-Heptadecenoic acid, Docosanoic acid, and Heneicosanoic acid, showed a decrease (log2FC < -1) across all varieties. In contrast, gamma-Linolenic acid and Arachidonic acid showed upregulation (log2FC > 0.5). For air-cured varieties, Samsun, Basma, and Cuba1 showed three shared downregulated differential fatty acids compared to K326, namely alpha-Linolenic acid, cis-11,14,17-Eicosatrienoic acid, and homo-gamma-Linolenic acid. Fatty acid clustering analysis (Fig. c) based on differences in fatty acid content showed that Samsun_AC and Basma_AC clustered together, while K326_AC and Cuba1_AC formed another group. The remaining samples each formed separate clusters. Overall, the measured fatty acids were divided into three categories: one category showed downregulation in flue-cured compared to air-cured, another showed upregulation, and the third showed upregulation in air-cured K326 and Cuba1 and flue-cured Basma, with downregulation in the other five sample groups. 

Tobacco growth can be broadly divided into two stages: the seedling stage and the field stage. The seedling stage spans from sowing to transplanting the seedlings into the field, while the field stage starts from transplanting and ends at harvest. Each stage can be further subdivided into smaller periods . Among the four key stages selected in this study, the seedling transplanting stage was conducted under stable conditions in a controlled artificial climate chamber. This provided more consistent light and temperature conditions compared to tobacco seedlings grown in greenhouses under actual production settings. Topping is a critical event in the field stage of tobacco, as there are significant changes in the main substances of the plant before and after topping . Therefore, both pre- and post-topping periods were chosen for analysis. Notably, Basma and Samsun varieties flower much earlier than K326 and Basma. According to the results of this study, after topping, the content of short-chain fatty acids in the middle leaves of K326 decreased, as did the levels of five long-chain fatty acids. Decanoic acid and hendecanoic acid were identified as the major fatty acids differentiating K326 from the other three varieties. Among the 49 fatty acids measured, there were no significant differences in the fatty acid content between Basma and Samsun before topping. This trend was also reflected in gene expression, where principal component analysis of differentially expressed genes between Basma and Samsun showed no major differences, and the number of differentially expressed genes was much lower compared to other varieties. After topping, Basma fatty acid content was similar to that of Cuba1. The synthesis of saturated fatty acids in tobacco primarily occurs in the cytoplasm, although fatty acid synthesis also takes place in the mitochondria and chloroplasts. The main precursor for synthesis is acetyl-CoA, which is converted into malonyl-CoA by acetyl-CoA carboxylase (ACCase). ACCase is one of the key rate-limiting enzymes in the de novo fatty acid synthesis pathway in plants, and in tobacco, ACCase exists in its prokaryotic form . Key enzymes in the fatty acid synthesis pathway include citrate lyase (ACL), acetyl-CoA carboxylase (ACCase), β-ketoacyl-ACP synthase I (KAS I), and acyl carrier proteins (ACP). In the seedling stage of this study, three genes Nta01g31980 , Nta08g22780 , and Nta23g11140 were enriched in the fatty acid biosynthesis pathway. These genes were annotated as acyl-activating enzyme, alcohol dehydrogenase, and enoyl-[acyl-carrier-protein] reductase, respectively. Expression of these genes was not detected in K326 during both periods, while they were expressed in the other three varieties. In the gene expression differences at the budding stage, genes such as Nta02g06150 (KAT1), Nta12g29320 (KAT1), Nta12g07350 (KCS11), and Nta11g08380 (KCS11) showed relatively higher expression levels in K326 compared to the other three varieties. This difference in gene expression may contribute to the distinct fatty acid profile in K326 during the corresponding periods. In long-term production and experimental research, it has been observed that during the growth and development of tobacco plants, as well as the maturation and senescence of tobacco leaves, the content of various fatty acids initially increases and then decreases. As the leaves mature and age, the content of saturated fatty acids decreases, while the content of unsaturated fatty acids gradually increases. Spatially, the fatty acid content is higher in the middle and lower leaves of the plant, and lower in the upper leaves. Liu et al. found that fatty acids and similar compounds decrease during the vegetative growth phase, but increase during the reproductive growth phase . In the upper leaves, the fatty acid content peaks around 75 days post-transplantation, just before early flowering, and continues to decrease in senescent leaves . These results are consistent with the observation in this study that the fatty acid content in K326 decreases after topping, which is likely due to the shift back to vegetative growth following bud formation. During the growth, development, and senescence of tobacco leaves, the total fatty acid content first increases and then decreasesz . During the air-curing process, the fatty acid content in tobacco leaves also decreases . According to the results of this study, the variation in fatty acid content in tobacco leaves can be attributed to differences in varieties, leaf positions, and processing methods. Therefore, when curing tobacco leaves, it is recommended to adjust the curing conditions and processing measures based on the fatty acid content to optimize leaf quality and flavor. 

Our study systematically analyzed and compared the fatty acid composition in four different tobacco varieties. Using metabolomics and transcriptomics, we performed differential analyses on tobacco leaves at four developmental stages and under two curing methods. A fatty acid metabolic profile was established for the K326 variety at the budding stage, and three fatty acid-related genes that were significantly different from those in other varieties were identified during the seedling stage. Gene expression profiles for three fatty acid-related pathways were constructed for this period. Additionally, we compared the fatty acid differences in tobacco leaves across the four varieties under flue-curing and air-curing methods. This study broadens our understanding of fatty acids in different tobacco varieties and developmental stages, providing a molecular basis for further research on fatty acids and the genes involved in their biosynthesis and metabolism. 

Supplementary Material 1: Table S1. Raw Data of Targeted Fatty Acid Metabolomics for K326 (Unit: ng/g). Supplementary Material 2: Table S2. Short-Chain Fatty Acids Information Table at the Seedling Stage. Supplementary Material 3: Table S3. Short-Chain Fatty Acids Information Table Before and After Topping. Supplementary Material 4: Table S4. qPCR Primer Information. Supplementary Material 5: RNA-seq analysis detailed methods. Supplementary Material 6: Table S5. DEG and pathway enrichment information 

