Combined metabolic phenotypes and gene expression profiles revealed the formation of terpene and ester volatiles during white tea withering process

As one of the six major tea types in China with a long history and unique manufacturing methods, white tea has a delicate and subtle flavor that is light, sweet, and mildly floral, with a hint of fruitiness. It is often described as having a mellow taste with a smooth, silky texture. Studies have indicated that white tea has a range of biological activities, including free radical scavenging^[1], anti-inflammatory^[2], and anticancer^[3] activities. As a consequence, white tea has been widely accepted by international consumers in recent years, and the white tea industry has been developing rapidly. White tea has the fewest processing steps compared to other tea categories, consisting of just a long withering and drying period without any enzymatic inactivation or fermentation treatment, and withering is the core step in forming quality^[4].

Aroma is an important attribute of tea flavor and has long been a concern of scholars from various countries. Terpenoids and esters are essential volatile metabolites in tea^[5,6]. Of these, terpenoids are important secondary metabolites in plants, with well-defined characterized biosynthetic pathways and accumulation patterns in many horticultural crops^[7]. Ester volatiles, the main precursors of which are fatty acids, together with terpene volatiles have a pivotal influence on the aroma evaluation of fermented and semifermented teas^[6,8]. It has been shown that a large amount of terpene and ester volatiles accumulate during the processing of different tea types^[9−11]. Li et al. showed that 2-heptanol, linalool and geraniol contribute significantly to the aroma profile of white tea^[12]. Feng et al. suggested that the significant contributors to the sought-after elegant orchid aroma of tea are methyl jasmonate and (Z)-methyl epi-jasmonate^[13]. Although 116 terpenoids and 83 esters were detected in spring and fall white teas by Zhang et al., their effect on aroma formation is still not fully elucidated^[9]. Furthermore, the metabolites of terpenes and esters have been found to interact with each other in the form of phytohormone signaling, such as the JA and salicylic acid pathways to help plants withstand biotic or abiotic stresses throughout growth and development^[14,15]. During white tea processing, the endogenous enzymes in tea leaves greatly accelerate the oxidation, hydrolysis and degradation reactions of precursor molecules, thus enhancing the aroma formation rate and aroma concentration^[16]. A previous study has suggested that endogenous enzymes in tea leaves are involved in and improve the aromatic properties of white tea products, with a significant increase in grassy, floral and sweet odors after β-glucosidase treatment^[17]. In contrast, other studies have shown that the aromatic substances accumulated in tea leaves during the spreading period (short-term withering) are mainly volatile terpenoids and fatty acids, which may originate mainly from ab initio synthesis rather than through the hydrolysis of glycosides^[18]. In short, the molecular mechanisms of fragrance formation and regulation in white tea have attracted widespread attention, but the effect of withering on fragrance formation is still controversial, and the key genes in this process, including members of the known synthetic pathways, still require further identification and elucidation.

In this study, we explored the dynamic patterns and key time points of terpene and ester aromas throughout the white-tea withering process by a combination of gas chromatography-mass spectrometry (GC/MS) detection and sensory evaluation. Subsequently, RNA sequencing (RNA-Seq) was performed on tea samples at crucial time points to further reveal the molecular mechanisms associated with terpenoids and esters. WGCNA of differentially expressed genes (DEGs) associated with aromatic substance formation at the withering stage was completed to identify key regulatory genes. It is hoped that these results can provide a new perspective and theoretical basis for the optimization of white-tea withering.

Tea cultivar of 'Fuding Dahaocha' (Camellia sinensis (L.) O. Kuntze) was chosen for manufacturing white tea. Ninety kilogram (kg) of the harvested leaves consisting of one bud with two leaves were processed for the tea samples at the 'Xiangheshan' production base of the Pinpinxiang Tea Industry Co. Ltd., located in Fuding City, Fujian, China. Upon arrival at the base, the fresh tea leaves were split equally into three replicates (30 kg each) and spread out in the withering chamber using circular bamboo sieves with mesh (90 cm in diameter, about 0.5 kg of the leaves per sieve). The fresh leaves were withered at 30 °C for 30 h with an air humidity of 65%−85%, noting that the withered leaves could not be randomly removed during this process. Tea samples were uniformly extracted from the withering leaves at 0, 6, 12, 18, 24 and 30 h throughout the process and preserved frozen in liquid nitrogen, while 500 g of the samples were retained for freeze-drying followed by sensory evaluation. Meanwhile, we also monitored the moisture content and temperature of the leaves at each time-stage using the previously reported method^[19].

Quantitative and qualitative evaluation of volatile extracts

