Duplication and selection in β-ketoacyl-ACP synthase gene lineages in the sexually deceptive Chiloglottis (Orchidaceace)

Keywords: Chiloglottis; chiloglottone; sexual deception; pollination; transcriptome; fatty acid; gene duplication; selection; evolution

Gene duplication is a prominent feature of plant genomes and many authors have hypothesized its importance in enabling evolutionary novelty within plants ([17]). Whole genome duplication (WGD) via an increase in ploidy, and tandem duplication via unequal crossing over between similar alleles, accounts for the majority of plant gene duplicates ([45]). Other mechanisms such as transposon-mediated duplication, segmental duplication and retroduplication are less prevalent. Despite their potential importance in the evolution of new functions, the majority of duplicated genes are often pseudogenized/lost due to the accumulation of mutations that render them non-functional ([17]). However, other duplicates are either retained by selection on existing functions (e.g. gene dosage effects and sub-functionalization), novel functions (e.g. neo-functionalization) or both ([45]).

Selection on a few critical residues of duplicated genes encoding plant volatile enzymes is thought to have enabled the evolution of the many plant volatiles involved in mediating plant–insect interactions ([21]; [52]; [50]). For example, the terpene synthase gene family is highly duplicated in plants ([51]). In a more specific example, the sesquiterpene (E)-α-bergamotene in the wild tobacco Nicotiana attenuata serves dual roles – being an indirect defence against Manduca sexta larvae in the leaves and a pollinator attractant of the adult moth in the flowers ([71]). Genetic evidence indicates that the sesquiterpene synthase gene (NaTPS38) responsible has evolved from a monoterpene synthase under positive selection ([71]).

In the sexually deceptive orchid genus Ophrys, blends of (Z)-9 and (Z)-12 alkenes are involved in the pollinator attraction of Ophrys sphegodes while (Z)-7 alkenes are involved in Ophrys exaltata ([65]). Phylogenetic analysis of orchid stearoyl-acyl carrier desaturase (SAD) genes suggests that Ophrys homologues (SAD1–6) evolved via gene duplication and the SAD2 clade in particular was under strong positive selection ([54]; [65]). Further work has shown that SAD2 diverged from an ancestral housekeeping desaturase enzyme activity by catalysing the introduction of specific double-bonds in intermediates, from which respective alkenes are hypothesized to be synthesized ([54]).

Although one of the most diverse plant families in the world, genomic resources for orchids still lag behind model and agriculturally important plant systems ([46]; [61]). Nonetheless, the genomes of four orchid species, namely Phalaenopsis equestris, Dendrobium catenatum, Apostasia shenzhenica and Gastrodia elata ([6]; [69], [70]; [68]), are now available. A growing number of orchid transcriptomes drawn from across diverse lineages are also becoming available ([55]; [12]; [2]; [9]; [64]; [70]; [25]; [59]). Together, the combination of genomes and transcriptomes provide an excellent resource to prioritize candidate genes and understand the evolution of gene families involved in the biosynthesis and regulation of plant volatiles in orchids via comparative genomic analysis ([61]).

A first step to elucidate the molecular basis of target plant volatiles often involve transcriptome comparisons between active (scent-producing) and non-active floral tissues/species. Such studies have indeed aided the identification of the following genes in orchids: (1) SAD2 and SAD5 genes involved in alkene biosynthesis in Ophrys as previously described above ([54]; [56]); (2) eugenol synthase genes involved in eugenol/isoeugenol biosynthesis in three Gymnadenia orchids ([26]); (3) geraniol synthase, alcohol dehydrogenase and geranial reductase genes involved in the complete biosynthesis of (S)-β-citronellol in Caladenia plicata ([64]); (4) a vanillin synthase in Vanilla planifolia ([20]); and (5) the regulation of floral mono/sesquiterpene biosynthesis in Phalaenopsis by a suite of transcription factors (TFs) including bHLH4, bHLH6, bZIP4, ERF1 and NAC1 ([10]).

The Australian sexually deceptive (SD) Chiloglottis orchids employ 2,5-dialkylcyclohexan-1,3-dione(s), with six known variants (named chiloglottones 1–6), to attract specific male wasp pollinators ([49]; [48]). Phylogenetic analysis using internal transcribed spacer and chloroplast DNA sequences revealed that Chiloglottis trapeziformis, which belongs to the Formicifera clade of the genus, employs chiloglottone 1 to attract its primary male pollinator, Neozeleboria cryptoides. The same compound is employed by Chiloglottis seminuda of the Reflexa clade to attract an undescribed species of Neozeleboria sp. (proxima2) ([49]). In all Chiloglottis species investigated, chiloglottone production occurs specifically in the densely clustered 'insectiform' calli structure (CAL) on the labellum. In addition, the glandular sepal tips (GS) are a further source of chiloglottone production in C. seminuda and other members of the Reflexa clade ([49]; [16]). Unexpectedly, chiloglottone production is also dependent on UV-B radiation, although the mechanism remains to be elucidated ([16]; [1]; [59]).

Chiloglottone biosynthesis is predicted to involve the fatty acid (FA) biosynthetic and/or degradative (β-oxidation) pathway intermediates ([18]; [4]). These precursor intermediates are ubiquitous in the plastids as acyl-carrier-protein (ACP) derivatives or the peroxisomes as coenzyme A (CoA) derivatives. To make these intermediates available for chiloglottone biosynthesis in an otherwise iterative pathways, mid-cycle termination of FA biosynthetic and/or β-oxidation may also be crucial ([4]). Thus far, evidence in support of this hypothesis has accumulated from studies of C. trapeziformis, the species in which chiloglottones were first discovered.

