An Integrated Analysis of Metabolomics and Transcriptomics Reveals Significant Differences in Floral Scents and Related Gene Expression between Two Varieties of Dendrobium loddigesii

Flower fragrance is one of the traits that holds important economical values in flowering plants. Extensive attention has converged on fragrance preservation in flower cultivation and breeding. Dendrobium loddigesii is an important species for cultivating aromatic Dendrobium orchid varieties for the long term due to its fragrance. Few studies focus on exploring related genes responsible for the aroma components in D. loddigesii. We analyzed flowers from two aromatic D. loddigesii varieties using high-throughput RNA sequencing and gas chromatography-mass spectrometry (GC-MS). The metabolomics results showed that the main volatile compounds responsible for the aroma formation of D. loddigesii were terpenes, especially monoterpenes. The de novo transcriptome assembly comprised 175,089 unigenes, and 24,570 unigenes of the genes were identified as differential expressed genes (DEGs) between the two varieties. Among these DEGs, 525 genes were mapped into seven pathways that related to the floral scent synthesis. Seventeen key genes were significantly correlated with volatile aroma metabolites, including geraniol, α-pinene, eugenol, and (Z)-3-hexenal. These results provide references for understanding the aroma biosynthesis and perfume formulations of D. loddigesii.


Introduction
Dendrobium loddigesii Rolfe, a perennial epiphytic herb in the Orchidaceae family [1], is a precious medicinal material in traditional Chinese medicine because of its secondary metabolites production [2,3]. Furthermore, D. loddigesii has a high economic value in ornamental horticulture due to its pleasing morphology and fragrance, and it is often used as the source material for essence and fragrance extraction in the cosmetic industry [4]. Floral scent is one of the most significant traits that constitute the ornamental value of flowering plants [5,6]. However, due to the complexities in the composition, structure, and biosynthesis process of floral scent, research findings are insufficient compared to other important ornamental characteristics such as flower shape and color [7,8]. Meanwhile, floral scent is not a breeding criteria during Dendrobium hybridization for the long term, which leads to the lack of fragrance in many commodity dendrobiums [9]. The selection of aromatic varieties has become one of the most urgent needs, and important trends in the breeding of Dendrobium, in which the role of D. loddigesii cannot be ignored.
Floral scent is composed of various low-molecular-weight volatile chemicals mostly belonging to three major groups: terpenoids, benzenoids/phenylpropanoids, and fatty 2 of 13 acid derivatives [10]. The floral scents of common ornamental plants such as rose and peony have been well studied [11,12], some research studies about aromatic compounds of Dendrobium also have been reported. In total, 33 volatile organic compounds (VOCs) were identified in D. chrysotoxum, including 15 terpenes, 7 esters, 6 alcohols, 2 aromatic compounds, 2 ketones, and 1 aldehydes [9]. Meanwhile, 41 VOCs were distinguished at the flowering stage of D. moniliforme, and the alkenes were the most important compounds, which affected its floral scent [13]. A recent study has identified esters as the main components of the aroma in D. loddigesii sampled from Yunnan province, China [14]. However, another research suggests that the main ingredients of floral scent in D. loddigesii are terpenes [15]. These findings suggested the complexity of volatile components in flowers of Dendrobium species.
Due to the invisibility and highly variable characters of floral scent and lack of efficient methods to identify genetic variations, research on Dendrobium has focused on the stage of volatiles identification for a long time. The genetic and molecular mechanisms responsible for the floral scent in Dendrobium species are poorly revealed. Thanks to the rapid developments of molecular biology and high-throughput sequencing technologies, recent research studies on floral scent have elucidated several biosynthetic pathways of fragrance compounds and some structural genes encoding key enzymes in D. officinal [16,17]. Nonetheless, there are few studies on the molecular mechanism of synthesis of D. loddigesii's aromatic ingredients.
To investigate the characteristic and molecular mechanisms of floral scent in D. loddigesii, the flowers of high-scented D. loddigesii (TH) originated from Guizhou province and light-scented D. loddigesii (TL) from Guangdong province of China were used as materials to screen differential metabolites (DMs). There were obvious differences in olfactory sense between the two varieties: flowers of TH had a strong floral smell, while TL has a grassy smell. Furthermore, we conducted the first de novo assembly of the transcriptome for the species by high-throughput RNA sequencing. The differential expressed genes (DEGs) and enriched pathways between two varieties were identified. Our findings provide important gene resources for further study of the biosynthesis pathway and regulation mechanism responsible for floral volatiles in D. loddigesii. This work holds merits for exploring the molecular mechanism of aroma biosynthesis of D. loddigesii and providing preliminary preparations for breeding new aromatic Dendrobium varieties and perfume blends.