The detection of volatiles was performed by MetWare Company (Wuhan, China) operating on the Agilent 7890B-7000D platform following the previously published approach^[20]. In brief, the frozen tea samples were mixed into liquid nitrogen and ground into a powder. An aliquot (1.0 g) of the powder was quickly moved to a 20 mL head-space vial with NaCl saturated solution to inhibit any enzymatic activity. The volatile substance [3,4,5,6-2H4,OD]-2-Methoxyphenol (CAS No. 20189-11-1) was applied as internal standard (IS) to 50 µg/mL final concentration. For solid phase microextraction (SPME), the vials were closed by gland using a TFE-silicone head-space septum and then each vial was incubated at 100 °C for 5 min, followed by exposure of 65 µm divinylbenzene/carboxen/polydimethylsilioxan fibre to the head-space of the tea sample at 100 °C for 15 min. Upon exposure, the volatiles in the fibers were desorbed at the inlet of 7890B GC apparatus to enable analyzing and identifying the volatile compositions by comparison with the MS from the data system libraries (NIST & MWGC) and linear retention index based on a previous study^[21]. The volatile compound analysis was also conducted on quality control (QC) samples prepared by mixing 15 g of each sample powder in triplicate. We identified the volatiles that were substantially moderated among groups by variable importance projection (VIP) ≥ 1 and absolute log2(fold change) ≥ 1. The following formula was used to compute the relative amounts of the volatile constituents: Relative concentration (µg/g) = [(Peak area of target/Peak area of IS) × 0.1 µg]/ (1 g)^[18].

Sensory evaluation on the aroma of tea samples

Referring to the study by Li et al.^[12], the professional team consisted of six members (four males and two females with an average age of 30 years) from the field of tea science at Fujian Agriculture and Forestry University to complete the aroma evaluation of tea samples at each withering stage in accordance with the Chinese methodology for sensory evaluation of tea (GB/T 23776-2018). The team discussed and then selected five exact terms for the white tea aroma characteristics and determined the total scoring criteria. The intensity of the aroma contribution was rated on a scale of 0 to 6, with the higher the score, the greater the intensity (or the better the total sensory), where 0 = no or no perceived intensity, 3 = moderate intensity, and 6 = high intensity. Sensory evaluation was performed at room temperature (25 °C). Six white tea samples were offered in a randomized order and it was ensured that the sample numbers were not known to the panelists. Then, 3 g of each tea sample was brewed with 150 ml of boiling water. And the aroma characteristics of the tea samples after brewing were evaluated after 5 min.

RNA-seq and analysis

Total RNA of the tea samples was isolated according to manufacturer's directions with the TaKaRa MiniBEST Plant RNA Extraction Kit (TaKaRa, Japan). After passing quality control, the whole RNA was transferred to Metware Biotechnology Inc (Wuhan, China) for cDNA library creation and transcriptome sequencing. Twelve libraries (0, 12, 24, and 30 h × 3 biological replicates) were assayed on the Illumina Hiseq platform, and then clean reads were mapped to the reference genome 'Shuchazao' with HISAT2^[22]. Levels of gene expression were counted as fragments per kilobase of transcript per million pieces mapped (FPKM). DESeq2 was applied for differential expression analysis, and the gene with |log2 (fold change) | ≥ 1 and adjusted p-value < 0.05 was accepted as DEG.

WGCNA develops a regulatory network for co-expression

Using the WGCNA R package, a co-expression network was built by genes that score FPKM > 1 and variation coefficient > 0.5^[23,24]. The FPKM values for the genes screened at withering stages is normalized and an adjacency matrix was generated using the default formula. Referring to the conclusion of the PickSoftThreshold command, the β value representing the soft threshold to the matrix in was determined to be 24 in order to better present the connectivity of every node. (Supplemental Fig. S1a). The genes were further analyzed by hierarchical clustering, with dynamically cut the clustering results to obtain different gene modules (Supplemental Fig. S1b). To recognize the modules, minModuleSize was set to 30 and cutting height threshold of combo modules was adjusted to 0.25. And then we found the module most relevant to the ester volatiles content by calculating the correlation coefficients between the ester volatiles and the eigenvectors of each module. Additionally, gene significance was statistically utilized to reflect the relevance of gene expression with each trait. Cytoscape was used to visualize the module networks (Version 3.9.1, Cytoscape Consortium, Bethesda, MA, USA).

Quantitative real-time PCR (RT-PCR) examination

The PrimeScriptTM RT reagent kit with gDNA Eraser (TaKaRa, Japan) was employed to reverse transcribe mRNA into cDNA for the qRT-PCR analysis, and SnapGene 6.1 software was utilized to generate the gene-specific primers (Supplemental Table S1). RT-PCR experiments were conducted on a qTOWER3 Real-Time PCR System (Beckman Coulter, Inc., USA), and the PerfectStartTM Green qPCR SuperMix (TransGen Biotech Co., Beijing, China) was utilized (Analytik Jena, Germany). The GAPDH gene was employed as an internal control, and gene expression degrees were determined by means of the 2^–ΔΔCᴛ technique^[25].

Data analysis

