CRISPR/Cas9-mediated mutagenesis of VvbZIP36 promotes anthocyanin accumulation in grapevine (Vitis vinifera)

Abstract Anthocyanins are plant secondary metabolites that have a variety of biological functions, including pigmentation. The accumulation of anthocyanins is regulated by both transcriptional activators and repressors. Studies have shown that bZIP family members act primarily as positive regulators of anthocyanin biosynthesis, but there are few reports of negative regulation. Here, we report that a grapevine (Vitis vinifera) bZIP gene from group K, VvbZIP36, acts as a negative regulator of anthocyanin biosynthesis. Knocking out one allele of VvbZIP36 in grapevine with CRISPR/Cas9 promoted anthocyanin accumulation. Correlation analysis of transcriptome and metabolome data showed that a range of anthocyanin biosynthesis genes were activated in VvbZIP36 mutant plants relative to the wild type, resulting in the accumulation of related metabolites, including naringenin chalcone, naringenin, dihydroflavonols, and cyanidin-3-O-glucoside. The synthesis of stilbenes (α-viniferin), lignans, and some flavonols (including quercetin-3-O-rhamnoside, kaempferol-3-O-rhamnoside, and kaempferol-7-O-rhamnoside) was significantly inhibited, and several genes linked to their metabolism were downregulated in the mutant plants. In summary, our results demonstrate that VvbZIP36 is a negative regulator of anthocyanin biosynthesis that plays a role in balancing the synthesis of stilbenes (α-viniferin), lignans, flavonols, and anthocyanins.