Plant Material and Sample Collection
Two varieties of D. loddigesii, TH and TL, were cultivated in the germplasm garden of the Institute of Horticulture, Guizhou Academy of Agricultural Sciences (Guiyang, China). Flowers of two varieties were collected in the blooming stage ( Figure 1). Nine replicates for each of the two varieties of flowers were sampled in total, including six replicates for volatile compounds and three replicates for transcriptomics analysis, respectively.
To investigate the expression patterns of screened genes, a series of samples from floral organs were collected for RNA extraction, including blooming flowers of two varieties and flowers of TH in different flowering stages ( Figure 1). All samples were immediately frozen in liquid nitrogen and stored at −80 • C until needed.

Gas Chromatography-Mass Spectrometry Analysis of Volatile Compounds in Flowers of Two Varieties
Headspace SPME was employed to collect the volatile compounds from flower tissues, which were absorbed by a 65 µm DVB/CAR/PDMS fiber for 1 h at 20 • C. After sampling, desorption of the VOCs from the fiber coating was carried out in the injection port of the GC apparatus (Model 7890B, Agilent, Santa Clara, CA, USA) at 250 • C for 5 min in the splitless mode.

Gas Chromatography-Mass Spectrometry Analysis of Volatile Compounds in Flowers Varieties
Headspace SPME was employed to collect the volatile compounds from flow sues, which were absorbed by a 65 µm DVB/CAR/PDMS fiber for 1 h at 20 °C. Aft pling, desorption of the VOCs from the fiber coating was carried out in the injecti of the GC apparatus (Model 7890B, Agilent, Santa Clara, CA, USA) at 250 °C for 5 the splitless mode.
The identification and quantification of VOCs was carried out using an Model 7890B GC and a 7000D mass spectrometer (Agilent) equipped with a 30 m mm × 1.0 µm DB-5MS (5% phenyl-polymethylsiloxane) capillary column. Heliu used as the carrier gas at a linear velocity of 1.0 mL/min. The injector temperatu kept at 250 °C and the detector was kept at 280 °C. The oven temperature w grammed from 40 °C (5 min), increasing at 6 °C/min to 280 °C, and then holding fo Mass spectra were recorded in electron impact (EI) ionization mode at 70 eV. Th rupole mass detector, ion source, and transfer line temperatures were set, respecti 150, 230 and 280 °C. Mass spectra were scanned in the range m/z 30-350 amu at 1 vals.
The identification of volatile compounds was achieved by comparing the ma tra with the database (NIST and Wiley). The relative contents of volatile compound calculated by the peak area normalization method (ratio of each peak area to th area). Differential metabolites (DMs) were defined as those fold changes (FC) ≥ 2 among them, with a variable importance for projection (VIP) value ≥ 1 and p-valu between two groups of samples.

