Melon ‘HT’ (Huai Tian No. 1, a climacteric cultivar of oriental melon, Cucumis melon var. Makuwa Makino) fruit were grown in a greenhouse on the Huaibei Normal University campus in Huaibei, China. Fruits were harvested at the ripe stage (30 days after pollination). Fruits of uniform size and without mechanical damage were selected, washed and dried at room temperature.

Fruit treatment was performed as described by Zhang et al. (2016). Briefly, melon fruits were divided into three groups: (A) stored at 4°C for 7 d (without chilling injury), and then transferred to room temperature (~22°C) for 1 d; (B) held at 4°C for 8 d until use (without chilling injury); (C) harvested on day 8 (8 d after the fruits of first two groups). Three replicates of each group were performed, and each replicate has three fruits.

A triangle test was performed according to the method described previously (Radovich et al., 2004). The test consisted of 58 untrained panelists, including twenty-eight males and thirty females, aged 20-45 years. Panelists were not trained but received instructions regarding the evaluation procedure. Panelists were presented with three samples of chilled (Group B) or unchilled (Group C) fruit, two were identical, and one was different. The samples were assigned blinding codes and randomized before testing. Panelists were instructed to identify the odd sample among the three by tasting all in the order presented. Panelists were provided with bottled water and sugar-free biscuits for palate cleansing, which they used between samples. Finally, they were asked to comment on the flavor and taste variation observed. The critical number of correct and incorrect responses for significance analysis was determined as described by Civille and Carr (2015). A binomial test was used to test for significance based on the correct number of accessors (Lawless and Heymann, 2010). The number of correct responses greater or equal to 27 of 60 panelists was considered a significant difference. Ambiguous responses were removed from the data. A plot of the frequency of the flavor descriptors for the chilled and unchilled fruit were obtained using the ggplot2 package in R (Wickham, 2016).

GC-IMS is a powerful technique for the separation and sensitive detection of VOCs and is widely used to analyze VOC changes during food storage (Wang et al., 2020). To investigate the volatile changes of HT1 melon cultivar during cold storage, the analysis of volatile fingerprints was performed on an HS-GC-IMS system (FlavourSpec^® Gesellschaft für Analytische Sensorsysteme mbH, Dortmund, Germany) as described by Yu et al. (2021). Three pooled fruit (three replicates for each group) were homogenized in a blender. For each replicate of melon sample, 10 g of melon juice was combined with 3 ml 200 mM ethylenediaminetetraacetic acid (EDTA) and 3 ml 20% CaCl₂ (m/v) in a 20 ml headspace vial, the samples were stored at -80°C and mixed well before use. The GC system is equipped with an automatic headspace sampler unit (CTC Analytics AG, Switzerland), which can collect samples from the headspace vials with a heated air-tight syringe (1 ml). Next, the prepared ‘HT1’ samples in a 20 ml headspace vial were incubated at 40°C for 20 min with an agitation speed of 450 rpm, and 500 μl of headspace was injected into the GC injector at 80°C using a heated gas-tight syringe (45°C) with an injection speed of 30 ml min^-1 and the GC program was as follow: 2 ml min^-1 for 2 min, ramped up to 100 ml min^-1 over 18 min, and held for 10 min. The separation was carried out on both a nonpolar FS-SE-54-CB capillary column (15 m × 0.53 mm × 0.25 μm, CS-Chromatographie Service GmbH, Durem, Germany) and a polar DB-Wax column (15 m × 0.53 mm × 0.25 μm, Agilent, Santa Clara, CA, USA) at isothermal condition of 60°C. Nitrogen was the gas at a flow rate of 150 ml min^-1. The drift tube is 5 cm in length, and it was operated at a constant voltage of 400 V cm^-1 at 45°C with a 150 ml min^-1 N₂ flow. Each IMS spectrum is an average of 12 scans. Each GC-IMS analysis was run in triplicate. To eliminate cross-contamination, the syringe was purged by an N₂ flow for 30 s before every analysis and another 5 min after every analysis. n-ketones C₃-C₉ (Sinopharm Chemical Reagent Beijing Co., Ltd., China) were employed as external references to calculate the retention index (RI) of each compound. All volatile compounds were identified by comparing RI and normalized drift time (RIP relative) with the authentic reference compounds in the GC-IMS library. Analysis was performed on three replicates of each group.

