The wild C. faberi collected in this study were cultivated in the greenhouse at Forest Orchid Garden in Fujian Agriculture and Forestry University (Fuzhou, Fujian Province, China) under natural light and temperatures. The temperature was about 20–25°C. Flowers, leaves, and pseudobulbs of C. faberi were sampled at the flowering stage. The buds about 1cm before anthesis, semi-open flowers about 3cm, and fully open flowers were also used in this study. All samples of C. faberi were frozen in liquid nitrogen for storage at −80°C until use.

Two domains – PF01397 representing the TPS N-terminal domain and PF03936 representing the TPS C-terminal domain from PFAM¹ – were used as queries to search the C. faberi protein database (El-Gebali et al., 2019). The C. faberi genome data will be published separately. An HMM search (built-in Tbtools) was used in this study with an e-value cut at 10⁻³. To avoid missing potential TPS genes, TPS sequences from A. thaliana in the TAIR database² were also used to screen the C. faberi protein database using BLASTP (built-in Tbtools; Chen et al., 2018). The candidate TPS genes were checked manually by Pfam to verify putative full-length TPS genes, and TPS genes lacking either PF03936 or PF01397 were excluded. The grand average of hydrophobicity (GRAVY), molecular weight (Mw), isoelectric points (pI), aliphatic index (AI), and instability index (II) of the TPS proteins were predicted by the ExPASy database (Artimo et al., 2012).³ Subcellular localization was predicted by Plant-mPloc (Chou and Shen, 2010),⁴ AtSubP (Kaundal et al., 2010),⁵ and Ploc-mPlant (Cheng et al., 2017).⁶ Terzyme⁷ and BLATP⁸ were used to predict gene function (Priya et al., 2018).

Conserved motifs in the C. faberi TPS sequences were employed and analyzed by MEME software⁹ with default parameters (Bailey et al., 2009). We identified 20 motifs in this study. The exon-intron structure of the sequences was determined using GSDS software (Hu et al., 2015).¹⁰

The transcriptomes of TPS sequences from P. equestris and A. shenzhenica were downloaded from their genome databases. TPS sequences from S. moellendorffii were downloaded from NCBI.¹¹ TPS sequences from Picea abies were downloaded from UniProt,¹² and TPS protein sequences from A. thaliana and Oryza sativa were downloaded from Phytozome.¹³ All these sequences were aligned with MAFFT (Rozewicki et al., 2019). The maximum likelihood (ML) method was used for the phylogenetic tree, which was constructed with RAxML on the CIPRES Science Gateway web server (RAxML-HPC2 on XSEDE; Miller et al., 2011). Bootstrap values were 1,000 replicates with the JTT model. The most appropriate protein evolution model for the alignment was predicted by ProTest (Darriba et al., 2011). The generated tree was redrawn and annotated by EVOLVIEW (He et al., 2016).¹⁴ The sequences of the CfTPS proteins used in this study are listed in Supplementary Table S7.

Tbtools software extracted the promoter sequences and 2,000bp regions upstream of 32 CfTPS genes (Chen et al., 2018). Afterward, the online software PlantCARE¹⁵ was used to identify the cis-active regulatory elements in the promoter regions of the CfTPS genes (Lescot et al., 2002).

Gene pairs with similar genetic relationships were selected based on the phylogenetic tree. DNAMAN software was used to select the gene pairs with a consistency greater than 60%. Tbtools software was then used to calculate Ka (non-synonymous rate), Ks (synonymous substitution), and Ka/Ks (evolutionary constraint) values. Divergence time (T) was calculated by using the formula T=Ks/ (2×9.1×10⁻⁹)×10⁻⁶ million years ago (Mya; Zhang et al., 2018). In general, Ka/Ks<1.0 represents purifying or negative selection, Ka/Ks=1.0 represents neutral selection, and Ka/Ks>1.0 represents positive selection (Zhang et al., 2006).

An RNA sequencing transcriptome database of leaves, pseudobulbs, petals, sepals, labellums, and gynostemium was established to study the expression patterns of CfTPS genes. Fragments per kilobase of transcript per million mapped reads (FPKM) values of CfTPS genes were used to evaluate translation abundance. DESeq was used to conduct gene differential expression analysis, and gene ontology (GO) classification analysis was performed based on the differentially expressed gene analysis. The heatmaps of CfTPS expression patterns were drawn by Tbtools software, and the color in the heatmap was expressed as the log2-transformed expression levels of each CfTPS gene (Chen et al., 2018).