RNA Sequencing
Total RNA was extracted from samples using Trizol (Invitrogen, MA, USA) ing to the manufacturer's protocol. The quality and quantity of RNA were checke a Qubit ® 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA) and analyzed in a 1 rose gel electrophoresis. Sequencing libraries were constructed using NEBNext ® U RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) following the m The identification and quantification of VOCs was carried out using an Agilent Model 7890B GC and a 7000D mass spectrometer (Agilent) equipped with a 30 m × 0.25 mm × 1.0 µm DB-5MS (5% phenyl-polymethylsiloxane) capillary column. Helium was used as the carrier gas at a linear velocity of 1.0 mL/min. The injector temperature was kept at 250 • C and the detector was kept at 280 • C. The oven temperature was programmed from 40 • C (5 min), increasing at 6 • C/min to 280 • C, and then holding for 5 min. Mass spectra were recorded in electron impact (EI) ionization mode at 70 eV. The quadrupole mass detector, ion source, and transfer line temperatures were set, respectively, at 150, 230 and 280 • C. Mass spectra were scanned in the range m/z 30-350 amu at 1 s intervals.
The identification of volatile compounds was achieved by comparing the mass spectra with the database (NIST and Wiley). The relative contents of volatile compounds were calculated by the peak area normalization method (ratio of each peak area to the total area). Differential metabolites (DMs) were defined as those fold changes (FC) ≥ 2 or ≤0.5, among them, with a variable importance for projection (VIP) value ≥ 1 and p-value < 0.05 between two groups of samples.

RNA Sequencing
Total RNA was extracted from samples using Trizol (Invitrogen, MA, USA) according to the manufacturer's protocol. The quality and quantity of RNA were checked using a Qubit ® 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA) and analyzed in a 1% agarose gel electrophoresis. Sequencing libraries were constructed using NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) following the manufacturer's recommendations. An AMPure XP system (Beckman Coulter, Beverly, MA, USA) was used to screen 250-300 bp cDNA, and polymerase chain reaction (PCR) amplification was performed. RNA-Seq was performed by MetWare (Wuhan, China) on the Illumina Novaseq platform.
Raw reads obtained from RNA-Seq were preprocessed using fastp, and adapters were trimmed. Low-quality sequences with N base content exceeding 10% or base (Q ≤ 20) content exceeding 50% were removed. All trimmed clean reads were assembled using Trinity, which generated contigs, transcripts, and unigenes [18]. All assembled unigenes were compared against databases including the NR (NCBI nonredundant), Gene Ontology (GO), SwissPro, euKaryotic Ortholog Groups (KOG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) for functional annotation.

Analysis of Differential Expressed Genes (DEGs)
The expression level of unigenes was calculated using fragments per kb of transcript per million (FPKM) method [19]. The p-value was corrected using the Benjamini and Hochberg method to get FDR (False Discovery Rate), a threshold of FDR < 0.05, and an absolute value of log 2 (fold change) > 1 were used to judge the significance of the DEGs between two varieties of flowers through DESeq2 [20]. ClusterProfiler R software was used to analyze the statistical enrichment of DEGs in the KEGG pathway and GO terms [21].

Quantitative Real-Time PCR (qRT-PCR)
Total RNA from different samples was isolated using the RNAprep pure Tissue Kit (Tiangen, Beijing, China). First-strand cDNA synthesis was performed using the Primescript RT reagent kit with a gDNA eraser (TaKaRa, Tokyo, Japan) according to the manufacturer's instructions. qRT-PCR was accomplished using the TB Green ® Premix Ex Taq™ (TliRNaseH Plus, ROX plus) (TaKaRa, Tokyo, Japan) on an ABI 7500 system (Applied Biosystems, Foster City, CA, USA).
The UBQ (ubiquitin) gene was used as an internal control, and each analysis had three biological replicates. The relative expression levels were calculated using the 2−∆∆Ct method [22], and the primers sequences are shown in Table S1.

Combined Analysis of DEGs and DMs
In order to gain an improved insight into links between genes and metabolites, DEGs and DMs were simultaneously mapped to the KEGG pathway. Pearson correlation coefficient analysis was used to calculate the correlation between expression of DEGs and relative contents of DMs in the same pathway.