Each tea sample had three independent biological replicates for GC/MS assay, RT-PCR and sensory evaluation. The identified volatiles' peak areas were normalized to the internal standard. Unsupervised principal component analysis (PCA) was completed with assistance from statistics function prcomp in R software (www.r-project.org, opened on 10 August 2021). The information on the fragrance properties of the volatile chemicals came from the literature and public sources (www.ncbi.nlm.nih.gov and www.thegoodscentscompany.com/search2.html, accessed on 10 August 2022). The Owuor's flavour index was calculated to assess variations in the fragrance quality of the tea samples^[26]. The Volcano, Venn, and heat maps were also conducted by means of statistics function ggplot2 within R software. Graph Prism 7.0 software (GraphPad, San Diego, CA, USA) was applied to conduct statistical analysis and draw the bar, line and pie graphs. Annotation of GO and KEGG were accomplished based on the previously published protocols^[19,27]. The tool (http://plantregmap.gao-lab.org/binding_site_prediction.php, accessed on 10 August 2022) was utilized to scanned for transcription factor binding sites from the upstream sequences of gene.

To explore the variation in aromatic volatiles with time, the aromatic volatile compound levels of fresh leaves were measured at 0, 6, 12, 18, 24 and 30 h of withering. We also measured the moisture content and temperature of the leaves to monitor the degree of withering. As the withering time extended, the fresh leaves gradually lost their luster, deepened in color and turned yellow‒green, and the leaf edges withered and curled (Fig. 1a). The moisture content measurements revealed that during the first 12 h of wilting of fresh leaves, the moisture decline was relatively gentle, with an average rate of decline of approximately 0.79%/h, while after 12 h, the water loss intensified, and the average rate of decline increased to 2.48%/h. However, as the moisture content of the fresh leaves continued to decrease, the leaf surface temperature increased. The temperature of the leaves increased significantly after 18 h of withering and continued until the end of the process (Fig. 1b). This finding shows that the withering of white tea is a simple but important process, i.e., the process of reducing the water content and increasing the temperature of tea leaves, which accumulates the material components to form the special flavor of white tea.

Figure 1. 

Tea plant leaves changes in water content, temperature and overview of nonvolatile com-pounds. (a) The phenotype of the picked tea leaves at 0, 6, 12, 18, 24, and 30 h. (b) The moisture content and temperature of tea leaves following withering for different time periods. Data represent the mean value ± standard deviation (n = 3), *, p < 0.05 compared with the sample of 0 h; **, p < 0.01 compared with the sample of 0 h. (c) The proportion of total content of various volatile substances in tea plant leaves during withering. (d) Principal component analysis score plot of tea leaves with different withering times (0, 6, 12, 18, 24, and 30 h). (e) KEGG classification of differential metabolites during withering.

All 288 differentially abundant volatile compounds with VIP ≥ 1 and absolute log₂(fold change) ≥ 1 were determined in the samples, consisting of terpenoids, esters, acids, alcohols, aromatics, aldehydes, hydrocarbons, ketones, and heterocyclic compounds (Fig. 1c). The total ion current (TIC) overlap diagram for MS detection of different QC samples and the distribution of the coefficient of variation (CV) for all groups indicated good reproducibility and reliability of the data recorded in this work (Supplemental Fig. S2a & S2b). QC samples were tightly clustered in the center of the PCA score plot, demonstrating the consistent reproducibility of the metabolite extraction and metabolomics analysis (Supplemental Fig. S2c). Then, PCA was applied to review the overall patterns of change in fresh leaf aroma components throughout white tea withering process. The score plot clearly showed stepwise alterations in volatile compounds from 0 to 36 h (Fig. 1d), indicating that the difference in white tea scent components became more significant with longer withering time. Since most of the volatiles were at low concentrations and the synthesis pathway of ester volatiles is intimately related to fatty acids, such as linolenic acid, we focused on the top three volatiles (Fig. 1c), including terpenoids (26.34% total volatiles), esters (24.16%) and fatty acids (20.31%). Several KEGG pathways of volatiles associated with terpene and ester volatiles were significantly enriched, e.g., monoterpenoid biosynthesis and alpha-linolenic acid metabolism. (Fig. 1e & Supplemental Table S2). Thus, the above results indicate that the synthesis and metabolism of terpene and ester volatiles are very active in leaves during wilting.

To reveal the variation pattern of terpene and ester volatiles throughout the processes, we further observed the content profiles of terpenes, esters, and fatty acids, including 17 terpenes, 36 esters, and seven fatty acids, and the concentrations of these volatiles presented significant variations in different stages of withering (Fig. 2a, Supplemental Table S3). In the initial stage of withering (0−12 h), the concentration of terpenoids in fresh tea leaves increased from 14.43 to 17.02 µg/g, esters from 12.12 to 16.68 µg/g, and fatty acids from 12.38 to 18.50 µg/g. After 18 h of withering, the concentration of terpenoids was 13.91 µg/g, that of esters was 11.69 µg/g and that of fatty acids was 6.92 µg/g. In the late stage of withering (24−30 h), the concentration of esters ranged from 15.47 to 21.58 µg/g, and fatty acids ranged from 12.39 to 17.47 µg/g, but there was no significant change in terpene concentrations. (Fig. 2b). The moderate withering treatment significantly increased the content of the three aroma types, but the extent of change varied, with the concentration of esters increasing the most throughout processing (78%), followed by fatty acids (41%) and terpenes (30.63%).

Figure 2. 

The variation rules of the main aroma substances and the assessment of aroma quality in tea leaves during the withering process. (a) Heat map of the types and contents of terpenes, esters, and acids volatiles detected in different samples. (b) Trends in volatiles of terpenoids, esters and fatty acids in tea leaves at different withering periods. (c) Aroma profile of different withering time points. (d), (e) The change in Owuor's flavor index during withering. The calculation of the Owuor's flavor index was based on the aroma compounds in Supplemental Table S3. Data represent the mean value ± SD (n = 3); Values with different letters (a–c) differ from each other significantly (different periods for the same component). **, p < 0.01 compared with the sample of 0 h.

Influence on aroma formation by volatile terpenes and esters

To investigate the effects of terpene and ester volatiles on scent formation in white tea, sensory evaluation was performed to evaluate the aroma grades of white tea with different withering times. Sensory evaluation is a strong method for identifying aroma characteristics and odor-related compoundsand is already generally used in white tea fragrance analysis^[28]. The main aroma characteristics of most terpene and ester volatiles are recognized as floral, fruity, woody and herbal in different tea categories^[29]. Therefore, we chose these four aroma characteristics to describe the aroma of withered leaves to determine the critical time when these two types of volatiles drive the aroma formation of white tea (Fig. 2c). The results of sensory evaluation suggested that the 0 h group had the strongest green aroma, close to the smell of a fresh leaf; and all samples after 6 h of withering had some floral aroma, while the 12 h group had the most prominent floral aroma; the 30 h group had the strongest fruity, woody and herbal aromas. 'Total' represents the integrated aroma features of the samples, and the higher its score is, the better its aroma. We found that 30 h of withering brought the best overall white tea aroma experience, followed by 24 h.

OFI was first applied to the aroma assessment of black tea, and subsequently, it has often been applied to assess the aroma characteristics of other tea types and positively correlates with tea tasters' results^[30]. OFI is the ratio of the total content of aromatic compounds with a sweet floral aroma (Group II) to the aromatic compounds with a grassy aroma (Group I); the higher the value of this ratio is, the better the aroma quality. To count OFI, we searched for the odor profiles of the terpene and ester aroma components found in this research and collected information on over half of them (Supplemental Table S3). As depicted in Fig. 2d, the sum of the aromatic components from Group II in withering tea leaves showed a significant increase at 12, 24 and 30 h compared to fresh leaves (0 h), with increases of 35.45%, 28.76% and 74.81%, respectively. It has been revealed that in Group II, methyl jasmonate (MeJA), geraniol, linalool, etc., are important terpenes and ester aromatic substances in tea that bring out strong floral and fruity aromas^[16]. The sum of aromatic compounds from Group I in fresh tea leaves did not vary significantly. The OFI fluctuated up and down but increased overall with withering time, with a peak at 30 h (Fig. 2e). Integrating the results including sensory evaluation and OFI suggested that the aroma evaluation of fresh leaves could be significantly improved by the withering process and that 12, 24 and 30 h were the critical timepoints at which the major aroma components (terpenes and esters) were significantly changed in withering leaves.

Transcriptome characteristics of white tea during withering.

Next, samples of 'Fuding Dahaocha' were analyzed by RNA sequencing at the four key timepoints (0, 12, 24, and 30 h) mentioned above. A total of 50,949,432 clean reads of 12 cDNA libraries were obtained, of which 75.88% were mapped to the 'Shuchazao' reference genome (Supplemental Table S4). The R² values (square of Pearson's correlation coefficient) indicated that the biological replicates of each point were reasonable and reliable (Supplemental Fig. S3a).

The numbers of DEGs increased with withering time (Supplemental Fig. S3b). However, compared with the 24–36 h stage, the rate of gene changes was greater in the 0–12 h and 12–24 h stages, and the number of downregulated genes was greater overall than the number of upregulated genes. (Fig. 3a). This result suggested that there was a dynamic response mechanism during withering that affected changes in gene transcription and that this mechanism was more active in regulating the transcriptome during the 0–24 h period. As there were more downregulated genes than upregulated genes, the former might also have essential effects on volatiles throughout processing, which is in line with previous findings^[28]. We then analyzed the number of overlapping DEGs in fresh tea leaves after 12, 24 and 36 h of withering, which totaled 3,644 overlapping DEGs (Fig. 3b). These overlapping DEGs converged mainly on two major types (biological process and molecular function) in the GO analysis results, including secondary metabolic process, serine-type carboxypeptidase activity, carboxypeptidase activity and phenylpropanoid biosynthetic process (Fig. 3c & Supplemental Table S5). Interestingly, the enrichment of DEGs in the serine-type carboxypeptidase activity and phenylpropanoid biosynthetic process pathways further indicated that the harvested leaves were under drought stress^[31]. Among the top 20 pathways in which the DEGs converged in the KEGG analysis results, several of them are involved in the biosynthesis of aroma precursors, including the biosynthesis of secondary metabolites, plant hormone signal transduction, sesquiterpenoid and triterpenoid biosynthesis, and alpha-linolenic acid metabolism (Fig. 3d & Supplemental Table S6a). This finding indicated that the biosynthesis of secondary metabolites was very active during the withering period and that the synthesis of terpene and ester volatiles might be regulated by drought stress and phytohormones.

Figure 3. 

Functional enrichment analysis of differentially expressed genes (DEGs). (a) The volcanic map was drawn by using the fold change in DEGs between the two groups of samples (including 0 h vs 12 h, 12 h vs 24 h, 24 h vs 30 h) and the p-value obtained by t-test, which was used to show the significant difference between the two groups. (b) Venn diagram of DEGs at different withering periods. (c) Top 20 significantly enriched GO pathways of the 3644 overlapping DEGs in (b). (d) Top 20 significantly enriched KEGG pathways of the 3644 overlapping DEGs in (a).

Expression patterns of key genes in pathways related to the synthesis of aroma components

According to the changes in aroma constituents observed during withering and the results of GO and KEGG pathway enrichment of the DEGs, we further examined the expression patterns of key genes associated with the synthesis pathway of aroma constituents, namely, the α-linolenic acid metabolism and terpenoid biosynthesis pathways (Fig. 4). Most key genes in the α-linolenic acid metabolic pathway except 13S-lipoxygenase (13-LOX) and allene oxide synthase (AOC), were differentially upregulated after withering, a trend in agreement with the changes in the content of ester volatiles (Fig. 4a). 13-LOX is the critical rate-limiting enzyme that catalyzes the oxidation of polyunsaturated fatty acids, and these gene members were able to maintain high activity in the early phase (0–12 h) of water stress, but were downregulated in the later phase (24–36 h), which may be related to the fact that lipid peroxidation is a damaging process^[32]. Thus, the contents of 3-hexen-1-ol acetate tended to increase, which may be affected by alcohol dehydrogenases (AHDs) and hydroperoxide lyase (HPL) (Fig. 4b). The aroma type of 3-hexen-1-ol acetate is sharp fruity-green and sweet, which may contribute to the fragrance formation of white tea. Of interest, among these signaling pathways, the AOS pathway is the only one capable of producing biologically active jasmonic metabolites^[33], and the expression level of the AOS gene (TEA001041) was significantly upregulated 57-fold after withering, indicating that the AOS gene has a major function in this process. Furthermore, the key rate-limiting enzyme in this pathway is acyl-CoA oxidase (ACX)^[34], the expression of which also increased. In addition, the expression of two jasmonic acid-amino synthetase (JAR) members (TEA009648 and TEA022080) increased markedly in late (36 h) withering, which may lead to the production of more jasmonoyl-isoleucine and thus increase the activity of monoterpenoid biosynthesis. Thus, during withering, phytohormone signaling associated with JA was positively regulated. The expression of the six jasmonate carboxyl methyltransferase (JMT) members (TEA015791, TEA032708, TEA031418, TEA028050, TEA028051, and TEA010054) also significantly increased over time and peaked at 12 h, which is nearly identical to the changes in methyl jasmonate (MeJA) content (Fig. 4b). Linolenic acid has been found to be oxidized autonomously or by drought stress to synthesize JA and its derivatives, which depend on its signaling to combat the stress^[35]. Moreover, jasmonic metabolites have a remarkable contribution to the aroma evaluation of tea leaves. For example, MeJA accounts for the strong floral aroma of oolong tea^[36]. The withering process might promote the aroma transformation of white tea by taking advantage of this same molecular mechanism in the harvested leaves.

Figure 4. 

Effect of apoptosis on gene expression in the α-linolenic acid metabolic pathway and terpene synthesis pathway and changes in related volatile content. (a) The expression patterns of key genes in α-linolenic acid metabolism pathway. 13-LOX, 13S-lipoxygenase; HPL, hydroperoxide lyase; ADH, alcohol dehydrogenase; AOS, allene oxide synthase; AOC, allene oxide cyclase; OPR3, 12-oxophytodienoate reductase 3; OPCL1, OPC-8:0 CoA ligase 1; ACX, acyl-CoA oxidase; MFP2, multifunctional protein 2; JMT, jasmonic acid car-boxyl methyltransferase, JAR, jasmonic acid-amino synthetase. (b) Changes in the content of partial fatty acid-derived volatile compounds during the withering process. (c) The expression regularities of key genes in terpene synthesis pathway. AACT, acetoacetyl-CoA thiolase; HMGS, 3-hydroxy-3-methylglutaryl-CoA synthase; HMGR, 3-hydroxy-3-methylglutaryl-CoA reductase; MK, mevalonate kinase; PMK, phospho-mevalonate kinase; MDD, mevalonate-5-diphosphate decarboxylase; IDI, isopentenyl diphosphate isomer-ase; PT, prenyl transferase; TPS, terpene synthase; DXS, 1-deoxy-D-xylulose 5-phosphate synthase; DXR, 1-deoxy-D-xylulose 5-phosphate reductase; MCT, 2-C-methyl-D-erythritol 4- phosphate cytidyl-transferase; CMK, 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase; MDS, 2-C-methyl-D-erythritol 2,4-cyclodiphosphate synthase; HDS, (E)-4-hydroxy-3-methyl-but-2-enyl diphosphate synthase; HDR, (E)-4-hydroxy-3-methyl-but-2-enyl diphosphate reductase. (d) The variation trend of total FPKM value of TPS family members at different withering time points. (e) Heat map of TPS family genes expression during withering. Values with different letters (a–c) differ from each other significantly; Data represent the mean value ± SD (n = 3).

The key genes of the terpene biosynthesis-related pathway showed various expression patterns during white tea withering (Fig. 4c). In the mevalonate (MVA) pathway used for sesquiterpene biosynthesis, two acetoacetyl-CoA thiolase (ACCT) members (TEA028475 and TEA013411), one 3-hydroxy-3-methylglutaryl-CoA synthase (HMGS) member (TEA007157), one 3-hydroxy-3-methylglutaryl-CoA reductase (HMGR) member (TEA012284) and one phosphomevalonate kinase (PMK) member (TEA022423) were markedly upregulated in expression after 30 h of withering. It has been indicated that mevalonate synthesis in plants is the first and most important stage in isoprene biosynthesis, and this reaction is mediated by HMGR (TEA012284). The overexpression of HMGR resulted in increased production of isopentenyl pyrophosphate (IPP)^[37]. Then, the expression of these genes in the methylerythritol phosphate (MEP) pathway of biosynthetic monoterpenes and diterpenes, including three 1-deoxy-D-xylulose-5-phosphate synthase (DXS) members (TEA012941, TEA012756, and TEA022165), one 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase (CMK) member (TEA033157), one 2-C-methyl-D-erythritol-2,4-cyclodiphosphate synthase (MDS) member (TEA027323), and one (E)-4-hydroxy-3-methyl-but-2-enyl diphosphate synthase (HDS)-related member (TEA010494), was also markedly upregulated by withering. One isopentenyl diphosphate isomerase (IDI)-related member (TEA000437) was similarly upregulated in expression. Since DXS (TEA012941) plays a critical regulatory role in the MEP pathway to catalyze isoprene biosynthesis^[38], increased expression of this gene promotes the synthesis of dimethylallyl diphosphate (DMAPP), a terpene precursor in tea leaves. Additionally, the JA signaling pathway regulates the MVA and MEP pathways by upregulating the expression of genes, such as HMGR and DXS, which encode key enzymes, leading to an increase in the production of IPP and DMAPP, which are then used for the biosynthesis of terpenes and other isoprenoids.

Moreover, Fig. 4d & e describe the expression pattern of terpene synthase (TPS). Up to 72 gene family members potentially associated with terpene biosynthesis have been found, but only 33 changed significantly during wilting, consistent with expression during green tea spreading^[18,39], indicating that the majority of TPS in postharvest leaves was not greatly affected. Statistics on their cumulative expression (total FPKM values) revealed that TPS gene expressions generally decreased after withering treatment, which also indicated that the biosynthesis of most terpenoids in postharvest tea leaves was significantly inhibited after 36 h of drought stress. This result might have also led to the weakening of floral aromas (of sensory evaluation) in white tea at 24–30 h of withering. However, we also observed significant increases in the expression of three neomenthol dehydrogenase-related members (TEA001376, TEA006904, and TEA011315), two farnesyl-diphosphate farnesyltransferase-related members (TEA000051 and TEA024261), and one terpene synthase 10 (TPS10, TEA009449) after withering. Linalool has been reported to be the major product of the TPS10 expressed at anthesis in Arabidopsis thaliana^[40]. The trend of TPS10 expression during withering was consistent with the changes in linalool content shown in Supplemental Table S3. In short, these TPS genes may exert an important function in the biosynthesis of terpene volatiles and thus promote the accumulation of aroma substances during white tea withering, and deserve further investigation.

To verify the data reliability of the RNA-seq transcriptome, we randomly selected 12 genes participating in the α-linolenic acid metabolism and terpenoid biosynthesis pathways for RT‒PCR analysis. In Supplemental Fig. S4, the RNA-seq data and RT‒PCR results both showed consistent expression patterns. This result indicates the reproducibility and reliability of the transcriptome data.

Expression patterns of transcription factors during withering

We further explored the transcription factors that may be involved in regulating aroma formation. From the 3644 DEGs at different withering periods (Fig. 3b), 64 transcription factors were identified by using high-quality binding motifs to scan for transcription factor binding sites in the promoter regions (2-kb sequences located just upstream of the transcription start site, p ≤ 1e-7) of the aroma synthesis-related genes mentioned above. All these transcription factors were differentially expressed in the samples at the four time points (0, 12, 24 and 36 h) (Fig. 5a). These transcription factors were further categorized into 24 gene families, including BES1-interacting Myc-like protein 2 (BIM2), B3 domain-containing proteins (B3), basic pentacysteine proteins (BPC), brassinosteroid signaling positive regulator family protein (BZR1), basic helix-loop-helix transcription factor (bHLH), double zinc finger protein (DOF), dehydration-responsive element-binding proteins (DREB), ethylene-responsive transcription factors (ERFs), LOB domain-containing protein (LBD), zinc finger protein (ZFP), and others (Fig. 5a & b). Figure 5a & c together depict the pattern of changes in these transcription factors and show that the cumulative expression (FPKM values) of transcription factors was significantly increased after 12 h of withering and that these upregulated transcription factors were mainly enriched in the ERF family.

Figure 5. 

Expression patterns of transcription factor families and co-expression modules identification based on weighted gene co-expression network analysis (WGCNA). (a) Heat map of TFs expression during withering. (b) Transcriptional regulation of α-linolenic acid metabolism and terpenoid biosynthetic genes. The color filled hexagons represent the structural genes, and the color filled spheres represent different transcription factors. The expression correlations between transcription factors and the structural genes are shown with colored lines (PCC, p ≤ 0.05), and the orange line indicates a positive correlation, and the blue line indicates a negative correlation. (c) The variation trend of total FPKM value of transcription factors at different withering time points. (d) Module-trait associations based on PCC. Each row represents one module, and the numbers of genes in each module were shown on the far left. Each column corresponds to one volatile metabolite. (e) protein co-expression network for the MEbrown module. The top 200 proteins based on topo-logical overlap measure (TOM) values were used to construct the coexpression network. (f) and (g) Protein co-expression net-works for the MEmagenta and MEpurple module. All proteins in the modules were used to construct the coexpression network.

Moreover, to further explore the contribution of these transcription factors in white tea processing, we analyzed the correlation (Pearson correlation coefficient, PCC) between the expression of the above key genes in the biosynthetic pathways of aroma substances and the expression of transcription factors (Fig. 5b & Supplemental Table S7). In the ERF gene family, there are five members (TEA003052, TEA013180, TEA026595, TEA013830, and TEA008270) that highly positively regulate the HPL gene member (TEA008699) and the ACX gene member (TEA016279) in the JA biosynthesis pathway (cor ≥ 0.8). ERFs mediating abiotic stress tolerance are induced by the JA signaling pathway^[41]. Interestingly, we also found a significant negative regulation of the LOX2 (TEA003727) by the bHLH gene family MYC2 (TEA010590) after 12 h of withering. It has been shown that MYC2 down-regulates its target gene LOX2 in Arabidopsis mutants subjected to external environmental influences^[42]. The expression of LOXs was decreased by the regulation of transcription factors in fresh leaves to maintain ROS homeostasis and prevent oxidative damage^[32], which in turn inhibits the biosynthesis of some ester aromas. We speculated that this result might be the reason for the weakening of the floral aroma of white tea in the late withering stage. Therefore, during the white tea withering process, the transcription factors were involved in the biosynthesis of aroma components of terpenoids and esters, and there was both positive and negative regulation.

Identification of aroma synthesis-related coexpression gene networks and key candidate genes

To further investigate aroma changes in tea leaves during withering, WGCNA was utilized to structure gene regulatory networks and to determine coexpressed gene sets. In total, nine coexpression modules (marked and highlighted with various colors), each containing between 34 and 2,770 genes, were identified by merged dynamic analysis of WGCNA. The previous analysis with aroma formation suggested higher levels of esters than terpenes and fatty acids and significant fluctuations throughout the withering process. Although esters have an important influence on the aroma and flavor of tea products, there are few studies on the molecular mechanisms of ester volatile compound formation during tea processing. Therefore, nine ester components, which were highly abundant in the volatile metabolome data and contributed to the enhancement of aroma quality, were selected for subsequent exploration.

Analysis of module-trait correlation with the nine ester volatiles as phenotypic traits demonstrated that the purple and magenta modules were positively related to the majority aroma accumulation (Fig. 5d), while the brown module was significantly negatively correlated. Thus, the genes in these three modules were used for further analysis. In Supplemental Fig. S5a & b, the different KEGG anabolic pathways were enriched in these three modules, and some overlap was observed. The KEGG species in the purple module were mainly clustered in plant hormone signal transduction, nitrogen metabolism and environmental information processing, while the magenta module was enriched in phenylpropanoid biosynthesis, signal transduction and biosynthesis of other secondary metabolites (Supplemental Table S6b). The KEGG enrichment in these two modules illustrates that the formation of ester aroma during white tea processing was also a result of the combined influence of the external environment, internal plant hormones and metabolic levels. Interestingly, the brown module, which was negatively correlated with ester aroma content, was enriched in KEGG categories such as transcription factors and metabolism of terpenoids and polyketides, suggesting that some genes may engage in the formation of esters and terpenes through transcriptional regulation patterns (Supplemental Table S6c).

A coexpression network was constructed from the top 200 topological overlap measure (TOM) values, in which genes with a high connectivity degree were screened as pivotal hub genes. Seven pivotal hub genes were detected in the brown module (Fig. 5e), including one serine/threonine-protein kinase gene (Nek2, TEA016476), one UDP-xylose transporter 3 gene (UXT3, TEA022766), one TRANVIA gene (TVA, TEA007818), one carboxylesterase 5 gene (CXE5, TEA029516), and three genes encoding proteins with unknown functions (TEA030989, TEA017052, and TEA026526). The CXE5 gene encodes a protein that functions in Arabidopsis thaliana to hydrolyze short-chain fatty acids^[43], and the expression of this gene is progressively downregulated with increasing tea leaf withering time, which might result in more volatile ester precursors eventually being retained. For the magenta and purple modules, a total of 11 genes were screened as pivotal hub genes using the same method as above. Notably, we found that the magenta module included an acyl-CoA thioesterase (ACOT, TEA020723) gene, and the purple module contained galactinol synthase 2 (GOLS2, TEA003328) and a chitinase (TEA019397). ACOT has the functional activity of hydrolyzing acyl-CoA to generate free fatty acids and coenzyme A^[44]. The total amount of fatty acids also increased (Fig. 2b). The accumulation of ACOT during white tea processing might also provide free fatty acids as part of the precursors for the formation of ester volatiles. GOLS2 has been demonstrated to be an important gene in response to moisture stress in poplar and rice^[45,46]. Thus, the increased expression of GOLS2 and the gradual loss of moisture during withering indicated that the postharvest tea leaves were under continuous drought stress. The gene TEA019397, which encodes an essential chitinase, engages in the ethylene/JA-mediated signaling pathway for acquiring systemic resistance^[47], and its expression level was also markedly increased during tea leaf withering. Therefore, combined with the results of transcription factor correlation analysis in Fig. 5b, we further speculated that when tea leaves experience drought stress due to withering, more free fatty acids would be released or retained to maintain leaf turgor pressure and act as signaling molecules (e.g., JA-mediated signaling pathway) to regulate various stress-responsive genes. Moreover, the coexpression network of the brown, purple and magenta modules suggests that these essential central genes are not only linked to numerous marginal genes but are also intensely interrelated with each other, suggesting that they cooperated or competed with each other in a gene network pattern in the accumulation of ester volatiles. In addition, we verified the accuracy of the RNA-seq data by randomly selecting 12 genes from the above modules for RT‒PCR analysis (Supplemental Fig. S6), and the results showed that the expression of these genes was relatively consistent under both assays.

In this research, the results of volatile metabolomics and transcriptomics analyses from postharvest leaves were integrated to reveal the role of withering on the aroma metabolism of white tea. Withering is a process in which the leaves of tea plants undergo drought stress. Unlike in oolong tea and black tea, withering in white tea takes a longer time and is the core process of character formation, which can significantly improve the aroma sensory evaluation. However, the degree to which different aroma components accumulated during the leaf withering process varied. Terpenoids, esters and fatty acids were the most abundant aroma volatiles accumulated in white tea during withering. By examining gene expression patterns in the pathways of α-linolenic acid metabolism and terpene biosynthesis, upregulated genes in the JA signaling pathway were found to be important for aroma formation and were induced by transcription factors. This observation is related to the mechanisms of response to abiotic stresses. Furthermore, we elucidated three coexpression module networks that were significantly associated with the changes in ester volatile components during withering, and revealed that the increase in free fatty acids might contribute to the accumulation of ester aroma precursors during the processing of white tea. Although the functions of the majority of pivotal hub genes in these networks remain unknown, some members have been demonstrated to be associated with the regulation of tea scent. Therefore, these genes can be considered critical candidates that deserve further investigation.