RNA was isolated from flowers of C. faberi at the flowering stage using the Biospin Plant Total RNA Extraction Kit (Bioer Technology, Hangzhou, China). First-strand DNA was synthesized with TransScript® All-in-One First-Strand cDNA Synthesis SuperMix for quantitative PCR (qPCR; TransGen Biotech, Beijing, China). TransScript® All-in-One First-Strand cDNA Synthesis SuperMix for qPCR was also used to remove genomic DNA. The real-time reverse transcription quantitative PCR (RT-qPCR) primers of CfTPS were designed by Primer Premier 5 software and can be found in Supplementary Table S12. Primer specificity was confirmed using Primer-BLAST on the NCBI website.¹⁶ PerfectStart™ Green qPCR SuperMix (TransGen Biotech, Beijing, China) was used for RT-qPCR analysis. The C. faberi reference gene GAPDH (GenBank Accession: JX560732; Tian et al., 2020) was used as the internal control and quantified by the 2^-△△CT method (Livak and Schmittgen, 2001). There were three biological replicates in the RT-qPCR analysis.

The ORF of CfPS18 was synthesized and ligated into pET28a vector. Then, the recombinant plasmid was transformed into Escherichia coli BL21 (DE3) pLysS cells (Transgen, China). Primers are given in Supplementary Table S12. The positive clones were incubated with shaking at 200rpm, 37°C until OD₆₀₀=0.6, and then at 37°C for 3h with shaking and 0.1mM IPTG. Recombinant CfTPS18 enzyme was induced at 16°C for 16h with 0.1mM IPTG. The precipitate was resuspended in extraction buffer [50mM NaH₂PO₄, 500mM NaCl, 10% (v/v) glycerol, and pH 7.0] and disrupted with a sonicator at 200W for 60s. The protein was purified with Ni-NTA agarose (Clontech). Purified CfTPS18 protein was examined by SDS-PAGE using Tris-HCl buffer (pH 7.5). The recombinant protein was performed in assay buffer [25mM HEPES, 10mM MgCl₂, 100mM KCl, 5mM DTT, 10% (v/v) glycerol, and 30μM GPP], and 30μg of CfTPS18 protein at pH 7.2 and 30°C for 1h²⁴.

The volatiles were exposed to SPME fiber with 50/30μm DVB/CAR/PDMS (divinylbenzene/carboxen/polydimethylsiloxane; Supelco Co., Bellefonte, PA, United States). The extract was analyzed using a gas chromatograph (Agilent 6,890N) and mass spectrometer (Agilent 5975B, Santa Clara, CA, United States) outfitted with a silica capillary column (DB-5MS; 60m×0.25mm×0.25μm). The temperature program was as follows: 55°C for 2min, 3°C min⁻¹ up to 80°C, 5°C min⁻¹ up to 180°C, 10°C min⁻¹ up to 230°C, and finally, 20°C min⁻¹ to 250°C. The ion source temperature was 230°C, and the electron energy was 70eV. The GC-MS interface zone temperature was 250°C, and the scan range was 50–500m/z. Reactions only added with GPP were used as the blank control. There were three biological replicates for the experiment. The retention time was compared with the NIST Mass Spectral Library.

To retrieve the TPS genes in C. faberi, two domains, PF01397 and PF03936 in the PFAM, were used to search the C. faberi protein database (El-Gebali et al., 2019). BLASTP (built-in Tbtools) was also used to screen the C. faberi protein database (Chen et al., 2018). After removing artefacts, 32 TPS genes were obtained. The deduced proteins of these genes were in the range of 115 for CfTPS1 and CfTPS2 to 902 amino acids for CfTPS12 and had predicted molecular weights (Mw) in the range of 10.24 for CfTPS12 to 103.89 KDA for CfTPS31. The theoretical isoelectric point (pI) values were in the range of 4.86 for CfTPS1 and CfTPS2 to 9.28 for CfTPS26, and the deduced grand average of hydrophilic (GRAVY) values were in the range of −0.358 for CfTPS26 to 0.073 for CfTPS6, suggesting that most CfTPS proteins were hydrophilic except for CfTPS1, CfTPS2, and CfTPS6. Additionally, the aliphatic index (AI) of CfTPS-deduced proteins was in the range of 87.09 for CfTPS21 to 107.22 for CfTPS6, and the instability index (II) was in the range of 28.62 for CfTPS1 and CfTPS2 to 52.36 for CfTPS5. To retrieve information on the subcellular localization of CfTPS proteins, three predictors were used in this study: AtSubp, Plant-mPLoc, and pLoc-mPLant (Chou and Shen, 2010; Kaundal et al., 2010; Cheng et al., 2017). The results showed that 16 CfTPS proteins were marked in the chloroplast, and 16 CfTPS proteins were marked in chloroplast or cytoplasm, indicating that these three predictors produce different results and need to be further analyzed. We also annotated 32 CfTPS genes using BLASTP¹⁷ and Terzyme software (Supplementary Table S4; Priya et al., 2018). The results showed that 11 CfTPSs were annotated sesquiterpene synthases, 15 CfTPSs were annotated monoterpene synthases, and six CfTPSs were annotated diterpene synthases. The secondary structure predicted by the SOPMA program revealed that the average of α-helices, random coils, extended strands, and β-turns comprised 70.86, 21.37, 4.35, and 3.63% of the structure, respectively (Geourjon and Deléage, 1995). The results showed that α-helices were predominant in all CfTPS proteins (Table 1).