Analysis of Metabolites between Two Aromatic D. loddigesii
The metabolomics analysis showed that 80 volatile components were detected by GC-MS from two different aromatic varieties of D. loddigesii ( Figure 2). The components with the highest number of species were aldehydes (15), followed by terpenes (14) and alkanes (11). Of the 80 metabolites, 75 metabolites were shared by both varieties. Four metabolites consisted of three terpenes (geraniol, p-menthene and o-cymene) and one aromatics (1-ethyl-4-methoxy-benzene) that were specific to the high-scented variety (TH), while only one phenol metabolite (4-propyl-phenol) was specific to the light-scented variety (TL).

Illumina Sequencing and De Novo Transcriptome Assembly
In order to explore the DEGs and related pathways that are involved in the difference of flower fragrance between two varieties of aromatic D. loddigesii, two groups of samples from high-scented and light-scented D. loddigesii were collected and submitted for RNA-seq using the Illumina paired-end sequencing technique.
In total, 271,241,852 sequencing data of raw reads were obtained. After quality control and data filtering, 127,791,342 and 130,146,106 high-quality clean reads were acquired, and the Q30 percentage of clean bases were 93.36 and 93.23%, respectively ( Table 1). All of these clean reads were employed for further de novo assembly; 230,951 transcripts and 175,089 unigenes were generated using Trinity software. The mean length of transcripts in the de novo transcriptome assembly was 952 bp with an N50 of 1627 bp. Furthermore, the mean length of unigenes was 1168 bp with a N50 of 1745 bp (Table 2).

Functional Annotation and Classification of Unigenes
In order to provide putative annotations for the assembled unigenes, the 175,089 unigenes were aligned against public databases including NR, SwissPro, KEGG, KOG, and GO databases using BLASTX or BLASTP (E-value ≤ 1.0 × 10 −5 ). A total of 121,096 unigenes (69.16% of the total unigenes) were matched to at least one of these databases. Then, 68.88% unigenes were matched to the proteins in the NR database, followed by 40.14% in the SwissProt database, 49.43% in the KEGG database, 36.65% in the KOG database, and 40.72% in the GO database (Table 3). Species homology distribution of D. loddigesii transcripts were determined based on NR annotation. The highest number of unigenes matched was D. catenatum (107,679 unigenes; 89.29%), Phalaenopsis equestris (5328; 4.42%), and Apostasia shenzhenica (778; 0.65%). The functional classification of the D. loddigesii transcriptome was performed by analyzing the BLAST results against the GO and COG databases. A total of 71,288 unigenes were assigned to 58 functional GO subcategories consisting of 28 in 'biological process', 12 in 'molecular function', and 18 in 'cellular constituent'. The most representative terms in the biological processes group were 'cellular process' (48,302; 67.76%) and 'metabolic process' (42,720; 59.93%); the top two terms in the cellular components groups were 'cell' (57,332; 80.42%) and 'cell part' (57,246; 80.3.0%), while the 'binding' (44,898; 62.98%) and 'catalytic activity' (42,720; 59.93%) were the dominant terms in the molecular functions groups ( Figure S1).

Functional Assessment and Verification of DEGs
The FPKM method was employed to compute the unigenes expression level. DEGs were screened out by comparison between the two samples, the parameters (FDR < 0.05 and |log 2 (fold change)| ≥ 1) were used as the threshold to judge the significance of the DEGs, which resulted in a set of 24,570 unique DEGs. Among them, 15,889 were up-regulated and 8681 were down-regulated in the high-scent (TH) relative to light-scent flowers (TL) (Figure 4). To determine the functions associated with the DEGs changes between two groups of samples, GO and KEGG pathway analysis were performed. groups ( Figure S1).