Total RNA from different treatments was extracted using TRIzol reagent (Life Technologies) and first-strand cDNAs were synthesized using the ReverTra Ace quantitative real-time PCR (qRT-PCR) Kit (Toyobo, Osaka, Japan) according to the manufacturer’s instructions. qRT-PCR was performed using the SuperReal PreMix Plus Kit (Tiangen, Beijing, China) in CFX Connect real-time PCR detection system (Bio-Rad, Hercules, CA, USA) according to the manufacturers’ instructions. Actin was used as the internal control (Primers used are listed in Supplementary Table S13 ). Total RNA samples were submitted for Next-Generation Sequencing (NGS) using Illumina NovaSeq (Shanghai Personal Biotechnology Co., Ltd, China). Briefly, the paired-end (PE) method was applied for RNA sequencing and the read length was 150 base pairs. The mapped reads for each sample ranged from 37,300,202 to 46,422,461 and the transcriptome coverage varied from 94.72% to 95.15%.

The procedure of analysis was performed as described by (Guo et al., 2021). Briefly, the adapter sequence and low-quality value (QV < 20) from the raw reads in FastQ format were removed using Cutadapt (Version 1.1) software (Kim et al., 2015). The pruned clean reads were aligned to the melon reference genomic sequence Melon (DHL92) v3.5.1 Genome database (http://cucurbitgenomics.org/organism/18) using HISAT2 software (v. 2.1.0). HTSeq (v. 0.1.1) was used to calculate the read count as the original gene expression. The transcript abundance was indicated by normalized Fragments Per Kilo bases per Million fragments (FPKM), and genes with an average value higher than 1 were considered as expressed genes, while the genes expressed at very low levels were removed. Pearson correlation coefficient was estimated to examine the correlation of gene expression levels among samples. Based on principal component analysis (PCA) in “DESeq” package in R software, samples in each treatment group were clustered according to relative expression. Absolute log2 fold change (FC) > 1 and adjusted p-value < 0.05 were differentially expressed genes (DEGs) in each treatment. Expression levels of DEGs in each group were visualized as a heat map using “Pheatmap” package. The DEGs in each treatment group were further submitted for functional enrichment analysis of gene ontology (GO) and transcription factor (TF). For GO enrichment analyses performed in “TopGO” package, DEGs were divided into GO molecular function (MF), biological process (BP), and cellular component (CC). The TF information and relative family classification were predicted according to Plant Transcription Factor Database (PlantTFDB).

Metabolite profiling was performed using a non-targeted metabolome method by Personal Biotechnology Co. Ltd. (Shanghai, China). The materials used for the transcriptomic analysis were the same as those in metabolite profiling. Extraction and analysis of metabolites were performed as previously described (Sangster et al., 2006; De Vos et al., 2007). 200 mg ( ± 1%) of melon fruit was accurately weighed in an EP tube, and then 0.6 ml 2-chlorophenylalanine (4 ppm) methanol (-20°C) was added, followed by vortexing, grinding, and centrifugation. Finally, 20 μl from each sample was taken to the quality control (QC) samples (These QC samples were used to monitor deviations of the analytical results from these pool mixtures and compare them to the errors caused by the analytical instrument itself), and the samples were used for LC-MS detection.

Chromatographic separation was accomplished in a Vanquish system (Thermo, USA) equipped with an ACQUITY UPLC^® HSS T3 (150 mm x 2.1 mm, Waters, USA) column maintained at 40°C. The temperature of the autosampler was 8°C. Gradient elution of analytes was carried out with 0.1% formic acid in water (A2) and 0.1% formic acid in acetonitrile (B2) or 5 mM ammonium formate in water (A3) and acetonitrile (B3) at a flow rate of 0.25 ml min^-1. Injection of 2 μl of each sample was done after equilibration. An increasing linear gradient of solvent B2/B3 (v/v) was used as follows: 0-1 min, 2% B2/B3; 1-9 min, 2%-50% B2/B3; 9-12 min, 50%-98% B2/B3; 12-13.5 min, 98% B2/B3; 13.5-14 min, 98%-2% B2/B3; 14-20 min, 2% B2-positive model (14-17 min, 2% B3-negative model). The ESI-MSn experiments were executed on the Thermo Q Exactive HF-X mass spectrometer with the spray voltage of 3.5 kV and -2.5 kV in positive and negative modes, respectively. Sheath gas and auxiliary gas were set at 30 and 10 arbitrary units, respectively. The capillary temperature was 325°C. The analyzer scanned over a mass range of m/z 81-1,000 for a full scan at a mass resolution of 60,000. Data dependent acquisition (DDA) MS/MS experiments were performed with HCD scan. The normalized collision energy was 30 eV. Dynamic exclusion was implemented to remove some unnecessary information in MS/MS spectra.

Co-expression networks were created using WGCNA (v1.29) package in R as described by (Langfelder and Horvath, 2008). The data used for WGCNA were performed on aroma components (alcohols, aldehydes, ketones, and esters) produced during storage. All genes and TFs were used as input to the signed WGCNA network construction. Modules were identified using “one-step network construction and module detection function” method with a soft thresholding power of 6 and a relatively large minimum module size of 30, and the threshold for merging modules was 0.25.

Promoter sequences (base pairs -1,500 to -1) were obtained from Melon (DHL92) v3.5.1 Genome (http://cucurbitgenomics.org/organism/3). Transcription binding motifs were analyzed with plantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) and PlantPAN 3.0 (http://plantpan.itps.ncku.edu.tw/promoter.php) (Lescot et al., 2002; Li et al., 2020).

Student t-test analysis was applied to compare the variation of volatile content and gene expression between Group A & B, B & C and A & C (Graphpad Prism 9). Letters and asterisks indicate significant differences (p-value < 0.05) for treatments at each measure time. PCA, DEGs, and GO enrich were described in each method section.

All accession numbers were obtained from the Melon (DHL92) v3.5.1 Genome database. CmLOX18 (MELO3C024348.2), CmADH1/2 (MELO3C023685.2/MELO3C014897.2), CmHPL1/2 (MELO3C010910.2/MELO3C018413.2), CmAAT1/3 (MELO3C024771.2/MELO3C024762.2), CmMGL (MELO3C013774.2), CmPDC1/2 (MELO3C009145.2/MELO3C007227.2), CmCXE1/2 (MELO3C010887.2/MELO3C014690.2), CmBCAT1 (MELO3C010776.2), CmArAT1 (MELO3C025613.2), CmAADC1 (MELO3C008357.2), CmACO1 (MELO3C014437.2), CmACS1 (MELO3C021182.2), CmNOR (MELO3C016540.2), CmMYB44 (MELO3C007586.2), CmbZIP1 (MELO3C013408.2), CmCBF1 (MELO3C006869.2).

To gain insight into how chilling affects flavor of melon fruit and consumer liking, we performed a triangle test to evaluate the flavor differences between chilled and unchilled fruit. Twenty-seven of fifty-eight panelists were able to identify significant differences in the flavor of the chilled compared to the unchilled (n > 26, p = 0.02 < 0.05) ( Figure 1A ). In addition, they described the unchilled melon as sweeter, juicier, and more aromatic than the chilled fruit, indicating that chilling had an adverse effect on flavor ( Figure 1B ).

Subsequently, VOCs extracted from melon fruit after 7 d chilling followed by 1-d recovery at room temperature (Group A), after 8 d chilling only (Group B), and harvested on day 8 (Group C) were quantified. Sixty known and nine unknown VOCs were detected in fruit from all three groups and used for further analysis ( Supplementary Table S1 ). These VOCs consist of nine alcohols, twelve aldehydes, six ketones, twenty-four esters, and eighteen others ( Figure 2A ). Among the VOCs, esters are the most abundant (> 70% of total VOC content) in all groups, and ~70% of the esters are VAEs. Based on the VOC content, melon fruit from the three different groups were separated by PCA, and biological replicates clustered together ( Supplementary Figure S1 ). 2D ion mobility spectrums of volatiles in fruit from these three groups were shown ( Supplementary Figure S2 ). In general, the concentration of volatiles in 8-d chilled fruit (Group B) was lower in comparison with unchilled fruit (Group C). Total amount of volatiles in 8-d chilled fruit decreased by ~73% relative to the unchilled. After 1-d recovery at room temperature (Group A), total VOC content increased, but was still lower compared to the unchilled ( Figure 2B ). In particular, content of VAEs was lower after chilling and recovery at room temperature for 1 d. However, total content of volatile alcohols or aldehydes was not substantially affected after cold storage compared to the unchilled ( Figures 2B, C ). Additionally, no significant differences were observed in fructose, glucose and citrate between the chilled and the unchilled fruit ( Supplementary Figure S3 ).

A gallery plot was performed to better visualize major chilling-induced changes to VOC profiles. A total of eighty-two detected signal peaks were characterized ( Figure 2D ). Although total content of alcohols or aldehydes was not altered after cold storage compared to the unchilled fruit, content of several alcohols and aldehydes was lower, including 3-methylpentanol, nonanal and (E, E)-2,4-nonadienal. In addition, total ketone content was also reduced in 8-d chilled fruit. Among the ketones, content of 3-hepten-2-one and 2-propanone was lower. Total content of VAEs in the chilled fruit was ~30% lower compared to that in unchilled fruit and those recovered after transfer to the room temperature. Content of multiple VAEs was reduced after cold storage, including propyl acetate, isobutyl acetate, pentyl acetate, hexyl acetate, and benzyl acetate. Chilling-induced decreases were also observed for other esters, including isopentyl isopentanoate, butyl butanoate, isobutyl butanoate, and isopentyl hexanoate. Content of the aforementioned VAEs increased after 1-d recovery. However, some VOCs were found to increase after cold storage compared to the unchilled. Interestingly, these VOCs mainly included methyl esters, such as methyl butanoate, methyl hexanoate, methyl 3-methylbutanoate, and methyl 2-methylpropanoate ( Supplementary Table S1 ). In addition, content of these methyl esters was reduced after return to room temperature. These results indicate that chilling affects the flavor of melon fruit, and chilling-induced changes to melon VOC profiles are mainly involved in flavor-associated esters, particularly VAEs.

The aforementioned VAEs are derived from fatty acid and amino acid pathways. In response to chilling, metabolic flux throughout these two pathways is critical for downstream VAE biosynthesis. Therefore, multiple intermediate metabolites were investigated, including linoleic acid, linolenic acid, 9(S)-hydroperoxy-10(E),12(Z),15(Z)-octadecatrienoic acid (9(S)-HPOT), and 13(S)-hydroperoxy-9(Z),11(E),15(Z)-octadecatrienoic acid (13(S)-HPOT) from the fatty acid pathway, in addition to isoleucine, L-leucine, valine, phenylalanine, 3-Methyl-2-oxovaleric acid and Ketoisocaprorate (KIC) from the amino acid pathway. According to the analysis of differentially accumulated metabolites, we found that content of linoleic acid, linolenic acid, and 13S-HPOT increased after chilling and was reduced after 1-d recovery at room temperature ( Figure 2E ; Supplementary Table S2 ). However, content of metabolites from the amino acid pathway showed inconsistent patterns in response to chilling. Content of L-leucine was higher in the chilled fruit than the unchilled and remained high after cold storage followed by 1-d recovery. Content of KIC, a ketoacid generated from L-leucine, was reduced after chilling and slightly increased after recovery. Phenylalanine content was lower after cold storage and remained at a low level after 1 d returning to room temperature ( Figure 2E ; Supplementary Table S2 ). These results suggested that chilling affects VAE biosynthesis via regulation of metabolic flux.

RNA-Seq data was generated for transcriptomic analysis. After discarding adaptor and low-quality reads, ~6.2 Gb clean reads (>90% of total reads) on average were produced ( Supplementary Table S3 ). Principle components analysis (PCA) of RNA-Seq data revealed the relationship between experimental samples, which exhibited clear separation for the samples of the aforementioned three groups ( Figure 3A ). Pearson correlation coefficient analysis indicated that, in general, replicates after chilling followed by 1-d recovery at room temperature (Group A) and harvested on day 8 (Group C) were less correlated with the replicates after chilling (Group B) ( Figure 3B ).

Analysis of Group A & B, and Group B & C identified 8553 (4272 upregulated and 4281 downregulated) and 8249 (4272 upregulated and 3977 downregulated) DEGs, respectively. However, only 1824 upregulated and 653 downregulated DEGs were identified via analysis of Group A and C ( Figures 3C, D ). Based on the expression levels of genes from each group, correlation among the samples was calculated and hierarchical clustering was visualized in a heatmap. Consistent with PCA analysis, gene expression clusters in Group B were distinguishable from that in Group A and Group C ( Supplementary Figure S4 ). These results indicate that expression levels of many genes are sensitive to temperature shift. The C-repeat binding factor (CBF)/dehydration responsive-element binding (DREB) TFs are known to be sensitive to low temperature (Thomashow, 2010; Zhang et al., 2016). In the present study, we found that transcript levels of multiple melon CBFs/DREBs were higher in the chilled melon and returned to unchilled levels after transfer to room temperature ( Supplementary Figure S5 ).

To understand the molecular basis of the reduction of melon flavor-associated VAEs during cold storage, we investigated the relevant biosynthetic pathways and genes involved in the biosynthesis of VAEs. As observed in the aromatic melon fruit that we used, content of volatile esters was much higher than other characterized VOCs, such as alcohols and aldehydes ( Figure 2B ; Supplementary Table S1 ). Alcohol dehydrogenase (ADH) and alcohol acyltransferase (AAT) are able to convert volatile aldehydes to their corresponding alcohols and esters ( Figure 4A ). The expression levels of two ADH family members (CmADH1 and CmADH2) were downregulated after 8-d cold storage and upregulated to higher levels after 1-d recovery at room temperature, and the transcript levels after recovery were higher than day 8. Similarly, the transcript abundance of two members of AAT family (CmAAT1 and CmAAT3) was lower in the chilled melon and recovered after 1 d at room temperature ( Figure 4B ; Supplementary Tables S4 , S5 ).

In the fatty acid pathway, biosynthesis of C5 and C6 green leafy volatiles depends upon LOX and HPL activities ( Figure 4A ). Here, we observed CmLOX18 expression level was lower after chilling and recovered after 1 d at room temperature. In addition, expression levels of CmLOX8 and a 9-lipoxygenase (CmLOX9) were also altered after chilling. Transcript levels of two CmHPLs were lower in the chilled melon. CmHPL2 expression was recovered after return to the room temperature, while CmHPL1 remained at low level. CmPDC1 involved in aldehyde biosynthesis through catalyzing straight-chain ketoacids was reduced after chilling and recovered to higher level after 1 d at room temperature ( Figure 4B ; Supplementary Tables S4 , S5 ). In the amino acid pathway, CmBCAT1 and CmArAT1 expression levels were reduced in the chilled melon fruit and upregulated after 1 d of recovery. Moreover, multiple aldehyde dehydrogenases (ALDHs) were also sensitive to low temperatures ( Figure 4B ; Supplementary Tables S4 , S5 ). Therefore, reduced production of VAEs after chilling cannot be directly explained by a single enzyme activity alone.

According to the expression patterns of DEGs among the three groups, clustering analysis was carried out, and DEGs were divided into nine clusters ( Supplementary Figure S6 ). Based on the Plant/Animal Transcription Factor Database (Plant/AnimalTFDB), thirty TF families were identified ( Supplementary Figure S7 ; Supplementary Table S6 ). To narrow the scope of differentially expressed TFs, we removed the genes without significance (p-value > 0.01), and 375 TFs were selected. 246 TFs belong to Cluster 1 and 5, indicating that their expression levels were higher after chilling and recovered after 1 d at room temperature. On the contrary, 129 TFs were clustered in Cluster 6 and 8, indicating that their expression levels were lower after chilling and increased after 1-d recovery ( Supplementary Table S7 ).

Clusters 6 and 8 identified melon NONRIPENING (CmNOR) which is a key ripening regulator, and its expression was reduced (~30-fold reduction) after chilling and upregulated (~40-fold increase) after 1-d recovery at room temperature, which is consistent with the pattern of ester levels and flavor intensity ( Supplementary Tables S4 , S5 ). In the same cluster with CmNOR, CmMYB44, an ester-associated apple MdMYB6 homolog, had lower expression during cold storage and increased transcript after 1-d recovery ( Supplementary Tables S4 , S5 ). In addition, Cluster 1 and 5 identified a bZIP family gene (named CmbZIP1, ester-associated MabZIP4 homolog in banana) with the highest expression during chilling and reduced transcript levels after recovery ( Supplementary Tables S4 , S5 ).

Ethylene also regulates various aspects of ripening, including VOC biosynthesis in fruit. Ethylene biosynthesis genes, including aminocyclopropane carboxylate synthase (ACS) and oxidase (ACO), were investigated. The transcript levels of CmACS1 and CmACO1 were the highest among the family members, respectively. The expression of CmACO1 was lower in the chilled melon and increased after 1 d at room temperature. Conversely, much higher transcript levels of CmACS1 was observed after cold storage and declined after returning to room temperature ( Supplementary Figure S8 ; Supplementary Tables S4 , S5 ).

To gain further insight into the regulation of ester biosynthesis after temperature shift, WGCNA was performed to investigate the co-expression networks of key DEGs. According to the similarity of expression patterns, genes were grouped into a total of eighteen co-expression modules ( Figure 5A ). Based on the top three coefficient values shown in the heatmap of module-trait correlations, content of esters was positively correlated with gene expressions in the blue, turquoise and green-yellow modules ( Figure 5B ; Supplementary Table S8 ). Subsequently, GO annotation analysis was performed to summarize functions of the DEGs identified in the modules, and DEGs in response to chilling were investigated. The functional categories of the DEGs included biological process, molecular function and cellular component. Top twenty GO terms indicated that cellular process in the biological process ranked the highest, and 753 DEGs were divided into the group of cellular aromatic compound metabolic process ( Figure 6 ; Supplementary Table S9 ).

Based on the analyses of transcriptome and metabolome, multiple genes involved in ester biosynthesis were preliminarily identified. Notably, WGCNA further narrowed the scope of DEGs correlated with ester biosynthesis. Two lipoxygenases (CmLOX8 and CmLOX18), one hydroperoxide lyase (CmHPL2), one alcohol acyltransferase (CmAAT1), one carboxylesterase (CmCXE1), and one amino acid decarboxylase (CmAADC1) were identified. These DEGs were confirmed by qRT-PCR using RNA from fruit after chilling followed by 1-d recovery at room temperature (Group A), chilled (Group B) and unchilled fruit (Group C) ( Figure 7 ). Additionally, WGCNA was also performed to investigate the correlation between differentially expressed TFs and ester content during chilling and 1-d recovery after chilling. Seventeen key TFs (GS > 0.85, GS < -0.85) were identified, including twelve positive regulators (one C3H, one CPP, one ERF, two GATAs, one GeBP, one LBD, one MYB, one STAT and three WRKYs) and five negative regulators (one ARR-B, two bHLHs, one DBB and one MYB-related) ( Table 1 ). These positive regulators exhibited lower expression levels in the chilled fruit and were upregulated after 1-d recovery at room temperature, which was consistent with the pattern of altered ester content. However, the negative regulators exhibited higher levels after chilling and were downregulated after 1-d recovery. Although the aforementioned three TFs, CmNOR, CmMYB44, and CmbZIP1 were not included in the list of top TFs, WGCNA revealed that the relative GS values were higher than 0.7 or lower than -0.7, indicating significant correlations between expression and ester content ( Supplementary Table S10 ).

Under chilling conditions, volatile aldehydes and alcohols were not substantially affected, and the last step of ester biosynthesis catalyzed by AAT seems to be an important control point. cis-regulatory elements in promoter regions containing TF binding sites play a critical role in regulation of gene expressions. Thus, 1.5-kb putative promoter sequences of CmAAT1 and CmAAT3 were investigated.

The analysis of CmAAT1 and CmAAT3 promoter regions identified motifs mainly involved in response to stress and hormone levels. For CmAAT1, stress responses included anaerobic induction (ARE), stress (STRE), wounding (WUN-motif), and hormone response consists of ethylene (ERE), MeJA (MYC and TGACG-motif) and salicylic acid (TCA-element). For CmAAT3, anaerobic induction (ARE), light (AE-box), wounding (WRE3) and ethylene (ERE), gibberellin (TATC-box), MeJA (MYC), salicylic acid (TCA-element) were found to be involved in stress and hormone responsiveness, respectively. In addition, various conserved motifs were found in the promoter region of CmAAT1, including G-box (bZIP binding site), MYB and ‘ACCGAC’ (DREB/CBF) and ‘TTA/GCGT’ (NAC). MBSI (MYB) and ‘TTA/GCGT’ (NAC) were found in the CmAAT3 promoter region.

After temperature shift, we found that expression levels of many TFs exhibited the same pattern with the CmAAT1 and CmAAT3. These TFs included Dof, ERF, WRKY and TCP, in addition to the aforementioned bZIP, MYB, MYC and NAC ( Table 1 ; Supplementary Tables S4 , S5 ). The putative binding sites of some selected TFs were also found in the 1.5-kb promoter sequence of CmAAT1 and CmAAT3 ( Figure 8 ), and all putative motifs were shown in Supplementary Table S11 and S12 . These observations indicated potential roles of the aforementioned TFs in regulation of VAE-related AAT expression during cold storage.