To understand the intron-exon structure of CfTPS genes, we analyzed TPS gene structure with GSDS software (Hu et al., 2015). The exons in CfTPS genes ranged in numbers from 2 to 15, and the results showed that most of the genes in the same category contained a similar intron-exon structure.

To further analyze the motifs of the CfTPS genes, we identified 20 motifs using MEME software (Bailey et al., 2009). The numbers of C. faberi TPS motifs ranged in length from 4 to 16. CfTPS3 and CfTPS23 contained the most motifs, with 16, while CfTPS1, CfTPS2 had only four motifs. Motif 3 can be found in all CfTPS proteins except CfTPS9 and CfTPS16. Motif 4 was also the most common CfTPS protein (28/32). Twenty-five CfTPS genes contained the RRX₈W motif (motif 1), and 30 CfTPS genes contained the DDxxD motif (motif 2; Figure 2C). Accordingly, different clusters have different forms of motifs. The same cluster’s CfTPS proteins generally contained similar motifs. These results of intron-exon structure and motif analysis verified the closeness of the phylogenetic tree in C. faberi (Figure 2).

A phylogenetic tree was constructed to analyze the evolutionary patterns of the CfTPS genes. Thirty-two CfTPS genes were used, and TPS protein sequences from six species (A. shenzhenica, P. equetris, O. sativa, A. thaliana, P. abies, and S. moellendorfii) were used. The maximum likelihood (ML) method was used for the phylogenetic tree, which was constructed with RAxML on the CIPRES Science Gateway web server (RAxML-HPC2 on XSEDE; Miller et al., 2011). Bootstrap values were 1,000 replicates with the JTT model. The phylogenetic tree indicated that TPS proteins belonged to seven categories. This classification result was consistent with a recent study: TPS-a, TPS-b, TPS-c, TPS-d, TPS-e/f, TPS-g, and TPS-h (Chen et al., 2011). Thirty-two CfTPS proteins belonged to three categories according to the phylogenetic tree (Figure 3): TPS-a, TPS-b, and TPS-e/f. Of these three categories, the TPS-a and TPS-b clades contained the most members and were the most expanded categories, with 13 and 15 genes, respectively, and were consistent with other plant species, such as A. thaliana, C. sinensis, V. planifolia, D. catenatum, P. equestris, and D. officinale (Aubourg et al., 2002; Yu et al., 2020; Zhou et al., 2020; Huang et al., 2021). The remaining four TPS genes belonged to the TPS-e/f subfamily.

We aligned the multiple sequences using MAFFT to analyze the TPS RRX₈W motif in the N-terminus, DDxxD, and NSE/DTE motifs in the C-terminus (Rozewicki et al., 2019). The alignment showed that almost all the CfTPS proteins in the TPS-a and TPS-b subfamilies had the highly conserved aspartate-rich motif DDxxD, except CfTPS9 and CfTPS16 in the TPS-a subfamily and CfTPS17, CfTPS27, and CfTPS28 in the TPS-b subfamily (Jiang et al., 2019). TPS genes in TPS-a and TPS-b had variations in the RRX₈W and RxR motif. However, in the TPS-a and TPS-b categories, the RRX₈W motif was not found in nine CfTPS genes. The RRX₈W motif was absent in the TPS-e/f subfamily, and the NSE/DTE motif was found in 16 CfTPS genes (Figure 4). Among them, DDxxD and NSE/DTE are important in the metal-dependent ionization of the prenyl diphosphate substrate, and the RRX₈W motif is essential in the cyclization of monoterpene synthase (Bohlmann et al., 1998; Chen et al., 2011; Jiang et al., 2019). In addition, the members in the TPS-a subfamily detected in both dicot and monocot plants encode only sesquiterpenes (Jiang et al., 2019). The secondary “R” in this family is not conserved (Martin et al., 2010). TPS-b encodes monoterpenes containing the conserved RRX₈W motif (Chen et al., 2011; Jiang et al., 2019).

