To investigate the role of NtGA3ox1, encoding a Gibberellin 3-beta-dioxygenase 1, in tobacco flowering, sgRNA (AATTACCCGAATCCCATGCATGG) was fused into the genome editing vector (pBUE-2gRNA-ZH) with T7E1 restriction sites (Xing et al., 2014). The recombinant vector was transformed into Agrobacterium GV3101, and genetic transformation of tobacco was performed in accordance with the leaf disc method (Gallois and Marinho, 1995). Next, CRISPR/Cas9-mediated NtGA3ox1 knockout mutant was obtained, and the floral phenotypes and plant architecture were evaluated. A late-flowering mutant HG14 was screened with significant delayed flowering and impaired plant development. Genomic DNA was extracted using the Plant Genomic DNA Kit (Tiangen Biotech, Beijing, China). Specific primer pairs (P-GSP5 + P-GSP3) were used to amplify sgRNA target sites and listed in Supplementary Table 1 . After purification, PCR product was linked to pEASY-T1 vector (Transgen, Beijing, China). The mutation site was identified by Sanger sequencing.

The ecotype Honghuadajinyuan (Hongda) and HG14 were floated in the greenhouse, with temperature at 26°C and the photoperiod was 16 h light/8 h darkness with 60% relative humidity. Tobacco seedlings with 35-40 days old were transplanted to the field. All the plants were field evaluated in Shilin (Yunnan province, 103°30’E, 24°81’N) in 2021 and 2022, and Xichang (Sichuan province, 102°06’E, 27°52’N) in 2022. The field experiment was conducted in a randomized block design, with 3 replicates, 2 rows in each plot, 60 plants for each replicate, with plant spacing of 0.5 m and row spacing of 1.2 m. Field management was carried out according to the field trial scheme.

We next detailed measured the critical growth period, such as budding and flowering stage, and morphological characteristics of Hongda and HG14. Budding stage: the date when blossom buds of a plant were fully visible. Flowering stage: the date when 50% of the test plants bud emergence and flowering, respectively. Natural plant height: the distance from the ground-level base to the first fruit. Topping plant height: the distance from the ground-level base to the top of the stem. It was measured when removing the flower and the two guard leaves after flowering. Stem girth: the circumference between the fifth and sixth leaf during the first green fruiting stage. Leaf number: the number of leaves on the plant after removing the bottom leaves. Leaf length and width: the biggest leaf in the center was chosen. The width was determined by measuring the length vertically between the main vein and the widest section of the leaf surface. Leaf length was measured from the stem joint to the leaf tip on the front side of the leaf. For each trait, at least 15 plants from each replication were measured.

Research on structural domains is important for protein functional analysis. Gene-specific primers P-GSP5 and P-GSP3, developed based on the acquired cDNA fragment, were used to get the full-length sequence of the NtGA3ox1. The primers were listed in Supplementary Table 1 . Then the amino acid sequence of NtGA3ox1 was used as query for BLAST alignment with E-value < 0.1 in NCBI (http://blast.ncbi.nlm.nih.gov/Blast.cgi), and similar amino acid sequences were retrieved by BLASTP. The phylogenetic tree was built using the neighbor-joining method in MEGA version 7 (http://www.megasoftware.net/), with 1000 bootstrap replicates (Kumar et al., 2016).

There is close relationship between the tissue expression pattern of a gene and its function. The temporal and spatial expression pattern of NtGA3ox1 were analyzed by qRT-PCR after NtGA3ox1 was demonstrated to be the functional gene controlling HG14 late flowering. The root, stem, leaf, and organ tissues from Hongda at the seedling, rosette, and flowering stages were harvested, respectively. Three individual plants were mixed as one replicate, three replicates for each. It was enveloped in foil film and frozen in liquid nitrogen immediately. Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, USA). RNA concentration and quality were determined using a Nanodrop 2000 instrument (Thermo, Waltham, USA). cDNA was synthesized by one-step gDNA Removal and cDNA Synthesis SuperMix (Genstar, Beijing, China) according to manufacturer’s instructions. qRT-PCR was performed on an ABI StepOnePlus instrument (Takara, Dalian, China) using an SYBR Premix Ex Teq II kit (Takara Bio., Kyoto, Japan). The RT-PCR program was as follows: melting at 95°C for 15 s, and amplification with 40 cycles of 95°C for 15 s and 60°C for 30 s. The relative expression level of the NtGA3ox1 was analyzed according to the 2^­ΔΔCT method. NtActin was used as reference gene. The primer sequences were listed in Supplementary Table 1 .