Transcriptome analysis of C. trapeziformis flowers has revealed two key gene expression trends linked to the distribution profiles of chiloglottone ([59]): (1) the majority of fatty acid biosynthesis and β-oxidation genes are exclusively expressed (or highly induced) in the active (chiloglottone-producing) tissue; and (2) some genes, particularly KETOACYL-ACP SYNTHASE (KAS), FATTY ACYL-ACP THIOESTERASE (FATB) and ACYL-COA OXIDASE (ACX), are also developmentally regulated in the active tissue during the transition from young buds to flowers. Furthermore, the specific KAS inhibitor, cerulenin, significantly reduces chiloglottone production in three chiloglottone-producing Chiloglottis species tested: C. trapeziformis, C. valida and C. aff. valida ([59]).

To advance our understanding of chiloglottone biosynthesis further, in this study we have strategically chosen Chiloglottis seminuda as an alternative model. This species is phylogenetically distinct from C. trapeziformis – belonging to a different clade of the genus (Reflexa vs. Formicifera) ([49]), but shares chiloglottone 1 as its active semiochemical. Furthermore, both the labellum callus and the sepal tips are the site of chiloglottone 1 production, thus providing us with an additional tissue contrast not available in C. trapeziformis.

Based on the hypotheses and accumulating evidence outlined above, we performed a combination of tissue-specific transcriptome analysis, enzyme inhibition experiments, and molecular phylogenetic and evolutionary analysis to address three specific questions. (1) Are FA metabolic pathway genes differentially regulated in the active (CAL and GS) compared to the non-active tissues in C. seminuda? (2) Does the inhibition of FA biosynthesis impact chiloglottone biosynthesis in C. seminuda in a tissue-specific manner? (3) Are there signatures of positive or relaxed purifying selection at genes potentially associated with the evolution of chiloglottone biosynthesis in Chiloglottis, and if so do they involve the FA metabolic pathway genes?

Whole C. seminuda plants with single flowers were sampled from a colony growing naturally near Mount Werong within the Blue Mountains National Park in NSW, Australia. To ensure total chiloglottone depletion prior to inhibition experiments, plants were collected in the field and kept in a growth chamber lacking in the UV spectrum under the following conditions: acclimatization period, ≥5 d; day–night cycle, 12 h; temperature, day (20 °C) and night (15 °C); light parameters, 300 µmol m−2 s−1 (<400 nm).

In the field, freshly picked open flowers were carefully dissected into the calli (active), glandular sepals (active), and a pooled sample of non-glandular sepals and petals (non-active), then quickly snap-frozen in liquid N2. Three biological replicates were used for each treatment (i.e. tissue type) in which each contains tissues pooled from eight individual plants. Subsequently, the samples were shipped to Michigan as frozen tissue in an N2 Vapour Shipper (MVE). Extraction of total RNA was conducted using the Spectrum Plant Total RNA kit (Sigma-Aldrich). Paired-end library construction was performed using the SMART cDNA synthesis kit (Clontech) according to the manufacturer's instructions. All libraries were pooled and sequenced on a single lane on an Illumina HiSeq2500 at the Purdue Genomics facility (Purdue University; http://www.genomics.purdue.edu/). All raw sequence reads have been deposited in the NCBI Sequence Read Archive under the BioProject accession PRJNA486025 and SRA study accession SRP157949.

Raw paired-end (100 nt) reads were trimmed and quality filtered to remove low-quality reads and contaminating sequences (adapters and primers) using Trimmomatic v0.35 ([5]) as previously described ([60]). Surviving reads were pooled and de novo transcriptome assembly was performed using a combination of single (Trinity v2.5.1; [24]) and multi-k-mer (SOAPdenovo-Trans v1.04; [63]) strategies. A default k-mer size of 25 was used for Trinity while k-mer values of 25, 31, 51, 71 and 91 were used for SOAPdenovo-Trans. Resolution of multiple assembly redundancy and the production of high-quality and accurate gene sets were made using the EvidentialGene tr2aacds pipeline (Don [22]). Briefly, removal of duplicate sequences, clustering of sequences and identification of regions of local similarity among sequences was achieved using exonerate ([57]), CD-HIT v4.6.8 ([37]) and various tools in BLAST v2.2.29 ([7]). Contigs (hereafter transcripts) classified as 'primary' sets were used as the final reference transcriptome.

Alignment of filtered reads towards the refence transcriptome was performed using HISAT2 v2.1.0 ([35]) using default parameters and counts summarized using FeatureCounts ([38]). Raw (un-normalized) counts were loaded into DESeq2 v1.2 ([41]) and a differential expression test between treatments was inferred by the use of negative binomial generalized linear models. An absolute log2 fold change (log2FC) > 0.5 and false discovery rate (FDR) < 0.05 indicates significant differential expression. The final transcript abundance was expressed as fragments per kilobase of transcript per million mapped reads (FPKM). Functional annotation (e.g. gene ontology, GO; InterPro domain, IPR; and MapMan ontology, BIN) of the reference transcriptome was performed using two plant-specific pipelines: TRAPID ([58]) and Mercator ([40]). The parameters for TRAPID were blast_db_type = CLADE, blast_db = Monocots, e_value = 10e-5, gf_type = HOM, func_annot = gf_besthit, while the parameters for Mercator were left at default, except having the 'ORYZA' settings enabled. Enriched GO and MapMan categories (FDR < 0.05) were determined based on a hypergeometric distribution adjusted with FDR for multiple hypothesis correction in R as previously described ([60]). Emphasis was given to high-level categories describing general plant biological processes (i.e. GO slim and BIN depth ≤ 2).