To retrieve the potential function of CfTPS genes, we obtained a 2,000-bp region upstream of the 32 CfTPS genes and analyzed them using the online software Plantcare (Lescot et al., 2002). In total, we found 784 cis-acting regulatory elements in the promoter areas of CfTPS genes. We classified them into three categories according to the function of these elements: plant growth and development, stress responsiveness, and phytohormone responsiveness. The plant growth and dependent category (166/784) contained nine cis-acting regulatory elements and consists of AAGAA motifs, As-1 elements, A1-rich elements, etc. Most of them had As-1 elements (66/166), which are related to shoot expression. The stress responsiveness category (216/874) contained ARE, DRE, LTR, ST-RE, ABRE, etc. STRE (48/216) is an essential element in the promoter that is related to stress. Interestingly, most of the cis-elements (401/784) were in the phytohormone responsiveness category, which contained CGTCA motifs, EREs, MYC motifs, TGAGG motifs, etc. Among them, most cis-elements were MYC (137/401), which is associated with MeJA responsiveness (Dombrecht et al., 2007). The results indicated that the CfTPS gene expression patterns might be regulated by MeJA treatment (Figure 5).

The Ka/Ks ratio can show positive selection (Ka/Ks>1), negative or purifying selection (Ka/Ks<1), and neutral selection (Ka/Ks=1) during the evolution (Zhang et al., 2006). In this study, 13 gene pairs with similar genetic relationships were selected for Ka/Ks calculation. The results showed that the Ka/Ks ratios of 12/13 CfTPS genes were less than one, indicated that most CfTPS genes underwent negative selection (Table 2). The divergence time of 13 CfTPS gene pairs was in the range of 0.75 for CfTPS1 and CfTPS2 to 60.46 for CfTPS3 and CfTPS10.

An RNA sequencing transcriptome database of leaves, pseudobulbs, petals, sepals, labellums, and gynostemium was established to study the expression patterns of CfTPS genes. Eighteen genes were expressed in both leaves and flowers, and 20 genes were expressed in the labellums. Twenty-five genes were found in sepals and 24 in gynostemium. Four genes were not found to be expressed in any of the tissues. CfTPS12, CfTPS15, and CfTPS23 exhibited a high level of expression in leaves and pseudobulbs. Notably, CfTPS3, CfTPS12, CfTPS15, CfTPS17, CfTPS18, CfTPS28, and CfTPS31 displayed high transcript abundance in floral organs, suggesting that these TPS genes might be related to flower scent in C. faberi (Figure 6).

Gene ontology classification analysis was performed based on the differentially expressed gene analysis. According to the classification results, molecular function contained most genes, and genes were mostly enriched in lyase activity, magnesium binding, and terpene synthase activity (Figure 7).

Real-time reverse transcription quantitative PCR was performed to investigate the expression patterns of CfTPS genes in flowers at three floral developmental stages. Four putative functional genes, CfTPS12, CfTPS18, CfTPS23, and CfTPS28, which are highly expressed in floral organs, were used. They were expressed differently among the three floral stages. Interestingly, these genes were mainly expressed in the full flowering stage (Figure 8). Notably, CfTPS18 had the highest transcript levels during the floral developmental stages, consistent with the expression patterns displayed in Figure 5. Gene annotation indicated that CfTPS12, CfTPS18, CfTPS23, and CfTPS28 likely encoded monoterpene synthases and clustered in TPS-b subfamily (Table 1).

To investigate the enzyme activity, the full-length ORF sequence of CfTPS18 was cloned to vector pET28a and ectopically expressed in E. coli. Recombinant CfTPS18 enzyme was induced with 0.1mM isopropyl-β-d-galactopyranoside and purified with Ni-NTA agarose (Clontech). SDS-PAGE analysis of CfTPS18 recombinant protein can be found in Supplementary Figure S1. After using GPP as substrate in the reactions, the products were analyzed by GC-MS. The result indicated that CfTPS18 could convert GPP into β-myrcene, geraniol, and α-pinene which were validated by the NIST Mass Spectral Library (Figure 9). The blank control with only GPP added to the reactions could not produce monoterpenes. Thus, CfTPS18 was considered a monoterpene synthase.