Functional Assessment and Verification of DEGs
The FPKM method was employed to compute the unigenes expression level. DEGs were screened out by comparison between the two samples, the parameters (FDR < 0.05 and |log2 (fold change)| ≥ 1) were used as the threshold to judge the significance of the DEGs, which resulted in a set of 24,570 unique DEGs. Among them, 15,889 were up-regulated and 8681 were down-regulated in the high-scent (TH) relative to light-scent flowers (TL) (Figure 4). To determine the functions associated with the DEGs changes between two groups of samples, GO and KEGG pathway analysis were performed. In total, 24,570 DEGs were classified into 57 GO terms; there are 28, 18, and 11 terms in the three major categories, respectively. Among them, the terms 'cellular process' and 'cell' are dominant in the biological process and cellular component categories, respectively, while 'binding' is dominant in the category of molecular function ( Figure S3). To further understand the function of DEGs in biological processes, 140 KEGG pathways showed hits with DEGs, of which 36 pathways were significantly enriched with DEGs. The top-20 enriched pathways were visualized on a dot bubble. The 'metabolic pathways' (ko01100) and 'biosynthesis of secondary metabolites' (ko01110) pathway were the top 2 remarkably enriched KEGG pathways ( Figure 5). In total, 24,570 DEGs were classified into 57 GO terms; there are 28, 18, and 11 terms in the three major categories, respectively. Among them, the terms 'cellular process' and 'cell' are dominant in the biological process and cellular component categories, respectively, while 'binding' is dominant in the category of molecular function ( Figure S3). To further understand the function of DEGs in biological processes, 140 KEGG pathways showed hits with DEGs, of which 36 pathways were significantly enriched with DEGs. The top-20 enriched pathways were visualized on a dot bubble. The 'metabolic pathways' (ko01100) and 'biosynthesis of secondary metabolites' (ko01110) pathway were the top 2 remarkably enriched KEGG pathways ( Figure 5).
The main components of floral scent are terpenoids, phenylpropanoids/benzenoids, and volatile fatty acid derivatives [10]. Based on these divisions, a total of 525 DEGs were annotated in seven pathways related to the floral scent synthesis (Table S3). Twelve key DEGs were identified and selected to verify the reliability of RNA-Seq data and explore the expression patterns of genes related to synthesis of floral scent using qRT-PCR analysis( Figure 6). The result revealed that in the terpenoids-related pathways, ECHS (limonene and pinene degradation, ko00903) and GPPS (terpenoid backbone biosynthesis, ko00900) had higher expression and FLDH (ko00900) displayed lower expression in the TH. In the phenylpropanoid biosynthesis pathway (ko00940), CCR and POX were up-regulated in TH, while the expression level of PAL, C4H, CAD, and TOGT were enhanced in TL. In the alpha-linolenic acid metabolism (ko00592), ACAA was up-regulated in TH, and TL displayed a higher mRNA level of AOC and LOX. Overall, qRT-PCR results indicated that The main components of floral scent are terpenoids, phenylpropanoids/benzenoids, and volatile fatty acid derivatives [10]. Based on these divisions, a total of 525 DEGs were annotated in seven pathways related to the floral scent synthesis (Table S3). Twelve key DEGs were identified and selected to verify the reliability of RNA-Seq data and explore the expression patterns of genes related to synthesis of floral scent using qRT-PCR analysis( Figure 6). The result revealed that in the terpenoids-related pathways, ECHS (limonene and pinene degradation, ko00903) and GPPS (terpenoid backbone biosynthesis, ko00900) had higher expression and FLDH (ko00900) displayed lower expression in the TH. In the phenylpropanoid biosynthesis pathway (ko00940), CCR and POX were up-regulated in TH, while the expression level of PAL, C4H, CAD, and TOGT were enhanced in TL. In the alpha-linolenic acid metabolism (ko00592), ACAA was up-regulated in TH, and TL displayed a higher mRNA level of AOC and LOX. Overall, qRT-PCR results indicated that the genes performed the same expression trends as the RNA-seq data, suggesting the high level of reliability of the RNA-seq data.   The main components of floral scent are terpenoids, phenylpropanoids/benzenoids, and volatile fatty acid derivatives [10]. Based on these divisions, a total of 525 DEGs were annotated in seven pathways related to the floral scent synthesis (Table S3). Twelve key DEGs were identified and selected to verify the reliability of RNA-Seq data and explore the expression patterns of genes related to synthesis of floral scent using qRT-PCR analysis( Figure 6). The result revealed that in the terpenoids-related pathways, ECHS (limonene and pinene degradation, ko00903) and GPPS (terpenoid backbone biosynthesis, ko00900) had higher expression and FLDH (ko00900) displayed lower expression in the TH. In the phenylpropanoid biosynthesis pathway (ko00940), CCR and POX were up-regulated in TH, while the expression level of PAL, C4H, CAD, and TOGT were enhanced in TL. In the alpha-linolenic acid metabolism (ko00592), ACAA was up-regulated in TH, and TL displayed a higher mRNA level of AOC and LOX. Overall, qRT-PCR results indicated that the genes performed the same expression trends as the RNA-seq data, suggesting the high level of reliability of the RNA-seq data.