Introduction
Flavonoids are a family of specialized plant metabolites that includes flavonols, flavones, proanthocyanidins (PAs; also called condensed tannins), and anthocyanins [1]. They accumulate in different organs and tissues of a wide range of plants and have a variety of biological functions [1]. For example, flavonols and PAs provide protection from UV radiation and fungal invasion, respectively [2], and anthocyanins contribute pigmentation to flowers and fruits, thereby attracting insects and animals and promoting pollen and seed dispersal [3]. Flavonoids have also been associated with enhanced resistance to abiotic and biotic stresses [4]. There is increasing evidence that dietary anthocyanins act as powerful antioxidants with beneficial effects on human health, such as the prevention of cardiovascular disease, obesity, diabetes, and cancer [5].
The biosynthesis of flavonoids is controlled by a complex network of structural and regulatory genes. Many transcriptional activators involved in regulating flavonoid biosynthesis, especially anthocyanins and PAs, have been identified in various species. These include AtMYB113 and AtMYB114 (A. thaliana) [14], MdBBX22 and MdHY5 (bZIP family genes from apple; Malus domestica) [15], FvbHLH9 (strawberry) [16], PpMYB10 and PpbHLH64 (pear) [17], PpNAC1 and PpMYB10.1 (peach) [18], and VvMYBA1, VvMYBA2, and VvMYBPA1 (grapevine) [11,19]. It has also been determined in multiple species that a transcriptional activation complex of R2R3-MYB, bHLH, and WD-repeat (WD40) proteins regulates the expression of flavonoid biosynthetic genes [4]. In addition to transcriptional activators, transcriptional repressors have been shown to play important roles in regulating flavonoid biosynthesis pathways. For example, heterologous overexpression of strawberry FaMYB1 in Lotus corniculatus inhibits the expression of LAR and ANR and significantly reduces the content of PAs in leaves [20]; the poplar (Populus) MYB182 protein significantly reduces the content of anthocyanins and PAs in poplar by reducing the expression of DFR1, ANS1, and ANR1 [21]; and PpMYB18 acts as a negative regulator, balancing the accumulation of anthocyanins and PAs during peach fruit ripening and at juvenile stages of development [13].
Grapevines (Vitis vinifera) are economically important fruit-producing plants rich in flavonoids [22]. The regulation of flavonoid biosynthesis by transcription factors (TFs) has been studied in grapevine, including analyses of VvMYBPA1/VvMYBPA2, VvMYBA1/VvMYBA2, VviMYB86, VvBBX44, and VvWRKY26 [1,11,19,[23][24][25]. For example, VvMYBPA1 or VvMYBPA2 activate the expression of VvANR and VvLAR and promote PA accumulation [1,19]. Overexpression of VvMYBA1 activates the transcription of Vv3AT in the white grape variety Chardonnay, causing the fruit to accumulate anthocyanins and turn red, whereas silencing of VvMYBA1 or VvMYBA2 in the red grape variety Shiraz results in a significant decrease in fruit anthocyanin content [11]. Notably, bZIP family genes, especially HY5, have also been shown to play key roles in the regulation of flavonoid biosynthesis in many species, including Arabidopsis [26], apple [27], tomato [28], pear [29], and strawberry [16]. However, there are only a few related studies of bZIP TFs in grapevine [30,31]. Characterization of the regulation of grapevine flavonoid biosynthesis has mostly centered on gain-offunction experiments, and there are few loss-of-function studies. Thus, much remains to be learned about the transcriptional regulation of flavonoids by bZIP TFs in grapevine, and in this regard, the analysis of new functional knockouts of genes of interest represents a promising experimental strategy.
In the current study, we used CRISPR/Cas9 to target the grapevine bZIP family gene VvbZIP36 and found that it contributes to the regulation of anthocyanin levels by balancing the synthesis of stilbenes (α-viniferin), lignans, flavonols, and anthocyanins.

Bioinformatic analysis of VvbZIP36
In previous studies, 47 bZIP TF-encoding genes were identified from grapevine (V. vinifera) [40]. Tissue-specific expression analysis showed that the bZIP gene VvbZIP36 has high expression in leaves, stems, flowers, tendrils, and especially fruits [40]. In this study, we found that VvbZIP36 is expressed in the roots, tendrils, stems, leaves, and fruit of grapevine, with the highest expression in mature stems, leaves at fruit set, and senescent leaves, as well as ripening fruit (Fig. 1A), similar to a previous study. According to a phylogenetic analysis using bZIP TF sequences from Arabidopsis and rice, we found that VvbZIP36 was classified in the K subgroup as a homolog of AtbZIP60 and OsbZIP74, which encode proteins with both a bZIP domain and a trans-membrane domain (TMD) [40]. Previous studies showed that both AtbZIP60 and OsbZIP74 play important roles in the endoplasmic reticulum (ER) stress response in Arabidopsis and rice, respectively [41,42]. To further analyze the group K bZIP TFs, we searched the National Center for Biotechnology Information (NCBI) (http://blast.ncbi.nlm.nih.go v/Blast.cgi) database for VvbZIP36 homologs and found 25 sequences from other species, all of which contained the bZIP domain and a putative TMD (Fig. S1). When analyzing their phylogenetic relationships, as shown in Fig. 1B, we found that VvbZIP36 was most closely related to Camellia sinensis bZIP18 and Ipomoea nil bZIP60.

Knocking out one allele of VvbZIP36 promotes the accumulation of anthocyanins in grapevine
In our previous study, we successfully knocked out one allele of VvbZIP36 in grapevine (Thompson Seedless) using the CRISPR/Cas9 genome editing system and obtained one mutant grapevine plant (KO45) [43]. No off-target events were identified in KO45, as determined  by whole-genome sequencing (WGS) [43]. Subsequently, through continuous identification of regenerated lines, 7 transgenic positive plants were identified from 22 regenerated plants (Fig. S2). Among them, the mutant line KO21 was found by Sanger sequencing. As shown in Fig. 2A, we designed four CRISPR/Cas9 sgRNAs for the first exon of VvbZIP36 [43], but only target 1 was mutated in KO21 and KO45. In addition, the KO21 and KO45 lines were found to contain monoallelic mutations, with a base insertion in KO21 and a base deletion in KO45.
We observed that the leaves of KO21 and KO45 turned red (Fig. 2B) and hypothesized that VvbZIP36 plays a role in regulating anthocyanin biosynthesis. We therefore collected leaves from the mutant lines (KO21 and KO45) and the wild type (WT) to determine the total anthocyanin contents. As shown in Fig. 2C, D, the anthocyanin levels in the leaves of both KO21 and KO45 were significantly higher than those in the WT.

Anthocyanin biosynthetic genes are differentially expressed between KO45 and WT plants
To further analyze the connection between VvbZIP36 and anthocyanin biosynthesis, we performed an RNAseq analysis of KO45 (no off-target events certified by WGS [43]) and WT leaves from plants grown in a growth chamber under normal conditions. The results of normalized cluster analysis and Pearson's correlation analysis show that the biological replicates of each group were highly correlated (Pearson coefficient > 0.9; Fig. S3). A total of 4103 DEGs were identified between KO45 and the WT, including 1965 upregulated and 2138 downregulated genes (Fig. S4). As shown in Fig. 3A, the gene ontology MapMan analysis revealed that most of the downregulated DEGs were related to "photosynthesis", "tetrapyrrole", and "cell wall". By contrast, DEGs involved in "lipid synthesis", "fatty acid synthesis", "amino acid metabolism", and "reactive oxygen" were both upregulated and downregulated. Finally, most of the DEGs related to phenolic compound metabolism were upregulated. This information is summarized in Supplementary Data S1.

Metabolome profiling indicates that VvbZIP36 regulates anthocyanin biosynthesis in grapevine
To further investigate the function of VvbZIP36, we compared the metabolomes of the leaves of transgenic (KO45) and WT plants. The results of normalized cluster analysis and Pearson's correlation analysis show that the biological replicates of each group are highly correlated (Pearson coefficient > 0.9; Fig. S6). A total of 347 differentially abundant metabolites between KO45 and the WT were identified and quantified; 185 were more abundant in KO45 than in the WT, and 162 were less abundant (Fig. S7). Of the top 20 metabolites showing higher levels in KO45, most were flavones, triterpene, flavanols, anthocyanidins, organic acids, and stilbenes, whereas the 20 that showed the greatest reduction were mainly flavonols, lignans, organic acids, amino acids and derivatives, and phenolic acids ( Table 2).

VvbZIP36 is important in anthocyanin biosynthesis regulation
To further understand the VvbZIP36-associated regulatory network involved in modulating anthocyanin biosynthesis, we performed a correlation analysis of the transcriptome and metabolome data. The expression of genes in the anthocyanin biosynthesis pathway (five PAL, three CHS, two CHI, two F3 H, two F3H, one DFR, and one ANS gene) was generally higher in KO45 than in WT, and the corresponding metabolites, such as naringenin chalcone, naringenin, dihydroflavonols, and anthocyanins, were more abundant to varying degrees (Fig. 6).
In other branches of the flavonoid biosynthesis pathway, such as flavanones and flavones, the levels of most metabolites were significantly higher in KO45 than in the WT. Correspondingly, the expression of VvF3 5 H2 (flavonoid 3 ,5 -hydroxylase 2: VIT_08s0007g05160) and two Type-I flavone synthase genes (VIT_16s0098g00860, VIT_13s0047g00210) was higher. The PAs and most of their precursor metabolites (flavanols) also accumulated to higher levels in KO45 than in the WT, accompanied by higher expression levels of VvLAR1 (leucoanthocyanidin reductase 1: VIT_01s0011g02960) and VvANR (anthocyanidin reductase: VIT_00s0361g00040). The content of some flavonols was higher in KO45 than in the WT, whereas the abundance of others was lower (Table 2 and Fig. 6), accompanied by higher expression of three FLS genes (VvFLS1: VIT_09s0002g08090; VvFLS3: VIT_ 08s0105g00380; VvFLS5: VIT_18s0001g03430) and lower expression of two FLS genes (VvFLS2: VIT_18s0001g03510; VvFLS4: VIT_03s0091g01080) and two FLR genes (flavonol-3-O-rhamnosyltransferase: VIT_06s0009g02010, VIT_ 06s0009g01990).
Interestingly, all the differentially abundant lignan metabolites were lower in KO45 than in the WT, whereas five PAL genes from the early general phenylpropanoid pathway were expressed at higher levels in KO45 (Fig. 6). Notably, we did not identify DEGs involved in lignan biosynthesis in the KEGG pathway enrichment analysis. In the stilbene biosynthesis pathway, two ROMT genes (trans-resveratrol di-O-methyltransferase: VIT_12s0028g01940, VIT_15s0045g01490) were identified,   one of which showed higher expression in KO45 than in the WT, whereas the other showed lower expression. Subsequently, we randomly selected 15 of the 30 DEGs involved in the phenylpropanoid pathway from Fig. 6 for confirmatory qRT-PCR expression analysis. As shown in Fig. 7, with the exceptions of VvF3H2 and VvLAR1, their expression patterns were consistent with the transcriptome data.

Discussion
The basic region/leucine zipper (bZIP) TFs are important regulators in plants and are involved in regulating many metabolic processes, such as energy metabolism, development, senescence, light signaling, anthocyanin biosynthesis, and response to biotic and abiotic stimuli [26,44]. In previous studies, 47 bZIP TFs were identified from grapevine [40]. Tissue specific-expression analysis showed that VvbZIP36 was highly expressed in leaves, stems, flowers, tendrils, and especially fruits [40]. Similar results were found in this study (Fig. 1A). Because previous reports indicate that bZIP TFs play an important role in the synthesis of anthocyanin [26][27][28][29][30][31], we speculated that VvbZIP36 might be involved in the regulation of anthocyanin biosynthesis. Knocking out one allele of VvbZIP36 did not affect expression ( Fig. 2A, Fig. S9) but caused 46.45% of the transcripts of this gene to be mutated (Table S1). This resulted in a significant accumulation of anthocyanins and a red leaf phenotype in white grapevine (Thompson Seedless) (Fig. 2B-D), suggesting that VvbZIP36 is a negative regulator. Previously, most bZIP genes reported to be involved in regulating anthocyanin biosynthesis have been described as positive regulators, including MdbZIP44 and MdHY5 in apple [27,45], SlHY5 in tomato [7,28], and PybZIPa and PyHY5 in pear [10,29]. Similar results were previously reported from studies of grapevine: VvibZIPC22 was observed to promote the accumulation of several flavonoids by regulating the expression of CHS and FLS [30], and VvHY5 promotes flavonol biosynthesis by activating the expression of an FLS [31]. These results indicate that VvbZIP36 may have a different regulatory framework from other bZIP TFs that control the anthocyanin biosynthesis pathway. Previously, most studies on K subgroup genes (including AtbZIP60 and OsbZIP74) have been related to endoplasmic reticulum stress [41,42]. This is the first report to show that K subgroup genes (such as VvbZIP36) participate in plant anthocyanin biosynthesis.
To further analyze the regulatory mechanism by which VvbZIP36 influences anthocyanin biosynthesis, we performed transcriptome analysis of the leaves of a CRISPR/Cas9-mediated VvbZIP36 mutant line (KO45) and WT plants. A total of 1965 and 2138 genes were expressed at higher and lower levels, respectively, in KO45 compared with the WT (Fig. S4). Genes expressed at lower levels in the transgenic leaves were mainly involved in photosynthesis, cell wall and lipid metabolism, and the fatty acid synthesis pathway (Fig. 3), suggesting that VvbZIP36 may positively regulate these synthetic pathways, consistent with similar studies. For example, in Arabidopsis, AtHY5 promotes photosynthetic pigment accumulation in response to light [46]. In soybean, GmbZIP04g and GmbZIP07g can transcriptionally regulate the expression of GmRCAα, a key photosynthetic gene, by binding to its promoter [47]. Moreover, overexpression of NsbZIP1 and HSbZIP1 significantly promotes the synthesis of lipids in Nannochloropsis salina [48] and Chlorella sp. HS2 [49], respectively. In addition, research shows that MabZIP93 can activate the expression of cell wall modifying genes during banana (Musa) fruit ripening [50]. Given this association of VvbZIP36 with photosynthesis and the central role that photosynthesis plays in plant growth and development [51], it is perhaps not surprising that we only obtained transgenic VvbZIP36 knock-out grapevine lines with single allele mutations ( Fig. 2A, Table S1): complete VvbZIP36 knockout may be lethal.
We identified an association between the DEGs that were expressed at higher levels in KO45 and the flavonoid metabolism pathway (Fig. 3), consistent with the observed phenotype (Fig. 2). Further analysis revealed that this subset included genes in the core anthocyanin biosynthesis pathway, including PAL, CHS, CHI, F3H, F3 H, DFR, FLS, and ANS (Table 1).
Metabolome analysis showed that a total of 185 metabolites were more abundant in KO45 than in the WT, and 162 were less abundant (Fig. S7). Those that were less abundant belonged to the lipid, organic acid, lignan, and amino acids and derivatives categories, whereas those that were more abundant belonged to the terpenoid, stilbene, tannin, and flavonoid categories (Fig. 4,5). Interestingly, we found that most of the genes related to terpenoids were downregulated (Fig. 3A), but the corresponding levels of most triterpene metabolites were higher in KO45 than in the WT ( Table 2). We speculated that the metabolism and regulatory network of triterpenes in mutant plants may be very complicated, or some terpenoid-related genes may not have been discovered, leading to this result.
Stilbenes, lignans, and flavonoids are all synthesized from the same precursor metabolite, p-coumaroyl-CoA [52]. Subsequently, correlation analysis of the transcriptome and metabolome data showed that many genes involved in the anthocyanin biosynthesis pathway, including PAL, CHS, CHI, F3 H, DFR, and ANS were expressed at higher levels in KO45 than in the WT, and many related metabolites (including naringenin chalcone, naringenin, dihydroflavonols, and the anthocyanins cyanidin-3-O-(6 -O-acetyl)glucoside, cyanidin-3-O-(6 -O-p-coumaroyl)glucoside and peonidin-3,5-Odiglucoside) were more abundant in KO45 (Fig. 6). Naringenin chalcone, naringenin, and dihydroflavonol metabolites are all important intermediates in anthocyanin biosynthesis, which is part of the flavonoid pathway [2,12]. Our data suggest that VvbZIP36 can suppress the expression of key anthocyanin biosynthesis genes; however, Qiu et al (2019) found that knocking out SlHY5 (a bZIP gene) in tomato can reduce anthocyanin levels and expression of these anthocyanin biosynthesis genes [28]. Similar studies have also been performed on grapevine with the same results [30,31]. This suggests that bZIP proteins can operate via different mechanisms in regulating anthocyanin biosynthesis. Interestingly, we found that the abundance of the anthocyanin compound petunidin-3-O-(6 -O-caffeoyl)glucoside was lower in KO45 than in the WT (Fig. 5), indicating that VvbZIP36 may be involved in regulating the balance of different anthocyanin compounds. The identification of genes involved in petunidin-3-O-(6 -O-caffeoyl)glucoside synthesis will be an interesting area of future study.
We found that when the anthocyanin biosynthetic pathway was induced, the abundance of various flavonoid components (flavonols, flavanols, PAs, flavanones, and flavones) increased significantly in the mutant lines (Fig. 6). However, we also found that some flavonols, flavanols, flavanones, and flavones were present at lower levels in the mutant lines (Fig. 5), especially flavonols (quercetin-3-O-rhamnoside, kaempferol-3-Orhamnoside, and kaempferol-7-O-rhamnoside) ( Table 2). Correspondingly, the expression of four genes from the flavonol synthesis pathway (VvFLS2, VvFLS4, and two VvFLR genes) was lower in the transgenic plants. These results suggest that VvbZIP36 is involved in balancing the synthesis of flavonols and anthocyanins. Similar results have also been found in the stilbene (α-viniferin) and lignan biosynthetic pathways, which share the same precursor metabolite, p-coumaroyl-CoA, as the precursor of flavonoids, naringenin chalcone. We infer from our results that the action of VvbZIP36 causes more p-coumaroyl-CoA to be used to synthesize stilbenes (α-viniferin) and lignans, and more dihydroflavonols to be synthesized into flavonols (quercetin-3-O-rhamnoside, kaempferol-3-O-rhamnoside, kaempferol-7-O-rhamnoside, etc.), thereby restricting the synthesis of naringenin chalcone and anthocyanins, ultimately leading to relatively low levels of anthocyanin accumulation. By contrast, more anthocyanin accumulation was seen in mutant plants because of the inhibition of flavonol, stilbene (α-viniferin), and lignan synthesis (Fig. 8). That knockout of one allele of VvbZIP36 caused such a remarkable phenotype is surprising. A reasonable explanation is that VvbZIP36 plays a complicated role in the phenylpropanoid pathway. It is an inhibitor for anthocyanins but an activator for metabolites such as flavonols and lignans. Single allele mutations lead to weakened inhibition of anthocyanin synthesis and promotion of stilbene, flavonol, and lignan synthesis. Finally, this caused significant accumulation of anthocyanins in leaves. Second, this may be caused by a cumulative effect. Single allele knockout only partially inhibits or activates downstream genes, but these genes will continue to function for a period of time after expression. This is consistent with the phenotype we observed, in which leaf redness gradually deepened with increasing leaf age (Fig. 2B).
In conclusion, the results reported here indicate that VvbZIP36 contributes to the regulation of anthocyanin biosynthesis by balancing the synthesis of stilbenes (α-viniferin), lignans, flavonols, and anthocyanins in grapevine. This study provides new insights into the regulation of anthocyanin biosynthesis that will be valuable for grapevine breeding programs.

Phylogenetic and expression profile analyses
To investigate the relationship between VvbZIP36 and other bZIP TFs, multiple amino acid sequence alignments were generated using DNAMAN software (Version 5.2.2.0, Lynnon Biosoft, USA) with default parameters, and a phylogenetic analysis was carried out using the neighborjoining (NJ) method and MEGA software (version 11.0.8), as previously described [53].

Plant materials and growth conditions
Vegetative propagation of mutant lines 21 and 45 (KO21 and KO45) and WT grapevine (V. vinifera L. cv. Thompson Seedless) was performed as described in a previous study [54]. Two-month-old transgenic and WT plantlets were transferred to plastic pots (10 × 10 × 8 cm) containing a mixture of soil (Pindstrup, Denmark) and vermiculite (1: 1, v: v). After one month of adaptation, the plants were transferred to larger plastic pots (diameter, 15 cm and height, 20 cm) filled with a mixture of soil and vermiculite (1: 1, v: v) and cultivated in a growth incubator (25 • C, 16h photoperiod, and a light intensity of 200 μmol m −2 s −1 ; #G2-1000D3, Youke, China) for four months.

Detection of on-target mutations by Sanger sequencing
Identification of mutant lines was performed as previously described [37,43]. Specific primers (NPTII-F: 5 -AGA GGC TAT TCG GCT ATG ACT G-3 ; NPTII-R: 5 -CAA GCT CTT CAG CAA TAT CAC G-3 ) were used to identify stable transgenic lines. Mutant type was identified using VvbZIP36 gene-specific primers (VvbZIP36-Target-F: 5 -ATG GAC GAT TTG GAA ATT GGG G-3 ; VvbZIP36-Target-R: 5 -TCA CAC CAA AAC TCC ATG AG-3 ). The PrimeSTAR Max DNA Polymerase kit (Takara, China) was used to amplify potential mutation regions in the first exon of VvbZIP36, and the PCR products were inserted into the pClone007 Simple Vector (TSINGKE, China) to select single clones for sequencing.

Determination of total anthocyanin contents
The anthocyanin contents of leaves from three sevenmonth-old transgenic (KO21 and KO45) and WT plants, grown in a growth chamber (#G2-1000D3, Youke, China), were measured after sampling the third to fourth leaf from the top of the plants. The total anthocyanin content was calculated according to Li et al. (2012) [55]. In brief, leaves (approximately 0.3 g) were incubated in 5 mL of 0.1% (v/v) methanol-HCl in the dark at room temperature, and the absorbance of each extract was measured at 530, 620, and 650 nm with a spectrophotometer (UV2600; Unico, China). The experiment was repeated independently three times.

RNA-seq and data analyses
The third to fourth leaves from the tops of seven-monthold transgenic (KO45) and WT plants were collected, and their transcriptomes were characterized by RNAseq analysis. Samples from three plants per line were pooled and used as one independent sample, and three independent samples were analyzed per genotype. Total RNA was extracted from the transgenic and WT leaves using TRIzol reagent (Invitrogen, USA). The RNA library preparation and sequencing were performed by Metware Biotechnology Co., Ltd. (Wuhan, China). In brief, RNAseq libraries were generated using the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer's instructions before sequencing on the Illumina HiSeq platform. The original data were filtered using fastp (v0. 19.3) to remove reads with adapter sequences. The resulting cleaned reads were mapped to the V. vinifera reference genome (PN40024) (http://pla nts.ensembl.org/index.html) using HISAT (v2.1.0). FPKM was calculated using featureCounts (v1.6.2), and differential gene expression was determined using DESeq2 (v1.22.1) with a threshold of corrected P-value < 0.05 and absolute log 2 FC (fold change) > 1.

Metabolome analysis
The third to fourth leaves from the tops of three sevenmonth-old transgenic (KO45) and WT plants were collected for a widely targeted metabolome analysis. Approximately 5 g of leaf tissue per line was used as one independent sample, and three independent samples were analyzed. Sample extraction, metabolite identification, and quantification were performed by Metware Biotechnology Co., Ltd. (Wuhan, China). An ultra-performance liquid chromatography-electrospray tandem mass spectrometry (UPLC-ESI-MS/MS) system (UPLC, SHIMADZU Nexera X2) was used to analyze sample extracts as described previously [57]. Metabolites with a variable importance in projection (VIP) ≥ 1 and an absolute log 2 FC (fold change) ≥ 1 were considered to be differentially regulated between groups.

Quantitative real time (qRT)-PCR analysis
For gene expression analysis, RNA was extracted from the leaves of seven-month-old transgenic (KO21 and KO45) and WT plants grown in a growth chamber using the E.Z.N.A. Plant RNA Kit (#R6827-01, Omega Bio-tek, USA). First-strand cDNA was synthesized using the HiScript II Q RT SuperMix for qPCR(+g DNA wiper) Kit (#R223-01, Vazyme) according to the manufacturer's protocol. qRT-PCR was performed using the SYBR qPCR Master Mix (#Q311-02, Vazyme) on a StepOnePlus RT-PCR instrument (Thermo Fisher Scientific), and the gene-specific primers used are listed in Supplementary Table S2.

Statistical analysis
Data analysis was performed using Excel 2016 software (Microsoft Corporation, USA). Data were plotted using SigmaPlot (v14.0, Systat Inc., CA, USA). Tests for significant differences were performed using SPSS statistics software (v25.0, IBM China Company Ltd., Beijing, China).