The leaves of HG14 and wild type Hongda at the budding, flowering, and mature stages were selected based on the spatial and developmental features of NtGA3ox1 gene expression. Six replicates for each. RNA was isolated using the TRIzol reagent (Invitrogen, Carlsbad, USA), and their qualities were determined using the Nanodrop2000, respectively. cDNA libraries were constructed using the NEBNext^® Ultra™ RNA Library Preparation Kit for Illumina Sequencing System (NEB, USA). The transcriptome sequencing was conducted by Wuhan Huada Gene Technology Co. Ltd. Following the FastQC quality assessment, clean reads were generated and matched to the tobacco reference genome (ftp://ftp.solgenomics.net/genomes/Nicotiana_tabacum/, Edwards et al., 2017) using Bowtie 2 (v2.3.4.3).

Phyper (https://en.wikipedia.org/wiki/Hypergeometric_distribution) was utilized to perform GO (http://www.geneontology.org/) and KEGG (https://www.kegg.jp/) enrichment analysis based on hypergeometric testing, with a cutoff of Q-value ≤ 0.05. In order to further validate the transcriptome sequencing results, DEGs enriched in plant hormone signal transduction, plant-pathogen interaction, GA and amino acid biosynthesis were chosen for gene expression verification by qRT-PCR. qRT-PCR was conducted as previously described, with three independent biological replicates and three technical replicates, primers were designed by qPrimerDB (https://qprimerdb.biodb.org/, Lu et al., 2018) and listed in Supplementary Table 1 . NtActin was employed as the reference gene.

For metabolome sequencing, the same samples for transcriptome analysis were employed. A freshly collected plant leaf was immediately frozen in liquid nitrogen and grounded into powder (30 Hz, 1 min). The leaf samples were weighed exactly at 50 mg and subsequently dissolved in 1 mL of MeOH/water/formic acid (15:4:1, v/v/v). After vortexing for 15 min, the mixture was centrifugated for 10 min (12000 r/min, and 4°C). The supernatant was transferred into sterile plastic microtubes and dried by evaporation. The residue re-suspended in 500 mL water containing 3.5% formic acid and 1 mL ethyl acetate, following by vortexing for 15 min, and centrifuged for 5 min (12000 r/min, and 4°C). 500 mL ethyl acetate was added to the re-extracted sample before the supernatant was mixed.

The combined ethyl acetate component was evaporation dried and then dissolved in ACN. 10 mL triethylamine and 10 mL 3-bromopropyltrimethylammonium bromide were added to the resultant solution. The reaction mixture was vortexed, heated at 90°C for 1 h, evaporated to dryness under nitrogen gas stream, redissolved in aqueous ammonia solution (90%), and then filtered through a 0.22 m membrane filter for further liquid chromatography mass spectrometry (LC-MS/MS) analysis. The targeted gibberellin metabolome sequencing was carried out by Metware biotechnology Co. Ltd. (Wuhan, China). A total of 18 substances including gibberellin A1 (GA1), GA15, and GA19, etc. ( Supplementary Table 2 ) were examined, with four biological replications for each.

Alkaloids, amino acids and other chemicals determine the economic value of tobacco leaf (Troje et al., 1997), which directly affects the cigarette quality. Based on the transcriptome profiling results, common DEGs enriched in Linoleic acid and fatty acid metabolism pathway. Therefore, alkaloids, chlorine, organic acids, amino acids, and total sugars in HG14 and Hongda were analyzed to explore the effect of NtGA3ox1 on chemical components. All samples were dried in an oven at 60°C for 30 min and then grounded thoroughly using a whirl wind grinding instrument. All testing procedures were in accordance with the Tobacco Industry Standard of People’s Republic of China. The content of routine chemical components of fresh tobacco leaves at the maturity stage were detected by a flow analyzer. The determination of routine chemicals, primary aromatic constituents, and inorganic nutrients were detected according to different analytical methods ( Table 1 ), respectively. In this study, each line was tested with three independent biological repeats, with three plants for each.

Statistical significance between wild type Hongda and HG14 was determined by one-way analysis of variance (ANOVA) or Student’s t-test in Excel. For DEG identification, a false discovery rate-adjusted p value (FDR) < 0.05 and fold change >4, DESeq2 (v1.4.5) was utilized. Significantly regulated metabolites between groups were considered at *P < 0.05 or **P < 0.01 by two-pairs Student’s t-test.

GA regulates a number of metabolic pathways to govern plant growth and development. To investigate the potential roles of NtGA3ox1 in tobacco growth and development, we generated NtGA3ox1 knock-out mutants by CRISPR/Cas9. After comparison with the wild type Hongda, a late-flowering mutant HG14 was identified, which exhibits significantly delayed flowering and stunting development, such as dramatically reduced biomass, shortened plant height, and increased glandular hairs ( Figure 1A ). According to the Sanger sequencing results, all clones had a T insertion in NtGA3ox1 target site, gaining a pre-stop ( Figure 1B ). Stable transgenic lines were used for sub-sequent analysis.

The agronomic characteristics of Hongda and HG14 were further evaluated in 2021 and 2022 under field conditions. Similar results were obtained from Shilin in 2021 and Xichang 2022 ( Figure 2 ). The development of HG14 continuously lagged behind that of wild type Hongda, and the deficiency grew as development progressed ( Figure 2A ). Take the results of 2022 in Shilin as an example, HG14 took 13.0 and 12.1 days longer than Hongda to reach the budding and flowering stages, respectively, at 72.4 and 83.8 days as opposed to 59.4 and 71.7 days ( Figure 2B ). HG14 has shorter elliptical leaves, narrower stem and leaf angles, thin main veins, and lower growth potential. However, leaf number did not significantly alerted ( Figure 2B ). As an example, plant height of HG14 was 106.4 cm, a reduction of 35.4% from Hongda’s (164.7 cm) in Shilin 2022 ( Figure 2B ). As for the topping plant height, it was much shorter for HG14 at 52.5 cm compared to 96.9 cm for Hongda (P < 0.01). When compared to Hongda, HG14 had significantly shorter and narrower leaves, averaging around 58.9 cm as opposed to 72.2 cm in length (P < 0.01), and 21.4 cm as compared to 32.4 cm in width (P < 0.01, Figure 2C ). Additionally, stem girth of HG14 was significantly smaller than Hongda, averaging 9.3 cm as compared to 12.2 cm (P < 0.01, Figures 2B, D ). In summary, it can be predicted that NtGA3ox1 should play a positive role in tobacco growth and development.

We cloned the NtGA3ox1 coding sequence from Hongda and searched the Pfam database (http://pfam.xfam.org) using the NtGA3ox1 protein sequence. As shown in Figure 3A , NtGA3ox1 contains two functional domains, non-haem dioxygenase N-terminal domain (44-146, DIPR026992) and Fe(II) 2-oxoglutarate dioxygenase domain (197-298, PS51471). In plants, Fe(II) 2OG dioxygenase domain enzymes catalyze the formation of plant hormones, such as ethylene, gibberellins, etc.

Gibberellin dioxygenase is a large gene family in plants, governs the bio-synthesis of GA and crucial to the whole life cycle of plants. Therefore, phylogenetic tree was built utilizing NtGA3ox1 and homologs of rice, Arabidopsis, wheat, and maize, etc. Gibberellin dioxygenase genes in plants exhibit evolutionary conservatism, as shown by their close relationship and sharing of two aforementioned dioxygenase domains. Interestingly, the phylo-analysis revealed that Nicotiana tabacum, Arabidopsis, and Gossypium hirsutum were grouped together, whereas Oryza sativa and Zea maize were clustered ( Figure 3B ), indicating the evolutionary distinctions between food and other crops.

We next investigated the expression pattern of NtGA3ox1 across tissues and growth stages in Hongda. Lower expression levels were observed in stem compared to root and leaf during the seedling stage ( Figure 4 ). Relatively higher expressions were detected in stem and leaf at the rosette stage. Flower has the highest expression level among the tested tissues, which further evidenced that NtGA3ox1 involved in controlling plant flowering. The expression level in root gradually declined during the course of development, but it was dramatically raised in stem, which might control the flow from root to stem and thereby manipulate plant growth.

36 cDNA libraries were created and sequenced to comprehensively assess the transcriptome alterations in Hongda and HG14 at budding, flowering, and mature stages. All clean reads were compared and matched to tobacco reference genome v4.5 (https://solgenomics.net/organism/Nicotiana_tabacum/genome), with clean reads Q20 > 98.11% and Q30 > 94.07% ( Supplementary Table 3 ), respectively. Pearson correlation analysis was used to assess the consistency of gene expression across samples, and high repeatability among biological replicates were observed for all time-points.

Volcano maps were created to visualize the distribution of DEGs in each group. As shown in Figure 5A , the red, green, and grey dots indicate the up-regulated, down-regulated, and non-DEGs, respectively. A heatmap cluster analysis was conducted to illustrate the results of biological replicates, rows correspond to genes and columns to samples ( Figure 5B ). Those genes with similar expression patterns were grouped into four major clusters. Among the three stages, the smallest number of DEGs was observed in FS-14 (flowering stage-HG14) vs. FS-CK (flowering stage-Hongda) (2147), with 1071 up-regulated and 1076 down-regulated, respectively ( Figure 5C ). In contrast, the largest number of DEGs was found in MS-14 (mature stage-HG14) vs. MS-CK (mature stage-Hongda) (6202) including 2289 up-regulated and 2278 down-regulated. There were 4449 DEGs, including 2604 up-regulated and 1845 down-regulated, in BS-14 (budding stage-HG14) vs. BS-CK (budding stage-Hongda). Furthermore, FS-14 vs. FS-CK had the fewest DEGs in comparison to the other two groups, illustrated that although Hongda and HG14 had distinct flowering times, most gene expression patterns were comparable because of their common genetic ancestry. The Venn diagram ( Figure 5D ) shows the overlap of DEGs among the three groups. Totally, 957 identical DEGs between BS and FS, 839 identical DEGs between MS and FS, and 1670 identical DEGs between MS and BS.

The 4449 DEGs in BS-14 vs. BS-CK were classified into three functional groups by GO word enrichment analysis: biological process, cellular component, and molecular function ( Figure 6 ). The two largest groupings in the biological process, as determined by the GO analysis, were “cellular process” and “metabolic process”. The DEGs were mostly engaged in cell, cell portion, membrane, etc. in the category of cellular component. “Binding” and “catalytic activity” genes predominated in the enrichment of molecular function. For the other two comparisons, MS-14 vs. MS-CK and FS-14 vs. FS-CK, the same was true. The KEGG pathway enrichment analysis revealed that a significant portion of DEGs for BS-14 vs. BS-CK group were mainly enriched in few pathways, including one plant-pathogen interaction pathway, three signaling pathways (plant hormone signal transduction, MAPK, and phosphatidylinositol), and porphyrin metabolism ( Figure 7 ). Plant-pathogen interaction pathway and porphyrin metabolism were also discovered in the FS-14 vs. FS-CK group. When comparing the MS-14 vs. MS-CK group, those previously identified KEGG pathways were included. Other metabolism pathways, such as biosynthesis of monobactam, metabolism of carbon, fructose and sucrose, starch and sucrose, were also enriched.

Totally, 436 comparable DEGs were identified among the three groups, implying that they play a continuous role throughout NtGA3ox1-regulated tobacco leaf development. 108 DEGs were up-regulated with loss function of NtGA3ox1 in all three stages ( Figure 8A ), and 30 DEGs were found to be enriched in known KEGG pathways ( Figure 8B ) - diterpenoid biosynthesis (KO00904), plant-pathogen interaction (KO04626), ribosome (KO03010), and plant hormone signal transduction (KO04075). Furthermore, 73 down-regulated DEGs were obtained in HG14 and mostly enriched in KEGG pathways, including those for plant-pathogen interaction and plant hormone signal transduction ( Figures 8C, D ), which corresponded to the up-regulated DEGs and these show their diverse roles in plant growth and associations with Gibberellin metabolism.

We further selected 5 genes in plant hormone signal transduction (Nitab4.5_0000008g0310, Nitab4.5_0000174g0290, Nitab4.5_0000245g0250, Nitab4.5_0000295g0090, and Nitab4.5_0002093g0090), 4 genes in plant-pathogen interaction pathway (Nitab4.5_0000036g0500, Nitab4.5_0000045g0460, Nitab4.5_0000108g0150, and Nitab4.5_0000258g0120), 3 genes in GA metabolism (Nitab4.5_0000076g0060, Nitab4.5_0000124g0040, and Nitab4.5_0001573g0060), 2 genes in amino acid biosynthesis pathway (Nitab4.5_0000412g0050 and Nitab4.5_0001769g0010), and 1 gene in organic acid metabolism pathway (Nitab4.5_0000137g0140) to verify the transcriptome sequencing results at the flowering stage. FPKM for these genes across the three time periods were presented in Supplementary Figure 1 , and the differential fold changes ranging from 0.09 to 68.76. The qRT-PCR results suggested that most genes were down regulated and showed consistent expression pattern in Hongda and HG14 with knock out of NtGA3ox1 ( Figure 9 ).

We focused on the gibberellin metabolites accumulation variation in Hongda and HG14 throughout various leaf development stages. At the budding and flowering stages, the overall content of the 15 detected gibberellins in HG14 was much lower than that in Hongda, whereas greater at the mature stage as a result of the large rise in GA8 in HG14 ( Figure 10 ). As shown in Figure 10A , knock out of NtGA3ox1 resulted in a substantial reduction in the levels of three GA metabolites, including a 50.22% reduction in GA12-aldehyde (P < 0.05), a 58.91% reduction in GA9 (P < 0.05), and a 58.39% reduction in GA29 (P < 0.01). And a substantial 53.91% rise in GA6 (P < 0.05) was also observed. Only GA9, transiently elevated at 180.62%, was discovered in the metabolic profiles of GAs from the two genotypes at the flowering stage ( Figure 10B ). As shown in Figure 10C , three GAs metabolites, GA20, GA7, and GA15, were found to be highly down-regulated at the mature stage, with a decrease of 42.56%, 33.05%, and 71.49%, respectively. Additionally, knocking out NtGA3ox1 caused the downregulation of GA51, GA4 and GA29 as well as leaf aging.