Combined Analysis of Transcriptome and Metabolome
Integrated analysis indicated that 17 DEGs and four DMs participated in four KEGG reference pathways together, including phenylpropanoid biosynthesis and the alphalinolenic acid metabolism pathway (Table S4). In the phenylpropanoid biosynthesis, one DM (eugenol) and nine DEGs (PAL, C4H, 4CL, POX, CCR, CAD, TOGT, COMT, and CCoAOMT) were mapped together, and there was a high correlation between eugenol and PAL, C4H, 4CL, CCR, CAD, COMT, and CCoAOMT. In addition, DMs (Z)-3-hexenal and six DEGs (LOX, AOS, AOC, OPR, ACAA, MFP2) were mapped to the alpha-linolenic acid metabolism pathway. Here, a strong correlation between (Z)-3-hexenal and LOX existed (Figure 7). olenic acid metabolism pathway (Table S4). In the phenylpropanoid biosynthesis, one DM (eugenol) and nine DEGs (PAL, C4H, 4CL, POX, CCR, CAD, TOGT, COMT, and CCoA-OMT) were mapped together, and there was a high correlation between eugenol and PAL, C4H, 4CL, CCR, CAD, COMT, and CCoAOMT. In addition, DMs (Z)-3-hexenal and six DEGs (LOX, AOS, AOC, OPR, ACAA, MFP2) were mapped to the alpha-linolenic acid metabolism pathway. Here, a strong correlation between (Z)-3-hexenal and LOX existed (Figure 7). In the phenylpropanoid biosynthesis pathway, the relative content of eugenol in TH was much higher than that in TL, and its major regulatory genes located in the downstream and alternative pathways showed the same tendency, such as CCR, CCoAOMT, and COMT. However, the main upstream structural genes associated with eugenol biosynthsis showed decreasing tendency compared with TL, including the significant members in the common phenylpropanoid pathway, such as PAL, C4H, and 4CL. In order to study the screened structural genes related to the eugenol synthesis, expression analysis by qRT-PCR was performed on PAL and C4H in the three flowering stages of high-scent D. loddigesii, and the results showed that the expression level of PAL was the lowest in the blooming period, and the expression level of PAL was higher, but it was not significantly differentiated between the other two flowering periods; C4H has the highest expression In the phenylpropanoid biosynthesis pathway, the relative content of eugenol in TH was much higher than that in TL, and its major regulatory genes located in the downstream and alternative pathways showed the same tendency, such as CCR, CCoAOMT, and COMT. However, the main upstream structural genes associated with eugenol biosynthsis showed decreasing tendency compared with TL, including the significant members in the common phenylpropanoid pathway, such as PAL, C4H, and 4CL. In order to study the screened structural genes related to the eugenol synthesis, expression analysis by qRT-PCR was performed on PAL and C4H in the three flowering stages of high-scent D. loddigesii, and the results showed that the expression level of PAL was the lowest in the blooming period, and the expression level of PAL was higher, but it was not significantly differentiated between the other two flowering periods; C4H has the highest expression in the bud stage and shows the lowest expression in the blooming period ( Figure 8). The biosynthesis of (Z)-3-hexenal was relatively simple, the relative content of (Z)-3-hexenal and the expression of its upstream structural gene LOX both showed a downward trend in TH.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 10 of 13 in the bud stage and shows the lowest expression in the blooming period ( Figure 8). The biosynthesis of (Z)-3-hexenal was relatively simple, the relative content of (Z)-3-hexenal and the expression of its upstream structural gene LOX both showed a downward trend in TH.