Chemical inhibition experiments were performed on cut flowers of growth-chamber acclimatized (chiloglottone-depleted) whole plants using 100 µm and 500 µm cerulenin (Cayman Chemical) and a solvent control (without cerulenin), as previously described ([1]; [59]). Briefly, stalks of cut flowers were immersed into either the treatment or control buffer, sealed with parafilm at the top of each test tube, and held in place for 3 d (three flowers per test tube). At the end of the 3-d incubation period, the inhibitor-treated and control plants were illuminated with UV-B light for 2 h using a custom made light-box to induce chiloglottone production. The calli and sepal tips were excised and assayed for chiloglottones as previously described ([16]). Briefly, excised tissues were washed in 100 µL of 5-methyl-1,3-cyclohexanedione-spiked (20 ng µL–1, internal standard) dichloromethane for 3 min. Each sample was analysed on an Agilent Technologies gas chromatograph (model 6890N) coupled with a Mass Selective Detector (Agilent Technologies) equipped with a SGE BP21 column (30 m × 0.25 mm × 0.25 μm) connected directly to the MS detector. For each sample, 4 µL of extract was injected splitless into the inlet at 250 °C, the column was held at 40 °C for 1 min, then programmed at 10 °C min–1 to 230 °C and held for 15 min. Chiloglottone abundance was calculated based on the corrected percentage areas relative to the internal standard using the Agilent Technologies Chemstation software. A least 11 flowers were used in each treatment or control group. An analysis of variance (ANOVA) with a priori comparisons between means assessed by the Tukey–Kramer HSD test was performed to evaluate the outcomes of inhibition treatments in R (https://www.r-project.org/).

Coding sequences of the four sequenced orchids and 18 others predicted from transcriptome de novo assembly were retrieved from the relevant publications as follows. Genomes: Phalaenopsis equestris ([6]), Dendrobium catenatum ([69]), Apostasia shenzhenica ([70]) and Gastrodia elata ([68]). Transcriptomes: Neuwiedia zollingeri, Apostasia wallichii, Cypripedium formosanum, Cymbidium ensifolium, Cymbidium sinense, Phalaenopsis bellina and Phalaenopsis modesta from [9]; Neuwiedia malipoensis, Paphiopedilum armeniacum, Cypripedium singchii, Hemipilia forrestii, Holcoglossum amesianum and Gastrochilus calceolaris from [12]; Orchis italica from [11]; Ophrys sphegodes from [55]; Oncidium 'Gower Ramsey' from [8]; Caladenia plicata from [64]; and Chiloglottis trapeziformis from [59]). Coding sequences of Arabidopsis (TAIR10) and rice (IRGSP-1.0) were retrieved from Ensembl Plants (http://plants.ensembl.org/index.html).

To build a phylogeny for use in the downstream selection analysis, putative single-copy orthologous genes from 25 species (23 orchids and two non-orchid species) were first identified using OrthoFinder v2.2.7 ([14]). The final concatenated alignment matrix of the single-copy orthologous genes was prepared with PosiGene and a maximum-likelihood (ML) tree was constructed using IQ-TREE v1.6.6 ([44]) with the following settings: -m MFP and -bb 1,000 arguments to obtain the best-fit model from a wide range of phylogenetic models tested ([34]) and perform an ultrafast bootstrap test ([27]) with 1000 replications, respectively. For the ketoacyl-ACP synthase (KAS) phylogeny, high-confidence KASI and KASII homologues of Orchidaceae and several non-Orchidaceae species (i.e. maize, rice, tomato, soybean and grape) was identified using Arabidopsis KASI and KASII protein sequences as the query from the OrthoFinder run. Prior to tree construction, sequences were aligned with MUSCLE with default settings ([13]). An ML tree was constructed using IQ-TREE ([44]) with the same settings as described above.

Using the coding sequences of the 25 species listed above, a transcriptome-wide analysis of genes under selection, based on the estimation of the non-synonymous to synonymous substitution ratio, ω (also dN/dS), was made using CodeML ([66]) as implemented in PosiGene ([53]). Briefly, filtered codon alignments with one isoform/gene/species were used for branch-site tests of positive selection. The anchor (-as) and reference species (-rs) was set to Apostasia shenzhenica and the target (-ts) and context (-cs) species set to Chiloglottis seminuda and Chiloglottis trapeziformis to determine genes under selection in the Chiloglottis orchids (foreground branch). Genes were deemed to be under positive selection when ω was >1, under neutral selection when ω = 1, and under purifying selection when ω was <1. Positive selection on duplicate genes is predicted to create the potential for the evolution of new functions, while strong purifying selection is expected to maintain protein-coding function of existing genes. Selection analysis of a KASI phylogeny was performed with CodeML using the free-ratio model implemented in PAML v4.9 ([66]).

In this study, sequencing of three floral tissue-specific transcriptomes was performed for the chiloglottone-producing (active) labellum callus (CAL), glandular sepal tips (GS) and the non-active sepal remains and petals (SP) (Fig. 1A). Three biological replicates were represented for CAL and SP tissues while only two biological replicates were available for GS.

Graph: Fig. 1. Analysis of Chiloglottis seminuda tissue-specific floral transcriptome. (A) The floral structure of C. seminuda consists of the callus (CAL), labellum (LAB), glandular sepal (GS), petal (PTL) and column (CM). Image adapted from [16]. (B) Principal component analysis of the callus (CAL), glandular sepals (GS), and the non-active remainder of pooled non-glandular sepals and petals (SP) transcriptomes. CAL, GS and SP samples are depicted as light red, blue and purple empty circles. A, active; N.A, non-active. (C) Summary of differentially expressed transcripts between active and non-active floral tissue – CAL vs. SP (TS1) and GS vs. SP (TS2). The Venn diagram depicts the commonalities (intersection) and differences of differentially expressed transcripts (upregulated, red; downregulated, blue) between TS1 and TS2. Summary of enriched (FDR < 0.05) high-level (D) MapMan BIN and (E) GO SLIM categories describing general plant biological processes for upregulated (red) and downregulated (blue) genes identified in each comparison shown in (C). Red and blue circle colours depict enrichments for upregulated and downregulated genes of each contrast. Circle size corresponds to number of transcripts in each enriched category.

The Illumina sequencing yielded over 91.7 million paired-end reads (100 bp in length). After the removal of low-quality sequences ~83.2 million reads (90.7 % of the total sequenced reads) was used for de novo transcriptome assembly (Table 1). Trinity assemblies (fixed k-mer of 25) yielded over 300 000 contigs while SOAPdenovo-Trans assemblies yielded using k-mer of 25 and 31 yielded >1 million contigs and higher k-mer assemblies yielded over 94 000 (k-mer of 51), 66 000 (k-mer of 71) and 20 000 (k-mer of 91) contigs. A final set of 28 778 primary transcripts (reference transcriptome) identified using EvidentialGene were used for read mapping, count summarization, and other downstream transcriptome analysis (Supplementary Data Table S1). Functional annotation of the reference transcriptome using TRAPID and mefisto revealed that 17 609 (GO biological process), 14 904 (InterPro domain) and 10 490 (MapMan BIN) transcripts were successfully assigned to at least one functional term in each ontology class (Table S2). Focusing on the MapMan BIN categories, we observed a greater representation of transcripts involved in protein metabolism (,411), RNA metabolism (1843), signalling (989), transport (751) and cell organization (569) categories among others. Similar patterns were also observed when considering GOslim (BP) terms (Fig. S1).

Summary of RNA sequencing analysis metrics

1 CAL, callus; GS, glandular sepals; SP, sepal and petal remains.

2 *Million.

Principal component analysis of the reference transcriptome showed a clear separation of the samples based on tissue specificity (PC1: 43 % and PC2: 32 %) with all samples (CAL, GS, SP) forming discrete clusters along both axes (Fig. 1B). Tissue-specific comparisons between the non-active SP and active CAL (TS1: CAL-SP) and GS (TS2: GS-SP) tissues revealed 2703 (up: 1374; down: 1329) and 1424 (up: 565; down: 859) transcripts that were differentially expressed (|log2FC| > 0.5, FDR < 0.05), respectively. In contrast, the number of differentially expressed transcripts were much lower (up: 538; down: 327) between the active CAL and GS (TS3: CAL-GS) (Table S3). In addition, 153 transcripts were consistently down-regulated in TS1 and TS2 while 394 transcripts were consistently up-regulated (Fig. 1C).

Enriched (FDR < 0.05) Mapman BIN categories in the two contrasts revealed that lipid (FDRTS1 < 3.92 × 10−4, FDRTS2 < 1.79 × 10−2), secondary (FDRTS1 < 5.61 × 10−6, FDRTS2 < 5.56 × 10−3) and protein (FDRTS1 < 1.60 × 10−5, FDRTS2 < 1.01 × 10−3) metabolism were consistently enriched in upregulated genes while photosynthesis (FDRTS1 < 2.74 × 10−13, FDRTS2 < 2.47 × 10−61) and cell organization (FDRTS1 < 3.61 × 10−5, FDRTS2 < 4.29 × 10−3) terms were enriched for downregulated genes (Fig. 1D; Table S4). Similar trends were also observed with genes that belong within the down- and upregulated intersects. Downregulated genes in TS2 were also enriched for many other categories such as tetrapyrrole synthesis (FDR < 3.08 × 10−11), redox (FDR < 1.29 × 10−2), RNA (FDR < 4.08 × 10−4) and signalling (FDR < 3.21 × 10−4). Enrichment analysis using GO Slim also revealed parallel trends to that of enriched Mapman BIN categories in TS1, TS2 and their intersects (Fig. 1E, Table S4).

Tissue-specific comparisons (TS1, TS2 and TS3) highlighted many genes putatively involved in FA biosynthesis and degradation (β-oxidation) that were co-ordinately induced in active (CAL and GS) compared to the non-active SP tissues. The genes putatively involved in FA biosynthesis and elongation included one KETOACYL-ACP REDUCTASE (KAR), three KETOACYL-ACP SYNTHASE I (KASI-1, KASI-2A, KASI-2B), one KETOACYL-ACP SYNTHASE II (KASII) and one FATTY ACYL-ACP THIOESTERASE (FATB) while genes putatively involved in FA β-oxidation were two LONG-CHAIN ACYL-COA SYNTHETASE (LACS1 and LACS2), two ACYL-COA OXIDASE (ACX1 and ACX4) and one MULTIFUNCTIONAL PROTEIN (AIM1-1) (Fig. 2).

Graph: Fig. 2. Fatty acid (FA) biosynthesis and degradation pathway gene expression in the active and non-active tissues of Chiloglottis seminuda flowers. Differentially expressed genes encoding predicted (A) FA biosynthesis and (B) degradation enzymes are depicted. Transcript abundance is represented as fragments per kilobase of transcript per million mapped reads (FPKM) ± s.e. Inset depicts subsequent FA elongation and desaturation (from C16:0) reactions prior to export from the plastid. Blue, red and green bars depict CAL, GS and SP gene expression, respectively. Points connecting any two conditions depict statistically significant (FDR < 0.05, |log2FC| > 0.5) differences between treatments as determined by DESeq2 using raw (un-normalized) count data. The numerator and denominator conditions of each comparison are represented as n and d, respectively. Red, blue and grey connecting points depict upregulation, downregulation and the absence of differential expression, respectively.

Some genes also exhibited distinctive transcriptional profiles in active tissues, indicating induction potential and/or tissue specificity. For example, FATB, which was the most highly upregulated transcript in active tissues (especially GS), was >64-fold upregulated when compared to SP. In addition, KASI-2B, KASI-1, and LACS1 were also highly upregulated in active tissues (especially CAL) compared to SP, albeit to a lesser extent than FATB. Meanwhile, other transcripts such as KASI-2A, KAR, ACX1 and AIM1-1 were expressed 1.5- to 3-fold more in active tissues compared to SP (Fig. 2). There were only four instances where we found significant differential expression between active tissues: KASI-1, KASII and LACS1 were significantly higher in the CAL compared to GS while LACS2 was significantly lower in the CAL (Table S3).

Following a 2 h exposure of depleted flowers to UV-B light, chiloglottone 1 levels in flowers incubated with the KAS inhibitor, cerulenin, were significantly reduced (P < 0.01) compared to the solvent controls (Fig. 3). In the CAL, the mean level of chiloglottone 1 in the control was 1338 ± 147 ng per callus. Treatment with cerulenin at 100 µm and 500 µm resulted in a 27 % (P > 0.05) and 57 % (P < 0.01) decrease in chiloglottone 1 compared to the control. In the GS, the mean level of chiloglottone 1 in the control was 1729 ± 226 ng per callus. Treatment with 100 µm and 500 µm cerulenin resulted in a significant decrease in chiloglottone 1 compared to the control, with 36 % (P < 0.05) and 49 % (P < 0.05) reductions, respectively as determined by Tukey's HSD test.

Graph: Fig. 3. Chiloglottone amounts in the callus (CAL) and glandular sepals (GS) of Chiloglottis seminuda flowers following cerulenin treatment at 100 µm (CER100) and 500 µm (CER500) concentrations compared to the controls (CTL). The number of flowers used for CRE treatment and controls in the CAL and GS are as follows: CALCER100 = 17, CALCER500 = 11, CALControl = 22, GSCER100 = 20, GSCER500 = 11, GSControl = 21. Bars represent ± s.e. Asterisks indicate significant differences between treatments at P < 0.05 (*), P < 0.01 (**) based on one-way ANOVA and Tukey's HSD test.

A total of 25 315 orthologous gene groups (orthogroups) were identified using the OrthoFinder pipeline using translated protein coding sequences (CDS) of C. seminuda as a reference and 24 other species (22 orchids and two non-orchid species). Of these, 222 were deemed to be single-copy orthogroups. After further stringent filtering to remove gaps and unreliable alignments, a final concatenated alignment matrix of 201 putative single-copy orthogroups consisting of 37 344 sites (parsimony-informative: 11 977, singleton: 5787, and constant: 19 580 sites) was used to produce an ML tree (Fig. 4A).

Graph: Fig. 4. The evolution of Chiloglottis orchids as revealed by comparative genomics analysis. (A) A maximum-likelihood tree of 23 Orchidaceae (including C. seminuda and C. trapeziformis) and two non-Orchidaceae species (as an outgroup) based on 201 putative single-copy orthologous genes. The black downward triangle indicates the evolutionary branch analysed for positive selection signatures at a transcriptome-wide scale. Bootstrap support for each branch is depicted as coloured circles: 0–50 % (blue), 51–75 % (cyan), 76–90 % (green), 91–99 % (orange) and 100 % support (red). Light yellow (Epidendroideae), blue (Orchidoideae), green (Cypripedioideae) and purple (Apostasioideae) boxes indicate the various orchid subfamilies. (B) Boxplot depicting the distribution of all valid non-synonymous to synonymous substitution rate ratio (ω or dN/dS) values in the last common ancestor of Chiloglottis orchids. Posigene-inferred ω for fatty acid pathway orthologous genes are indicated as blue dots.

The ML phylogeny revealed four distinct clades, beginning with the Apostasioideae (four species: e.g. N. zollingeri, A. shenzhenica) as the basal clade, followed by the clade formed by Cypripedioideae (three species), Orchidoideae (six species) and Epidendroideae (ten species). As expected, the three Australian orchid species represented, Chiloglottis seminuda, Chiloglottis trapeziformis and Caladenia plicata, which belong to the monophyletic Diurideae ([28]), clustered together within the Orchidoideae.

A total of 2689 high-quality, filtered codon alignments with one isoform/gene/species were used for branch-site tests of positive selection on the terminal branch that corresponds to the Chiloglottis orchids (C. seminuda and C. trapeziformis) using PosiGene. The median non-synonymous to synonymous substitution rate ratio (ω or dN/dS) for all transcripts with alignments was 0.183 (Fig. 4B), thus indicating purifying selection for the majority of transcripts. However, 110 genes showed footprints of positive (ω ≥ 1) selection on the branch leading to Chiloglottis orchids. Positively selected genes were represented by many different functional categories such as cell wall metabolism, carbohydrate metabolism, lipid metabolism, protein metabolism, cell organization and transport, among others (Table 2, Supplementary Data Table S5).

List of selected genes under positive selection grouped by MapMan BIN functional category

• 3 The inferred ω value for each gene candidate is given in parentheses. Total PSG, total number of genes under positive selection in each BIN category; positively selected genes that are differentially expressed in at least one tissue contrast are indicated in bold.
• 4 *DE, number of PSGs that are differentially expressed.

The substitution ratios for the majority of FA metabolic pathway genes in Chiloglottis indicated most of these genes to be under purifying selection (Fig. 4B), and thus this negative selection should maintain the status quo functions of these genes. However, the ratio for FATA indicated strong positive selection (ω = 11.08) while the ratio for KASI-2 indicated some relaxation from strong purifying selection (ω = 0.38) when compared to the median ω observed (ω = 0.18) (Table S5), potentially enabling new functions to evolve. Beyond this intriguing selection analysis result, we have already noted that the expression of KAS genes was high in the active tissues (Fig. 2). Therefore, we investigated the phylogeny and evolution of KAS genes in more detail (see below).

An ML phylogeny constructed using KASI and II sequences from 14 Orchidaceae, maize, rice, tomato, soybean and grape revealed that KASI and II sequences formed two distinct clades (Fig. 5A). Furthermore, two distinct lineages were revealed within the KASI clade. One clade, which we have named KASI-1, contained both Orchidaceae and all non-Orchidaceae KASI sequences. The other clade was an orchid-specific lineage (11 species) showing evidence for two rounds of duplication. Thus, we have named these clades KASI-2A and KASI-2B, according to their phylogenetic placement. Both Chiloglottis seminuda and C. trapeziformis, as well as Caladenia plicata, possess all three KASI duplicates. In other orchids, only two KASI copies were found in A. wallichii, N. zollingeri, O. sphegodes and P. bellina while three were found in O. italica and Cypripediumformosanum (Fig. 5B). Conversely, nearly all orchids analysed contain one KASII copy, with no distinct splits observed within that clade.

Graph: Fig. 5. The evolution of KETOACYL-ACP SYNTHASE (KAS) in the Orchidaceae. (A) A maximum-likelihood tree of KASI and KASII protein sequences including 13 Orchidaceae and five non-Orchidaceae species. Purple and blue boxes indicate plant KASI and KASII clades, respectively. Red, green and blue bars depict KASI-1, KASI-2A and KASI-2B clades, of which the former two are Orchidaceae-specific. Dashed rectangles indicate the KASI-1 (blue), KASI-2A (red), and KASI-2B (green) branch used for selection analysis. Bootstrap support values for each branch are depicted as coloured circles. All sequence information used for phylogeny construction is available from Ensembl Plants and Orchidstra 2.0. (B) Distribution of KASI sequences in the 13 Orchidaceae species assayed. (C) Boxplot illustrating the ω (from CodeML free-ratio model) distribution of branches in Orchidaceae KASI-1, KASI-2A and KASI-2B clades. Outlier ω values (without either non-synonymous or synonymous sites) were removed from statistical testing. Asterisks indicate significant differences between observed ω values at P < 0.05 (*) based on Welch's t-test. n.s. not significant. (D) ω values of selected KASI-2A and KASI-2B branches. Black upward triangle and open downward triangle indicate positive and relaxed purifying selection, respectively.

In light of these new phylogenetic results, we propose a revision of the existing KASI nomenclature used in our earlier studies ([59]). As such, the previous C. trapeziformis KASI-1, KASI-2 and KASI-3 annotation is synonymous with the KASI-2A, KASI-2B and KASI-1 sequences discussed in this work, respectively.

Testing for signatures of selection across the three Orchidaceae KASI clades revealed evidence for strong purifying selection within the KASI-1 and KASI-2A clades, whereby the median observed ω of branches within them were 0.1 and 0.12, respectively (Fig. 5C). By comparison, the KASI-2B clade showed evidence for relaxation of purifying selection, with the observed ω ratios generally larger compared to those observed in the KASI-1 and KASI-2A clades (P < 0.05). Furthermore, we note: (1) a signal indicative of strong positive selection (ω7) at the split of KASI-2B, suggesting strong selection for a new function; (2) evidence for relaxed purifying selection (ω5) in the Chiloglottis KASI-2B clade and strong purifying selection (ω1) in the Chiloglottis KASI-2A clade; (3) no indication of positive or relaxed purifying selection was found at the split (or after) of KASI-2A (ω1–ω4); and (4) in Caladenia plicata, an Australian non-chiloglottone-accumulating species, purifying selection was indicated for both KASI-2A (ω1) and KASI-2B (ω6) sequences (Fig. 5D).

By applying a combination of transcriptomics, comparative genomics and molecular inhibition experiments, we have addressed three specific questions with the goal to progress our understanding of the biosynthesis of chiloglottones. Beginning with the construction of a comprehensive (three tissue pools: CAL, GS, SP) and reliable floral transcriptome of C. seminuda, we obtained a final assembly containing 28 778 transcripts with key assembly statistics (e.g. N50 score) closely mirroring other previously assembled orchid transcriptomes within the Orchidstra 2.0 database ([9]) and C. trapeziformis ([60]).

Comparison of tissue-specific floral transcriptomes revealed that the active CAL and GS tissues are both highly specialized (Fig. 1B) sharing a common theme of enriched lipid metabolism pathway genes (Fig. 1D, E). Such observations are congruent with our previous studies showing the tissue-specific expression of FA pathway genes in the active chiloglottone 1-emitting callus compared to the non-active labellum of C. trapeziformis ([59]). Despite a consistent trend of differential expression of both biosynthetic and degradative FA pathway genes compared to non-active tissues, the relative magnitude of expression did vary between active CAL and GS tissues of C. seminuda (Fig. 2A, B). In part, this is probably due to our inability to separate the sepal tissue from the active glandular structures on the surface.

Despite reaffirming the general involvement of FA pathway genes in chiloglottone formation across Chiloglottis orchids ([4]), there are differences in transcriptional trends between the active tissues of C. seminuda and C. trapeziformis. The extent of coordinated regulation of the entire pathways is lower in C. seminuda compared to C. trapeziformis ([59]). Only one isoform of KAR, ACX and AIM1-1 was induced in the CAL and GS of C. seminuda while multiple isoforms (four to six) were induced in C. trapeziformis calli. In addition, KASIII and KAT were not differentially expressed in any C. seminuda tissues while the latter two transcripts were induced in C. trapeziformis calli. Therefore, these genes may not be directly involved in chiloglottone formation.

Further supporting evidence for the role of FA biosynthesis as the major biosynthetic route for chiloglottones has been recently provided by experiments using inhibitors of KASI ([59]). Using cerulenin, one of several known compounds that irreversibly inhibit KASI ([42]; [33]), we demonstrated that chiloglottone production in the callus was significantly reduced upon induction with UV-B compared to controls in three Chiloglottis species, namely C. trapeziformis, C. valida and C. aff. valida.

The KASI (KASI-1, KASI-2A and KASI-2B) and KASII sequences of C. seminuda all contained the highly conserved catalytic cysteine–histidine–histidine triad (active site) common to KAS enzymes, thus rendering them potentially susceptible to fatty acid synthase inhibitors such as cerulenin (Supplementary Data Fig. S2). Treatment with cerulenin at 100 µm (as in [59]) resulted in significant inhibition (36 %) in the GS, while at higher concentrations (500 µm), chiloglottone 1 production was inhibited by 57 % and 49 % in the CAL and GS, respectively (Fig. 3). Thus, we confirmed that cerulenin had the same effect on chiloglottone inhibition in C. seminuda although the level of inhibition was lower and was tissue-dependent.

Comparative genomics using predicted single-copy orthologues is becoming a feasible approach for phylogenetic construction within a specific genus or across multiple orchid genera ([12]; [70]; [25]). This phylogenetic tree can then be used to guide the identification of genes that evolved under positive selection (ω or dN/dS > 1) at a transcriptome-wide scale ([53]). In this study, an ML tree was constructed using a concatenated alignment matrix of 201 putative single-copy genes shared across 25 species (after stringent filtering) to infer the positioning of Chiloglottis orchids in a wider phylogenetic context (Fig. 4). The constructed ML tree was largely congruent with previous Orchidaceae molecular phylogeny ([12]; [23]; [70]), placing Apostasioideae as sister to all orchids, followed by the Cypripedioideae as sister to the Orchidoideae and Epidendroideae.

Positive and/or relaxed purifying selection acting on just one or a few amino acid changes is predicted to drive the evolution of new substrate specificity/preference. In turn, this enables the evolution of novel reaction products such as new classes of plant volatile-producing enzymes ([54]; [29]; [26]; [71]). In the orchid genus Gymnadenia, relaxed purifying selection acting on eugenol synthase genes appear to underpin the evolutionary switch from the enzymatic production of eugenol (in G. odoratissima and G. conopsea) to eugenol and isoeugenol (in G. densiflora) ([26]). Evidence for strong positive selection on SAD2 genes in some Ophrys orchids is linked to the novel production of (Z)-9 and (Z)-12 alkenes ([54]; [65]). In the orchid genus Cypripedium there is evidence for positive selection at lipid metabolism pathway genes (e.g. CAC1, FAB2) across multiple species and such events may have been associated with the evolution of Cypripedium floral scent differentiation ([25]).

By using genome-wide selection analysis ([53]), we have identified two lipid metabolism pathways genes, FATA and KASI-2B, that are candidates for selection for new gene function in Chiloglottis (Table 1). To date, many plant FATA enzymes have been characterized (mangosteen, Arabidopsis, soybean, sunflower, etc), some of which are known to preferentially act on long-chain unsaturated acyl-ACPs, i.e. 18:1 ACP over saturated ones, e.g. 18:0 and 16:0 ACP ([32]). Notwithstanding the lack of differential expression, changes in critical residues can also alter substrate specificity of FATA and alter acyl-CoA pools ([15]; [43]), and thus FATA remains a candidate gene for precursor supply for chiloglottone biosynthesis. KASI-2B is also a prime candidate showing evidence of relaxed purifying selection and strong tissue-specific induction in the active tissues.

Phylogenetic analysis of plant KAS enzymes revealed that the Orchidaceae-specific KASI-2 lineage diverged before the monocot–eudicot split in the KASI-1 lineage (Fig. 5A). We hypothesize that an ancient KASI duplication once common to angiosperms was retained in the Orchidaceae but was lost in other plant lineages (e.g. eudicots and Poaceae) during speciation. This is not unexpected as hundreds of gene families are estimated to have undergone differential contraction, expansion and retention across plants and the Orchidaceae ([6]; [70]). Nonetheless, this raises the question regarding the role of these potentially retained KASI-2 duplicates in Orchidaceae.

While most Apostasioideae and non-Apostasioideae orchids assayed had KASI-2A sequences, fewer orchid lineages had KASI-2B sequences (Fig. 5B). Therefore, we hypothesized that duplication of KASI genes in select lineages by mechanisms other than WGD (e.g. tandem and segmental duplication) is likely to underpin the expansion of the KASI-2B clade as the relevant (and closely related) species of this clade are often diploid and have little variation in chromosome numbers (2n = 38–42) ([47]; [3]; [39]). In Phalaenopsis orchids, one recent study showed that segmental and tandem duplications contribute to the expansion of the chalcone synthase gene family ([36]9).

In light of the selection analysis outcomes (Fig. 5C, D), we hypothesize that the Orchidaceae-specific KASI-2A duplicates may have conserved functions common to plant KASI-1 while Chiloglottis KASI-2B genes may have evolved novel function along a two-tiered evolutionary path – strong positive selection for a new function shared with all other orchids within the clade followed by relaxed purifying selection unique to Chiloglottis. These events may have been the determining factor for enabling the predicted novel enzymatic functions probably required for chiloglottone biosynthesis (Fig. 6). For example, we speculate that Chiloglottis KASI-2B genes may possess short chain length acyl-ACP substrate specificities as seen with KASI duplicates in FA-accumulating plants ([67]) and activated short chain FA precursor condensation activities found in certain 2,5-dialkylcyclohexane-1,3-diones-accumulating bacteria ([19]).

Graph: Fig. 6. Predicted biosynthetic steps for chiloglottone 1 biosynthesis. Blue and green compounds indicate activated precursors. X indicates coenzyme A or acyl-carrier-protein. As an example, condensation of 3-ketohexanyl-ACP and 2-hexenyl-CoA to form 2-ethyl-5-propylcyclohexan-1,3-dion-4-carboxylate, the penultimate precursor to chiloglottone 1, may involve the biochemical activities of Chiloglottis KASI-2B proteins.

While morphology has now been shown experimentally to play an important role in the deceptive pollination of Chiloglottis orchids ([30], [31]9), nonetheless, chiloglottones hold the key for long-distance pollinator attraction as well as pollinator discrimination in the genus ([49]). The biosynthesis of chiloglottones has been predicted to involve FA metabolism pathway genes ([18]; [4]). In a previous study with C. trapeziformis, differential expression analyses of developmental stage transitions (e.g. very young buds vs. very mature buds vs. flowers) and specific tissues (calli vs. labellum) implicated a small subset of four FA pathway genes – CtKASI-2B, CtFATB2, CtACX2/3 and CtACX4 – as strong candidates for involvement in chiloglottone biosynthesis ([59]).

Here, we demonstrate that tissue-specific differential expression of many FA biosynthesis and β-oxidation pathway genes are also evident in the flowers of C. seminuda. These results demonstrate that such transcriptional trends are conserved across Chiloglottis orchids, thus further implicating these genes in chiloglottone biosynthesis. We further demonstrate that cerulenin significantly inhibits chiloglottone production in C. seminuda, adding to the growing evidence that chiloglottone production is likely to involve FA biosynthetic enzymes ([59]).

As a complementary approach to differential expression analysis, we used selection analysis to determine whether FA pathway genes are under positive or relaxed purifying selection in Chiloglottis orchids. We found that the majority of FA biosynthesis and β-oxidation pathway genes are under strong purifying selection, including those that show strong and highly specific induction in the active tissues (e.g. FATB2, KASI-1, KASII and LACS1). By contrast, Chiloglottis KASI-2B was the exception – showing evidence for relaxed purifying selection, a condition potentially enabling the evolution of a new function, as well as having strong tissue-specific induction in active tissues and across species.

Therefore, our combined analysis places KASI-related genes, especially KASI-2B of C. seminuda and C. trapeziformis, as prime candidates for chiloglottone 1 formation, and thus warrant further functional analysis. These findings also reinforce the prediction that the crucial 'starting point' for the formation of chiloglottones in Chiloglottis species involves FA biosynthesis and degradative pathways. We predict that our findings will have wider implications for other orchid systems such Arthrochilus and Paracaleana that are also known to produce chiloglottones ([49]). Our approach, which combines strategic transcriptome, comparative genomics, phylogenetic and selection analysis, is widely applicable for prioritizing candidate genes in any non-model plant system, thus making downstream biochemical investigations both feasible and rewarding ([61]; [59]9).

Supplementary data are available online at https://academic.oup.com/aob and consist of the following. Table S1: Transcript sequences (in FASTA format) of the final representative C. seminuda floral transcriptome assembly. Table S2: Comprehensive annotation for the representative C. seminuda floral transcriptome assembly. Table S3: Differential expression statistics generated by DESeq2 for treatments TS1 (CAL vs. SP), TS2 (GS vs. SP) and TS3 (CAL vs. GS). Only differentially expressed transcripts (|log2FC| > 0.5, FDR < 0.05) are listed. Table S4: Detailed enrichment results for MapMan BIN and GO categories for upregulated (pos), downregulated (neg), and intersect of pos and neg groups in TS1 and TS2. Table S5: List of all genes under positive selection (PSG) in Chiloglottis orchids (Chiloglottis branch) and their associated selection pressures (ω or dN/dS). Selection pressures for a subset of fatty acid pathway transcripts (FA) are also shown. Fig. S1: Distribution of transcripts in respective MapMan BIN and GO functional categories. Fig. S2: Alignment of C. seminuda and C. trapeziformis KASI (KASI-1/-2A/-2B) and KASII sequences. Red arrows indicate the cysteine–histidine–histidine triad (active site) typical of KAS enzymes.

We thank the Australian National Botanic Gardens, Canberra, for permission and permits to take plant material from the gardens and the New South Wales (NSW) National Parks & Wildlife Service for a scientific licence to take plant material from NSW. The photographs of Chiloglottis orchids were taken by Rod Peakall. This work was supported by Australian Research Council projects DP1094453 and DP150102762 to R.P. and E.P.