Identification and Comparison of Metabolites between Two Varieties
Terpenoids constituted one of the most diverse class of volatile metabolites in plants [23]. Terpenes were the dominant components of floral scent in D. chrysotoxum, 15 ter-

Identification and Comparison of Metabolites between Two Varieties
Terpenoids constituted one of the most diverse class of volatile metabolites in plants [23]. Terpenes were the dominant components of floral scent in D. chrysotoxum, 15 terpenes were identified as scent components in its flowers [9]. Similarly, the aromatic composition of D. hancockii and D. officinale was mainly terpenes [16]. In the research, among the top ten volatile metabolites with the highest level in high-scented D. loddigesii (TH), 3-carene (11.87%) had the highest relative content with sweet smell. 3-Carene was a kind of monoterpene, which was identified widely in aromatic constituents of Dendrobium species, such as D. chrysotoxum, D. hancockii, and D. officinale [16,24]. Therefore, we inferred that 3-carene was the base content of floral scent in these Dendrobium species, and it was also the basic constituent of high-scented D. loddigesii aromatic scent. Geraniol was a specific aroma component of TH with rose smell, the relative contents and fold changes of geraniol and its derivative geranyl acetate in TH were significantly higher than those in TL [25]. Eugenol was a phenylpropane compound with a strong clove smell; its relative content and log 2 FC value were both top 10 in TH. Therefore, we speculated that geraniol, geranyl acetate, and eugenol were the main components of the characteristic aroma of high-scented D. loddigesii. Comprehensive results of relative content and log 2 FC value of volatile organic compounds (VOCs) showed that furans and (Z)-3-hexenal were important volatile components of TL, the aroma of these substances were mainly vegetable, bean, and grassy smell [25].
The relative content of terpenes (21.85%) played a dominant role in the volatile metabolites of TH, which was followed by esters (16.69%) and alcohol (8.68%). In contrast, heterocyclic compounds (6.84%), terpenes (6.43%), and aldehydes (6.06%) were the top three species with the highest relative content in light-scented D. loddigesii (TL). Meanwhile, the total relative content of the top ten VOCs in TL was only 22.41%, which was much lower than the corresponding parameter in the TH (52.63%). It could be seen that the relative content of terpenes in TH was much higher than that of other metabolites, but the relative content of different metabolites in TL were more balanced. In addition, three important monoterpenes species (geraniol, p-menthene, and o-cymene) were missing in TL components compared with THs. Perhaps the difference in the relative content and types of terpenes determined its weak flower fragrance in the light-scented D. loddigesii.

Association between Differential Metabolites and Differential Expressed Genes
GPPS was generally considered as a structural gene in terpenoid backbone biosynthesis; it was responsible for the head-to-tail condensation of one dimethylallyl diphosphate (DMAPP) with one isopentenyl diphosphate (IPP) molecule to form geranyl diphosphate (GPP), which is the precursor of monoterpenes [10]. In the TH, the expression trend of GPPS was consistent with the abundance of monoterpenoids in floral scent. Nevertheless, the integrated analysis indicated that the downstream pathway of GPPS, monoterpenoid biosynthesis (ko00902), showed no strong correlation to DMs except for geraniol. Further analysis found that many screened monoterpenoids DMs from two varieties of aromatic D. loddigesii could not be mapped into the KEGG reference pathway, because their related biosynthesis and metabolism pathways have not been studied clearly, including 3-carene, p-menthene, and o-cymene. In addition, the types of terpenoids in D. loddigesii were very different to those common components seen in other orchids, such as limonene, myrcene, linalool, and ocimene [26,27].
LOX is a key branchpoint enzyme leading to the biosynthesis of short-chain aldehydes and alcohols in alpha-linolenic acid metabolism (ko00592). In the pathway, linoleic acid could be formed to 9-hydroperoxy and 13-hydroperoxy intermediates under the oxidation of LOX, and these intermediates were further converted to c6/c9 aldehydes and alcohols such as (Z)-3-hexenal with a 'fresh green' smell; this main branch pathway is called the 'LOX pathway' [28]. Meanwhile, 13-hydroperoxy intermediates could be catalyzed by a series of enzymes such as AOC, allene oxide synthase (AOS), and 12-oxophytodienoate reductase (OPR) to generate jasmonate and its derivatives; this branch pathway was called β-oxidation [29]. Combining the results of metabolomic and transcriptomic analysis, the relationship between (Z)-3-hexenal and LOX could be established. The expression of LOX and AOC genes was all increased, and the downstream metabolites of LOX, (Z)-3-hexenal, and (E,E)-2,4-nonadienal were accumulated more in TL. Moreover, (Z)-3-hexenal has not been reported in the study of volatile components in Dendrobium species before. It indicated that the alpha-linolenic acid metabolism pathway and its related genes played a significant role in the formation of volatile components in light-scented D. loddigesii.
Phenylpropanoid compounds were the second-largest class of plant VOCs; the first committed step in the biosynthesis of the majority of phenylpropanoids is catalyzed by-PAL, which deaminates phenylalanine (Phe) to trans-cinnamic acid (CA) [10]. Then, CA is catalyzed by a series of methyltransferases and acyltransferases to form the volatile phenylpropenes, such as eugenol and isoeugenol. PAL and its downstream C4H, 4CL genes constituted the common phenylpropanoid pathway, which played a significant role in the synthesis of phenylpropanoids [30]. Nevertheless, the expression change trends of PAL, C4H, and 4CL were opposite to that of the relative content of eugenol in TH. To further explain this inconsistency, the expression level of PAL and C4H in the three flowering stages of the TH sample were analyzed using a real-time PCR method. The results showed that the expression levels of PAL and C4H in the flowering stage were much lower than those of the bud stage and half-open stage. It implied that the biosynthesis and release of phenylpropanoid volatiles might not be synchronized in TH.

Conclusions
The aroma compositions and characteristics of high-scented and light-scent D. loddigesii were determined by GC-MS based metabolomics, 38 differential metabolites (DMs) were identified, and several components including terpenes, especially monoterpenes, were suggested to be responsible for the aromatic difference between two Dendrobium varieties. Transcriptomic analysis was performed on the flowers of two aromatic varieties, 24,570 DEGs were obtained, and seven of the KEGG pathways were involved in the development of aromatic differences. Association analysis indicated that the most significant pathway was 'phenylpropanoid biosynthesis', which enriched the most DMs and DEGs. The expression trends of DEGs and mapped metabolites suggested that the biosynthesis and release of phenylpropanoid volatiles might not be synchronized in high-scented D. loddigesii. These results provide references for exploration of the perfume formulation and biosynthesis of aroma in D. loddigesii flowers.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/app12031262/s1, Figure S1: GO classification of unigenes in D. loddigesii floral transcriptome, Figure S2: KOG categorization result of unigenes in D. loddigesii floral transcriptome, Figure S3: GO assignments of differential expressed genes in D. loddigesii floral transcriptome, Table S1: Sequences of primers related to floral fragrance synthesis and internal reference of D. loddigesii, Table S2: Differential metabolites between two D. loddigesii species, Table S3: Statistic of aroma-related KEGG pathways enrichment, Table S4: Associated conditions between differential metabolites and differential expressed genes.