Genes involved in sex pheromone biosynthesis of Ephestia cautella, an important food storage pest, are determined by transcriptome sequencing

Abstract

Background

Insects use pheromones, chemical signals that underlie all animal behaviors, for communication and for attracting mates. Synthetic pheromones are widely used in pest control strategies because they are environmentally safe. The production of insect pheromones in transgenic plants, which could be more economical and effective in producing isomerically pure compounds, has recently been successfully demonstrated. This research requires information regarding the pheromone biosynthetic pathways and the characterization of pheromone biosynthetic enzymes (PBEs). We used Illumina sequencing to characterize the pheromone gland (PG) transcriptome of the Pyralid moth, Ephestia cautella, a destructive storage pest, to reveal putative candidate genes involved in pheromone biosynthesis, release, transport and degradation.

Results

We isolated the E. cautella pheromone compound as (Z,E)-9,12-tetradecadienyl acetate, and the major pheromone precursors 16:acyl, 14:acyl, E14-16:acyl, E12-14:acyl and Z9,E12-14:acyl. Based on the abundance of precursors, two possible pheromone biosynthetic pathways are proposed. Both pathways initiate from C16:acyl-CoA, with one involving ∆14 and ∆9 desaturation to generate Z9,E12-14:acyl, and the other involving the chain shortening of C16:acyl-CoA to C14:acyl-CoA, followed by ∆12 and ∆9 desaturation to generate Z9,E12-14:acyl-CoA. Then, a final reduction and acetylation generates Z9,E12-14:OAc. Illumina sequencing yielded 83,792 transcripts, and we obtained a PG transcriptome of ~49.5 Mb. A total of 191 PBE transcripts, which included pheromone biosynthesis activating neuropeptides, fatty acid transport proteins, acetyl-CoA carboxylases, fatty acid synthases, desaturases, β-oxidation enzymes, fatty acyl-CoA reductases (FARs) and fatty acetyltransferases (FATs), were selected from the dataset. A comparison of the E. cautella transcriptome data with three other Lepidoptera PG datasets revealed that 45 % of the sequences were shared. Phylogenetic trees were constructed for desaturases, FARs and FATs, and transcripts that clustered with the ∆14, ∆12 and ∆9 desaturases, PG-specific FARs and potential candidate FATs, respectively, were identified. Transcripts encoding putative pheromone degrading enzymes, and candidate pheromone carrier and receptor proteins expressed in the E. cautella PG, were also identified.

Conclusions

Our study provides important background information on the enzymes involved in pheromone biosynthesis. This information will be useful for the in vitro production of E. cautella sex pheromones and may provide potential targets for disrupting the pheromone-based communication system of E. cautella to prevent infestations.

Background

Pheromone-based methods of insect control are essential components of integrated pest management practices worldwide. The pheromones of over 2,000 insect species are now known, and The Pherobase is an updated compilation of pheromones and other behavior-modifying chemicals found in insects [1]. Common biosynthetic pathways have also been well-cited in many scientific publications over the last two decades, leading to production of species-specific pheromone compounds [2–5]. The female pheromones of almost all moth species are multicomponent blends of long hydrocarbon chains (10 to 18 carbons long), unbranched alcohols, and acetates or aldehydes, and are synthesized in the modified epidermal cells (pheromone-producing cells) from C16 or C18 fatty acid precursors [4, 6, 7]. A typical moth pheromone biosynthetic pathway begins even before the adult eclosion by releasing pheromone biosynthesis activating neuropeptide (PBAN) from the brain and transporting it to the pheromone gland (PG), which in turn activates functional group modification enzymes [3, 4, 8, 9] or acetyl-coenzyme A (CoA) carboxylase (ACC) [10]. As the first step in pheromone biosynthesis, carboxylation of acetyl-CoA to malonyl-CoA is catalyzed by ACC [10]. This is followed by fatty acid synthase (FAS) activity to produce saturated fatty acids (C18:0 and C16:0) using malonyl-CoA as the substrate. Later, the fatty acyl desaturases (DESs) introduce double bonds in the acyl chains, and then, specific β-oxidation enzymes shorten the chains. Once specific unsaturated pheromone precursors are formed, the terminal carboxyl group is modified to form one of the functional groups, alcohol, aldehyde or acetate ester (OH, CHO or OAc, respectively), and is catalyzed by fatty acyl reductase (FAR), aldehyde reductase (AR) or fatty acetyltransferase (FAT), respectively [3–5, 10]. A variety of desaturases, which introduce double bonds into the acyl at the ∆6 [11], ∆9 [12–14], ∆10 [15], ∆11 [13, 16, 17] and ∆14 [18] positions, have been cloned and functionally expressed from many moth species [2–5, 10]. Great progress has also been made in the functional characterization of FARs since their discovery in Bombyx mori [19] through detailed studies of pheromone evolution and the FARs of nine Ostrinia spp. [20–22], Yponomeuta spp. [23], Helicoverpa spp. and Heliothis spp. [24]. However, the molecular characterizations of other critical enzymes in the pheromone biosynthetic pathway, such as ACC, FAS, and several β-oxidation and acetylation enzymes, have not been characterized at the enzymatic level in insects.

Female moths typically start releasing sex pheromones a few days after emergence. In male moths, the pheromone molecule binds to odorant receptor (OR) proteins (in the antenna), signals are transmitted to the central nervous system where they are processed and identified by the brain, messages are then passed to the effector neurons, and finally the behavioral response is elicited. The expression of OR proteins is necessary and sufficient for odor detection in insects [25]. At first, volatile odors are bound to odorant-binding proteins (OBPs), a family that includes two sub-families, the pheromone-binding proteins (PBPs) and the general odorant-binding proteins (GOBPs) [26, 27]. Other important soluble secreted proteins that are found within the sensillum lymph include chemosensory proteins (CSPs) and the antennal binding protein X (ABPX) [28]. Finally, odorant molecules bind with ORs located in the dendritic membrane of receptor neurons [27, 29]. Sensory neuron membrane proteins (SNMPs) are another class of proteins involved in pheromone reception at the olfactory receptor neuron (ORN) [29–31]. Later, the signal termination is accomplished by the odorant-degrading enzymes (ODEs, also known as pheromone-degrading enzymes) [26, 32]. Knowledge of the olfactory communication system at the molecular level in insects is still in its early stages.

The tropical warehouse moth (almond moth), Ephestia cautella (Lepidoptera: Pyralidae) is a destructive polyphagous storage pest of wheat flour, dried figs, dates, nuts, chocolate, dried fruits, grain and associated processed food products worldwide. The control of these pests has depended exclusively on methyl bromide; however, methyl bromide was reported as facing an international phase-out by the year 2015 [33]. In this context, pheromones hold great potential in insect pest management [34]. In the last few decades, the elucidation of pheromone biosynthetic pathways, and the molecular characterization and functional gene expression of pheromone biosynthesis enzymes (PBEs) and OR proteins increased greatly [2, 10]. Most recently, through a synthetic biology approach, transgenic Nicotiana benthamiana plants with insect desaturases, FARs and FATs produced pure multi-component pheromone compounds [35]. Such in vitro production technology (green technology) could be cost effective and produce isomerically pure compounds that should be identical to chemically synthesized compounds. This research requires complete knowledge of the specific pheromone biosynthetic pathway and the functional characterization of enzymes (genes) involved in pheromone biosynthesis.

The rapid progress over the last decade resulted from the convergence of modern techniques from different areas of science has enriched our knowledge of the genetics of pheromone-based communications and olfactory communication systems. Transcriptome sequencing strategies are efficient for identifying a large number of expressed genes in specific tissues; thereby, providing information on the physiological, as well as molecular, properties of the tissue. Over the last few years, next-generation sequencing (NGS) techniques have provided easy and effective methods for the discovery of novel genes. These approaches are particularly relevant when no genomic data are available for the target species [36]. Over the past 5 years, RNA sequencing data on the insect pheromone gland amassed rapidly [37–41]. In the present study, using the Illumina sequencing approach, we constructed the transcriptome dataset of the PG of E. cautella and identified genes with putative roles in pheromone biosynthesis, transport and degradation. We combined the transcriptomic datasets with the E. cautella female sex pheromone precursors characterized through GC-MS studies, identified specifically or highly abundantly expressed genes in the PG and proposed roles for them in pheromone biosynthesis, binding, transport and release.

Results and discussion

Sex PG extraction and fatty-acyl precursor analysis

Analysis of the E. cautella PGs excised at the calling period (2-day-old, at mid-scotophase) showed the presence of the compound (Z,E)-9,12-tetradecadienyl acetate (Z9,E12-14:OAc) by their GC retention times (18.19 min) and mass spectra [ion fragment of m/z 61, a characteristic of acetate compounds (CH3COOH2+) and diagnostic ion at m/z 192] in comparison with those of authentic pheromone samples (Fig. 1). Our results were consistent with the earlier reports of the E. cautella female sex pheromone [34, 42]. Many studies have reported geographical variations and host-induced changes in the sex pheromone compounds and pheromone blend ratios in moths [2, 20, 21]. We isolated E. cautella (dried date fruit strain) sex pheromones to identify such differences. Date Palm (Phoenix dactylifera L) has been cultivated in Middle Eastern countries since ancient times, and E. cautella is native to Saudi Arabia where it infests dried date fruits in storage houses. To determine sex pheromone differences in the native moth strain, we studied its pheromone biosynthetic pathway as follows.

Fig. 1

Pheromone compound analysis of PG extract from E. cautella by GC-MS

Fatty acid methyl esters (FAMEs) were made from the total lipid extract of E. cautella PG to determine the corresponding fatty acid precursors of Z9,E12-14:OAc. The PG extracts contained unsaturated and saturated FAMEs, such as methyl hexadecanoate (16:COOMe, related abbreviations used hereafter for similar FAMEs), 14:COOMe, 15:COOMe; Z9-16:COOMe; E9-16:COOMe; Z11-16:COOMe; E11-16:COOMe; 17:COOMe; 18:COOMe; Z9-18:COOMe and Z9,Z12-18:COOMe. The FAMEs were identified on the basis of their retention times relative to those of the authentic standards, as well as on their mass spectra [characterized by parent ions and an intense m/z = 74 and 87, M+, M+ −31, M+ −32, M+ −74 (example, C16: Me = 270, 239, 238, 196 and C14: Me = 242, 211, 210, 168, respectively)]. The GC-MS analysis of a methanolyzed gland lipid extracts showed the corresponding precursors, C16:acid, C14:acid, E14-16:acid, E12-14:acid and Z9,E12-14:acid (Fig. 2a). The mono- and di-unsaturated precursor GC retention times and mass spectra matched those of the authentic standard samples (Fig. 2b). When comparing the relative proportions of the derived acids, E12-14:acid appeared to be more abundant. There were also large amounts of other FAMEs identified as those of E9-16:COOMe and Z9-16:COOMe, as well as small amounts of others tentatively assigned as E11-16:COOMe and Z11-16:COOMe (Fig. 2a) (diagnostic ions at m/z 242 for 14:COOMe, m/z 240 for E12-14:COOMe, m/z 270 for 16:COOMe, m/z 268 for Z9-, E9-, E11- and Z11-16:COOMe, m/z 252 [M+ −32] and 284 [M+] for E14-16:COOMe and m/z 206 [M+ −32] and m/z 238 [M+] for Z9,E12-14:COOMe).

Fig. 2

Fatty acid analysis of PG extracts from E. cautella by GC-MS. The fatty acids methyl ester (FAME) were identified based on the retention time (RT) of the authentic standard compound and mass spectra analysis

Based on the identified pheromone precursors, the putative sex pheromone biosynthetic pathway of E. cautella was predicted as shown in Fig. 3. In E. cautella five major pheromone precursors, C14:acid, C16:acid, E14-16:acid; E12-14:acid and Z9,E12-14:acid, were identified using the FAME analysis of the PG. Thus, it is rational to propose a pheromone biosynthetic pathway in which the saturated fatty acid precursor of the E. cautella sex pheromones is palmitic acid (16:0) that is desaturated by Δ14-desaturase to form the pheromone precursor E14-16:acyl-CoA, which in turn has its chain shortened by β-oxidation to E12-14:acyl-CoA. A unique Δ9-desaturase uses the E12-14:acyl-CoA to produce Z9,E12-14:acyl-CoA that is reduced and acetylated to form Z9,E12-14:OAc, the final pheromone compound (Fig. 3). An alternative pathway can also be proposed that involves the chain shortening of C16:acyl-CoA to C14:acyl-CoA, which is later desaturated by Δ12-desaturase to produce E12-14:acyl-CoA, and then a unique Δ9-desaturase uses the E12-14:acyl-CoA to produce Z9,E12-14:acyl-CoA (Fig. 3). This is reduced and acetylated to form Z9,E12-14:OAc, the final pheromone compound of E. cautella. Based on the FAME analysis, the first pathway appears more fitting; however, further studies using in vivo labelling are required to test the hypothesis. We compared the E. cautella pheromone biosynthetic pathway with two Spodoptera spp., S. exigua and S. littoralis, that use Z9,E12-14:OAc as a sex pheromone compound (see Additional file 1: Figure S1). In the present study, we isolated the sex pheromone, Z9,E12-14:OAc from E. cautella infesting dried date fruit and identified a major pheromone precursor E12-14:acid. Hence, the proposed pheromone biosynthetic pathway (Fig. 3) appears to be more appropriate. The common biosynthetic pathway leading to the production of a moth sex pheromone compound, based on the activity of a desaturase with a strict regio- and stereo-selectivity, produced different pheromone precursors, which are characteristic of different species (Additional file 1: Figure S1) [3, 11–16].

Fig. 3

Proposed pheromone biosynthetic pathway leading to the sex pheromone of E. cautella, (Z,E)-9,12-tetradecadienyl acetate

Illumina sequencing and de novo assembly

Illumina sequencing of a cDNA library prepared from mRNA of the E. cautella PG produced 237,048,152 raw reads with an average length of 101 base pairs (bp). After trimming adaptor sequences and eliminating low quality reads, there were 231,851,937 reads (227,994,544 sequences in pairs and 3,857,393 single sequences) with an average length of 100 bp (Additional file 2: Table S2). The raw reads were deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database with the accession number SRX646348. After assembly, with scaffolding, 83,792 transcripts with an average length of 590 bp were obtained, having a maximum length of 19,518 bp. Most transcripts had lengths that ranged from 376 to 760 bp. The whole transcriptome size was 49.5 Mb, and the N50 size was 760 bp, with 10,856 sequences longer than 1 kb. This Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession GBXH00000000. In comparison with previously reported PG transcriptome/EST data [37–41], this pooled assembly of E. cautella PG sequences has the second greatest data volume and sequence lengths (Additional file 2: Table S2).

Functional annotation

The assembled transcripts were used as query in a BLASTx against the non-redundant (nr) NCBI protein database, UniProtKB, Flybase and KEGG, all with an e-value cut-off of 10E − 5. Most of the sequences had an e-value between 1.0E − 4 and 1.0E − 10 (Additional file 3: Figure S3A). The similarity between E. cautella PG sequences and those of the databases ranged from 36 % to ~100 % (value: 3,434) with a peak at 65 % (value: 12,152) (Additional file 3: Figure S3B). A BLAST2GO analysis of the 83,792 transcripts of the E. cautella PG resulted in 30,582 transcripts with blast hits, 53,210 without blast hits, 4,217 with mapping results and 20,615 annotated sequences (Additional file 4: Figure S4A). The sequences without blast hits may have low similarities to functionally similar genes in the database, novel genes or parts of the 5′ or 3′ UTR regions. The PG transcript of E. cautella produced the most significant hits to B. mori sequences, followed by Danaus plexippus sequences (Additional file 3: Figure S3C). The evidence code distribution for the BLAST hit chart indicates an over-representation of Inferred Electronic Annotation (IEA), followed by Inferred by Mutant Phenotype (IMP) and Inferred by Direct Assays (IDAs) (Additional file 5: Figure S5A). The maximum evidence code for the individual sequences was through IEA, IMP and lastly IDA (Additional file 5: Figure S5B). The majority of functional predictions from the coding sequences were obtained from UniProtKB followed by FlyBase (FB) (1,017,318 and 107,608, respectively) (Additional file 5: Figure S5C).

GO terms were assigned by BLAST2GO through a search of the nr database, and INTERPRO was searched using INTERPROSCAN, resulting in ~34,953 transcripts from INTERPRO and 48,838 transcripts ‘without INTERPRO’ that had GO-annotation average lengths of 590 bp. Using this method, 12,455 unigenes were assigned to one or more GO terms. ANNEX was run after BLAST, and INTERPROSCAN results were annotated with the following results: 105,242 total original annotations, 7,810 new annotations, 1,007 original annotations replaced by new annotations due to specificity, and 3,853 confirmed annotations.

As shown in Table 1, 30,097 genes, 35 % of all transcripts in nr, 14,036 genes in UniProtKB, and 20,615 enzymes encoded in the Kyoto Encyclopedia of Genes and Genomes (KEGG) returned cut-off blast hits > 1.0E − 5. A KEGG metabolic pathway analysis revealed 5,762 transcripts could be assigned to generate 130 predicted pathways (Additional file 6). The major enzyme commission (EC) classes included oxidoreductases (964 transcripts), transferases (2,576 transcripts), hydrolases (2,125 transcripts), lyases (200 transcripts), isomerases (125 transcripts) and ligases (391 transcripts). The KEGG pathway map revealed the presence of a large number of PG transcripts involved in fatty acid biosynthesis (42 transcripts, 8 enzymes), fatty acid elongation (34 transcripts, 6 enzymes), fatty acid degradation (97 transcripts, 15 enzymes) and most importantly, the biosynthesis of unsaturated fatty acids (Additional file 7: Figure S6) and genes (enzymes) that may participate in pheromone biosynthesis (pathway: see Fig. 3). In total, 65 transcripts encoding five enzymes, DES, FARs, FATs, Acyl-CoA oxidases and dehydrogenases, were assigned functional annotations (Additional file 7: Figure S6).

Table 1 Annotation of a pooled assembly, representing the E. cautella PG transcriptome

GO for the genes expressed in the E. cautella PG

Based on the matches to INTERPRO proteins, the E. cautella PG transcriptome was GO-annotated. The annotation results and distribution, GO-level distribution, number of GO-terms for E. cautella sequences with a specific length (x-axis), annotation score distribution and the percentage of E. cautella sequences with a specific length (x-axis), are depicted in Additional file 4: Figure S4. Among the total transcripts with BLAST results, 63 % were assigned GO terms, 32 % were unannotated proteins that had no matches in the GO database and 5 % were sequences assigned as predicted uncharacterized proteins (Fig. 4). The proteins with associated GO terms, such as “molecular function”, “biological process” and “cellular component” were grouped and recorded at different match levels (Fig. 4). The “cellular process” (10,559) and “metabolic process” (9,305) GO categories had the most abundant transcripts within the “biological process” GO ontology (Fig. 4). In the “cellular components” the most abundant transcripts were in “binding” (9,941) and “catalytic activity” (8,479) (Fig. 4). In the “cellular components” the transcripts were mainly distributed in “cell” (7,207) and “cell part” (7,123) (Fig. 5). In the “molecular function” ontology, 9,941 transcripts with “binding” functions were annotated, as were 8,479 that had “catalytic activity” (Fig. 5c). Of the proteins that had matches to the nr database, the most abundant protein class was the binding proteins. Other highly abundant proteins included oxidoreductase proteins, kinases, peptidases, cytoskeletal proteins, ribosomal proteins and proteins involved in other major functional categories (Fig. 5c). Of the direct GO counts identified for “biological process”, lipid metabolic process (including pheromone biosynthesis) and reproduction were among the first 20 dominant terms (Fig. 5a). Of the categories enriched for the direct GO counts identified as “cellular component”, the protein complex, nucleus and cytoplasm were the largest groups (Fig. 5b).

Fig. 4

Pie and Stack chart showing the percentage of E. cautella predicted genes as annotated proteins, predicted proteins and unannotated proteins

Fig. 5

Distribution of enriched functions in a) Biological Process (BP), b) Molecular Functions (MF) and c) Cellular Component (CC)

Transcript abundance in the E. cautella PG

The highly expressed transcripts in the E. cautella PG are summarized in Table 2. The most abundant transcripts included vitellogenin and vitellogenin precursor (Total read count: 1,646,398 and 148,469, respectively), a major reproductive protein and its precursor, respectively, in insect egg production [43]. The results were consistent with a previous report of transcriptionally abundant proteins in Agrotis ipsilon’s PG [39]. The acyl-CoA desaturases, contigs 349 and 1,286, were highly expressed in the PG with 3,656 and 4,967 reads per kilobase per million reads (RPKM), respectively, indicating their roles in pheromone biosynthesis. Other highly abundant transcripts were of contigs 88 and 106 with 7,575 and 3,711 RPKMs, respectively, encoding CSPs that exhibited a 61 % identity with Sesamia inferens (Genbank: AGY49267) [40] and 62 % with A. ipsilon (Genbank: AGR39573) [39], respectively. The major housekeeping genes, such as elongation factor, cytochrome c oxidase subunit I and III, and circadian clock-controlled protein (period gene), were highly expressed in the PG of E. cautella (Table 2).

Table 2 The most abundant mRNAs in the E. cautella PG

Comparative analysis of PG transcripts in Lepidoptera

By comparing E. cautella PG transcripts with those of B. mori and H. virescens from the NCBI database of differentially expressed transcripts and A. ipsilon from the SRA database, a large number of PG transcriptome sequences were found to be homologous. After assembly, we obtained 17,508 unigenes from A. ipsilon and 11,001 and 13,612 ESTs from B. mori and H. virescens, respectively.

We selected the first 10 bidirectional hits for each transcript from E. cautella, A. ipsilon, H. virescens and B. mori (producing a total of 240,753, 134,988, 46,912 and 46,940 blast hit results, respectively) for the comparative analysis. When comparing the PG transcripts pairwise using the bidirectional blast hit results, we found that between E. cautella and the three other Lepidoptera, 45 % of the blast hits were shared, and 65 % of the blast hits were unique to E. cautella (Fig. 6). The comparison between E. cautella and A. ipsilon showed that 20 % of the blast hits were shared, and 65 % of the blast hits were unique to E. cautella (Fig. 6). Similarly, a comparative analysis of blast hits of E. cautella, A. ipsilon and B. mori showed that 3.5 % of the blast results were shared. Comparative blast hits of E. cautella, A. ipsilon and H. virescens showed that 3.9 % had homologous hits, while between E. cautella and H. virescens there was 1.9 %, and between E. cautella and B. mori there was 1.4 % shared blast hits (Fig. 6). A large portion of the E. cautella transcripts (65 %) had no homologous hits in the available PG transcriptomes/ESTs of the other three species. This may have been because of the larger data set (83,792 transcripts) for E. cautella and the lower coverage in the other studies (Fig. 6). The high number of E. cautella blast hit results, which did not match A. ipsilon, B. mori or H. virescens may be due to novel genes with unique functions or highly conserved genes.

Fig. 6

Venn diagram showing the comparative analysis of the E. cautella PG transcriptome with those of A. ipsilon, B. mori and H. virescens

Identification of candidate genes involved in pheromone biosynthesis

In the present study, the E. cautella pheromone compound identified was Z9,E12-14:OAc, and the pheromone biosynthetic pathway is likely to be similar to those in other Pyralid moths (or type I pheromone biosynthesis), which include fatty acid synthesis (ATP-dependent carboxylation and decarboxylation condensation with several malonyl moieties), including the actions of desaturases and β-oxidation enzymes, followed by modifications of the carboxyl group by reductases and acetyltransferases [44]. Using BLASTx searches, we identified members of gene subfamilies in the E. cautella PG transcriptome putatively involved in Z9,E12-14:OAc pheromone production (Table 3). These include two PBAN receptor isoforms, five fatty acid transport proteins (FATPs), six ACCs, 12 FASs, 22 DESs, 28 FARs, 18 FATs and 11 ARs (Table 3). Additionally, 87 transcripts encoding putative β-oxidation enzymes, including 28 acyl-CoA dehydrogenases, 17 acyl-CoA oxidases, 13 enoyl-CoA hydratases, 17 L-3-hydroxyacyl-CoA dehydrogenases, eight 3-ketoacyl-CoA thiolases, three delta-3, delta-2 trans-enoyl-CoA isomerases and a delta(3,5)-delta(2,4)-dienoyl-CoA isomerase, were identified (Additional file 8: Table S7). There were also 36 transcripts encoding putative pheromone degrading enzymes (Additional file 9: Table S8), three transcripts encoding putative ABPs, 17 transcripts encoding putative OBPs, seven candidate CSPs, two transcripts encoding PBPs, 21 candidate ORs, two candidate sensory neuron membrane proteins and three candidate ionotropic receptors (IRs) (Additional file 10: Table S9 and Additional file 11: Table S10). Their abundance levels, based on RPKM values, in the PG transcriptome are shown in Table 3.

Table 3 Putative pheromone biosynthesis enzymes (PBEs) in the E. cautella PG

The PBAN receptor

Previous studies concluded that the sex pheromone biosynthetic machinery of Lepidopteran PG cells is regulated by PBAN, which is released from the brain, goes to the hemolymph and binds to the PBAN receptor in the membrane of pheromone producing cells, triggering pheromone production [8, 9]. We found two transcripts, EP_contig_27375 and EP_contig_24961, encoding proteins highly homologous to PBAN receptor isoforms C and A, respectively (Table 3). They have very low abundance levels in the E. cautella transcriptome (4.38 and 0.8 RPKM) (Table 3) but high identities (60–63 %) to the O. nubilalis PBAN receptors C and A in GenBank (AGL12068 and AGL12066, respectively) [45]. The PBAN receptors functionally characterized from O. nubilalis [45] and H. virescens [46] include isoforms A and C, and in the present study we identified PBAN isoforms A and C from E. cautella PG, which should be involved in pheromone production. The sequence identity (93 %) of E. cautella PBAN isoforms A and C indicate that they are likely produced by alternative splicing at the 3′-end of the receptor gene as reported in other moths, generating multiple receptor subtypes [46]. We also found a G-protein-coupled receptor (EP_contig_34693) that shows homology (63 %) to the diapause hormone receptors of H. zea and B. mori (AGR34305 and NP_001036913, respectively). The diapause hormone receptor is a G-protein gamma-subunit homolog, which is hypothesized to interact with the PBAN receptor, and has been reported in the PG transcriptomes of A. segatum [37] and H. virescens [38].

Fatty Acid Transport Protein (FATP) [EC:6.2.1.-]

FATPs belong to an evolutionarily conserved family of membrane-bound proteins that facilitate the uptake of extracellular long-chain fatty acids (LCFAs), and/or very LCFAs, and catalyze the ATP-dependent esterification of these fatty acids to their corresponding acyl-CoA derivatives [47]. The important role of FATPs in pheromonogenesis has been demonstrated in B. mori [47] and in O. scapulalis [48]. In E. cautella, we found five FATP isoforms in Unigene_3 (RPKM 140) with high transcript abundance levels and a high identity (80.8 %) to those of O. scapulalis (GenBank: BAJ33524) (Table 3).

Acetyl CoA Carboxylase (ACC) [EC:6.4.1.-]

Pheromone biosynthesis begins with an ACC catalyzing the production of malonyl-CoA from acetyl-CoA in the first committed biosynthesis step [49, 50]. In the E. cautella PG we found six transcripts encoding ACCs. ACC partial sequence EP_unigene_1, 2 and EP_contig_14940 showed more than 90 % identity with A. ipsilon ACC (GenBank: AGR49308). EP_unigene_3_ACC showed 93 % similarity with B. mori ACC (GenBank XP_004930758) (Table 3). Based on their RPKM values (286), EP_unigene_1, 2 and EP_contig_14940 were relatively highly expressed in the E. cautella PG (Table 3).

Fatty Acid Synthase (FAS) [EC:2.3.1.-]

In moth pheromone biosynthesis, FAS is supposed to catalyze the conversion of malonyl-CoA and NADPH to produce saturated fatty acids (16:acyl in E. cautella) [49]. We found 12 FAS-like partial transcripts in the E. cautella PG, which produced five different BLASTx hits in NCBI. Thus, we are proposing the existence of five FAS-like genes in E. cautella (Table 3). Partial sequences of EP_contig_284 and 8286 showed high similarity levels (<80 %) to A. ipsilon FAS (GenBank: AGR49310), whereas EP_contig_1101 showed a high similarity to D. plexippus FAS (GenBank: EHJ78836). The details of other FAS transcripts and BLASTx hit similarities are given in Table 3. Based on the RPKM value (110), EP_contig_1101 was highly expressed in the E. cautella PG (Table 3).

Desaturases (DES) [EC:1.14.19.-]

The desaturases introduce a double bond into the fatty acyl carbon chain, with strict regio- and stereo-selectivity. The desaturases characterized thus far include enzymes that act on saturated and monounsaturated substrates, which include Δ5 [51], ∆6 [11], Δ9 [12–14, 52, 53], Δ10 [15, 54], Δ11 [13, 16, 17, 55, 56] and Δ14 [18, 57]. Desaturases are characterized by having three histidine boxes containing eight histidine residues, which are used for binding essential metal complexes used in the enzyme reaction, and acyl-CoA desaturases introduce unsaturated bonds into fatty acids that are bound to CoA [58].

In E. cautella, four major pheromone precursors, C14:acid; E14-16:acid; E12-14:acid and Z9,E12-14:acid, were identified using a FAME analysis of the PG. In the E. cautella sex pheromones’ biosynthesis, a two-step desaturation process is proposed, involving ∆14, or ∆12, and ∆9 desaturases (Fig. 3). In the E. cautella PG transcriptome, 22 transcripts encoding desaturases have been identified (Table 3). EP_unigene_9 and EP_contig_349 are the highly expressed desaturases in the E. cautella PG (4,967 and 3,657 RPKMs, respectively), followed by EP_unigene_4, and EP_contigs_145 and 343 (3,391, 1,752 and 1,739 RPKM, respectively) (Table 3). Contigs_349 and 45 are closely related to the Ctenopseustis obliquana desaturase (GenBank: AER29852), which has ∆9, ∆11 and ∆14 fatty acid desaturase activities [12]. EP_unigene_9_1286 and EP_contig_343 showed high similarities to Amyelois transitella (GenBank: AGO96562) and O. furnacalis Z/E11 desaturases (GenBank: AAL32060), respectively, which have ∆11 and ∆9 desaturase activities [15]. EP_contig_349 and Contig_145 may have ∆14 and ∆9 desaturase activities (multifunctional) and could be involved in the formation of E14-16:acyl-CoA and Z9,E12-14:acyl-CoA, and EP_unigene_9_1286, or EP_contig_343, may have ∆9 desaturase activity and could be involved in Z9,E12-14:acyl-CoA synthesis from E12-14:acyl-CoA or multifunctional ∆12 and ∆9 desaturase activities (Fig. 3). Further studies on the functional gene expression levels of these desaturases in a transformed yeast (Saccharomyces cerevisiae) expression system are in progress.

The phylogenetic analysis of E. cautella desaturases with other moth desaturases is shown in the Fig. 7. Based on the phylogenetic tree, three possible candidate desaturases have been identified, EP_contig_349/EP_unigene_4, EP_unigene_9/EP_contig_343 and EP_contig_70932/EP_contig_71065, which form a clade with ∆9, ∆11, ∆12 and ∆14 desaturases (Fig. 7). EP_contig_349/EP_unigene_4 is closely related to O. nubilalis and O. furnacalis ∆11 and ∆9 desaturases (GenBank: AAL35331 and AAL320660, respectively). EP_unigene_9/EP_contig_343 is in the clade with Spodoptera littoralis desaturases, which have Z9 and E10,12 desaturase activities (GenBank: AAQ74259). The third putative desaturase type includes EP_contig_70932 and EP_contig_71065 and is closely related to the ∆14 desaturases of O. furnacalis and O. nubilalis (GenBank: AAL35746 and AAL35330, respectively) (Fig. 7). The one and only ∆14 desaturase reported so far is from a Pyraloidea moth, O. furnacalis [18, 57], although a later study showed several cryptic ∆11- and ∆14-desaturase genes exist in the O. nubilalis genome [59]. Further studies on the functional expression of desaturases in E. cautella will provide more insights into the origin and evolution of the ∆14 desaturases.

Fig. 7

Phylogeny of Lepidopteran desaturase genes

β-oxidation enzymes

Once the ∆14 desaturase introduces a double bond in palmitate it forms E14-16:acid (Fig. 3), which is later subjected to chain shortening by β-oxidation, resulting in the fatty acyl pheromone precursor, E12-14:acid (Fig. 3). In the alternative pathway, it is involved in the 16:acyl chain being shortened to 14:acyl by β-oxidation. β-oxidation is the action of a series of enzymes, working sequentially and forming a reaction spiral [35].

First, by the action of acyl CoA oxidase (ACO) (in peroxisomes) and acyl-CoA dehydrogenase (ACD) (in mitochondria), acyl-CoA is converted into E2-enoyl-CoA. There was an earlier report of four different ACDs, short-chain, medium-chain, long-chain and very-long-chain ACDs, depending on the fatty acyl chain-length specificities [60]. However, there is no report characterizing the ACDs involved in moth pheromone biosynthesis. It is possible that medium-chain ACDs could be more active because they act on hexanoyl-CoA, whereas long-chain ACDs preferentially act on octanoyl-CoA and longer chain-length substrates. We found many candidate genes of ACDs and acyl-CoA oxidases in the PG of E. cautella. In particular, EP_unigene_3_ACD [EC:1.1.1.211] and EP_unigene_2_ACD are the most abundant ACDs (271 and 219 RPKMs, respectively), and EP_unigene_1_ACO [EC:1.3.3.6] and EP_unigene_6_ACO are the most abundant acyl-CoA oxidases in the PG (192 and 127 RPKMs, respectively) (Additional file 8: Table S7). Moreover, we found two unigenes (EP_unigene_9_ACD [EC:1.3.3.6] and EP_contig_419_ACD) of isovaleryl coenzyme A dehydrogenase, which is specific to the metabolism of branched-chain fatty acids [61] (Additional file 8: Table S7).

The next step of β-oxidation involves E2-enoyl-CoA, which is reversibly hydrated by enoyl-CoA hydratase to L-3-hydroxylacyl-CoA. Two kinds of enoyl-CoA hydratases have been identified in mitochondria, one specialized for crotonyl-CoA (4C) and the other one being a long-chain enoyl-CoA hydratase, which effectively hydrates medium and long-chain substrates [62]. We found many candidate genes for enoyl-CoA hydratases and, among these, the EP_unigene_4_ECH [EC:4.2.1.17] is the most abundantly expressed in the PG of E. cautella (RPKM: 323). It shows a 93 % amino acid identity with the PG of Papilio xuthus (GenBank: BAM18079) (Additional file 8: Table S7).

The third reaction involves a reversible dehydrogenation of L-3-hydroxyacyl-CoA to 3-ketoacyl-CoA catalyzed by L-3-hydroxyacyl-CoA dehydrogenase. There are three different kinds of L-3-hydroxyacyl-CoA dehydrogenases that have been reported in mitochondria, long-chain, medium-chain and short-chain L-3-hydroxyacyl-CoA dehydrogenase (active with long-, medium- and short-chain substrates, respectively) [55, 62]. In the E. cautella PG, we found many candidate genes for long-, medium- and short-chain L-3-hydroxyacyl-CoA- dehydrogenases. EP_unigene_4_HCD [EC:1.1.1.35] is highly abundant in the PG (RPKM: 542), followed by EP-Unigene_2_HCD (RPKM: 271), and they both show high amino acid identities with D. plexippus and B. mori hydroxyacyl-CoA-dehydrogenases (GenBank: EHJ72407 and NP_001040132, respectively) (Additional file 8: Table S7).

Finally, 3-ketoacyl-CoA is cleaved by a thiolase between its α- and β-carbon atoms, producing the two carbon shorter substrate (E14-16:acid to E12-14:acid or C16: acid to C14: acid). Three kinds of thiolases exist in mitochondria, acetoacetyl-CoA thiolase or acetyl-CoA acetyltransferase (specific to acetoacetyl-CoA), 3-ketoacyl-CoA thiolase or acetyl-CoA acyl transferase (acts on C4-C16 unsaturated fatty acids), and long chain 3-ketoacyl-CoA thiolase (component enzyme of the membrane-bound tri-functional β-oxidation), where the first two kinds of thiolases are components of the soluble matrix enzyme complex [63]. We found eight candidate thiolase genes in the PG of E. cautella, with EP_contig_2704_KCT [EC:2.3.1.16] having the most highly abundant transcript in the PG (RPKM: 108) (Additional file 8: Table S7) and a 94 % amino acid identity with D. plexippus (Additional file 8: Table S7) (GenBank: EHJ7447).

The degradation of unsaturated fatty acids requires auxiliary enzymes, such as delta-3, delta-2 trans-enoyl-CoA isomerase and 2,4-dienoyl-CoA reductase, to modify the structure of the double bond during the β-oxidation process [64]. In the E. cautella PG, we found four delta-3, delta-2 trans-enoyl-CoA isomerases, two mitochondrial and two peroxisomal, and among these EP_unigene_1_TECI [EC:5.3.3.8] has the most abundant transcripts (RPKM: 177) (Additional file 8: Table S7). Additionally, we found a delta(3,5)-delta(2,4)-dienoyl-CoA isomerase [EC:5.3.3.-], which is specialized for processing odd-numbered double bonds [61] (Additional file 8: Table S7).

Moth pheromones generally consist of 10C-16C compounds synthesized from C16-C18 fatty acid moieties, involving many chain-shortening reactions [44]. Previous research has mostly been related to the desaturases and functional group modification enzymes, while research on the chain-shortening enzymes involved in pheromone biosynthesis has been meager. In the present study, we found many promising candidates that may be involved in β-oxidation, and further research on their heterologous expression, or RNAi, could reveal their significance in E. cautella pheromone biosynthesis.

Fatty Acyl Reductase (FAR) [EC:1.2.1.-]

FAR enzymes catalyze the reduction of fatty acyl precursors to fatty alcohols in a reaction that is dependent upon NADPH as a cofactor [19, 22–24]. FAR genes have been shown to function in pheromone biosynthesis in moth species directly through the production of an alcohol that confers species specificity or indirectly through the biosynthesis of precursor compounds [10, 21]. The number of FAR genes per genome can vary greatly between organisms. In vertebrates, there are two reductase genes present in the genomes, whereas there are more than a dozen present in the moth O. scapulalis [22]. The FAR gene family undergoes birth- and death-related evolution [65]. Even though the evolutionary origins of this gene family are not well understood, it has been assumed, based on protein sequence similarity, that the acyl-CoA synthetase, acyltransferase and oxidoreductase gene families are close relatives of this family, and thus form a superfamily [65]. In the E. cautella PG transcriptome pooled data, we identified 28 FAR-like genes, which included partial and full-length sequences (7), and BLASTx results identified them as putative FAR-like genes (Table 3). Based on the sequence assembly, multiple sequence alignment and BLASTx hit results, we named them uniquely; however, they may represent partial sequences of the same FARs. We took care to avoid duplications; however, if the full-length sequence was not available in our transcriptome dataset or in the NCBI database, then it could be present in partial sequences. EP_unigene _9 (2,019 bp) is the most highly expressed FAR in E. cautella’s PG (RPKM: 772). It showed a 67.6 % identity with D. plexippus (GenBank: EHJ72233), followed by EP_unigene _1 (RPKM: 742) and EP_unigene _14 (RPKM: 536). It shares a 77.3 % amino acid identity with O. nubilalis and a 59.1 % with D. plexippus (GenBank: ADI82778 and XP_004926010, respectively). All other FARs, except two, EP_unigene_2 and EP_unigene_3 (502 and 487 RPKMs, respectively), have a low abundance, with RPKM values of less than 100, in the PG transcriptome (Table 3). Further studies on the tissue specificity of each FAR to determine the PG-specific FARs is in progress.

The phylogenetic analysis of moth FARs is shown in Fig. 8. Based on the phylogenetic tree, three possible candidate FARs were identified, EP_contig_72742, EP_contig_45618 and EP_contig_79681, which formed a clade with moth pgFAR (Fig. 8). Until now pgFAR had been characterized from B. mori [19], O. scapulalis [22], nine Ostrinia spp. [20, 21], three Yponomeuta spp. [23], as well as Helicoverpa and Heliothis [24]. In the E. cautella PG, EP_contig_72742 and EP_contig_45618 FARs formed a cluster with the Ostrinia pgFARs (Fig. 8), and EP_contig_79681 formed a cluster with the Yponomeuta, Helicoverpa and Heliothis pgFAR clade (Fig. 8). Further studies on tissue-specific expression and heterologous gene expression in a yeast system are in progress.

Fig. 8

Phylogeny of Lepidopteran fatty acyl reductase (FAR) genes

Fig. 9

Maximum Likelihood (ML) tree of the Fatty Acetyltransferase (FAT) from the E. cautella PG and various FATs from the A. ipsilon and H. virescens PGs

Aldehyde Reductase (AR) [EC:1.1.1.-]

ARs are members of the aldo-keto reductase superfamily and can reduce long-chain acyl-CoA to form aldehyde intermediates [44]. In insects that use both an alcohol and an aldehyde as part of their pheromone, it is unclear how the production of both components occurs. Even though we did not identify any aldehyde precursors in the PG of E. cautella during the FAME analysis, we are still discussing ARs in this study because ARs and FARs share an evolutionary history, and the gene families are closely related [65]. In the E. cautella PG we identified 11 transcripts with homology to the aldo-keto reductases of A. ipsilon, Papilio xanthus, B. mori, Chilo suppressalis and D. plexippus (Table 3). The derived protein sequences of these 11 transcripts showed a 68–82 % amino acid identity with their homologs in other insects. All of the AR transcripts were present at a low abundance (less than 50 RPKM) in the PG transcriptome (Table 3); therefore, we assumed that AR does not have a role in E. cautella pheromone biosynthesis.

Fatty Acetyltransferases (FAT) [EC:2.3.1.-]

To produce the acetate ester pheromone components, most moths use an acetyl-CoA: fatty alcohol acetyltransferase that converts fatty alcohols to acetate esters [2–5]. The genes involved in this step have not been characterized from any insects [2–5, 10]. However, different acetyltransferases, which have a characteristic motif (HXXXD) and a conserved region (DFGWG), have been cloned from plants [66]. In the E .cautella PG transcriptome, 18 FAT-like genes were identified, showing homology to D. plexippus, B. mori, O. scapulalis and A. ipsilon (Table 3). Except for two FATs, they showed a greater than 70 % amino acid identity, the highest being EP_contig_7673 at 95 %, with A. ipsilon (GenBank: AGQ45625). We searched the conserved domains (CD) within the protein or coding nucleotide sequence databases at NCBI, but none of the E. cautella PG acetyltransferases belonged to the plant category FATs, suggesting that E. cautella may not express this gene family or that they have undergone substantial evolutionary changes. Nevertheless, most of the E. cautella FATs had hits to members of the N-acyltransferase (NAT) superfamily with CD accession no. cl17182 (Table 3). The CD accession nos. of E. cautella PG acetyltransferases are given in the Table 3 (fatty acetyltransferase). All FAT transcripts were present at low abundance levels (less than 120 RPKM) in the E. cautella PG transcriptome (Table 3); therefore, we could not predict which, if any, had roles in pheromone biosynthesis. However, phylogenetic analysis showed that EP_unigene_7 and EP_unigene_10 clustered with yeast (S. cerevisiae) alcohol acetyltransferase [EC:2.3.1.84] (GenBank: BAA05552.1), which catalyzes the esterification of isoamyl alcohol by acetyl coenzyme A (Fig. 9). There are several candidate E. cautella FAT transcripts (EP_unigene 3, 4, 6, 8, 9, 13, 14 and 15, and EP_contig 45366) that did not form a clade with any other FAT transcripts of the A. ipsilon or S. inferens transcriptome datasets (Fig. 9).

Candidate pheromone degrading enzymes in the E. cautella PG

Pheromone molecules would be potentially harmful to insects if they remained on the ORs after they had stimulated the ORNs. Many studies emphasize that there are mechanisms to protect the ORNs using ODEs [67], including esterases [67–69], aldehyde oxidases [70–72], cytochrome P450 [73–75], carboxyl esterases (cxe) [67] and glutathione S-transferase (GST) [76], which occur in major chemosensory tissues, including the terminal abdominal segment [77]. In general, the esterase gene family consists of three major groups: intracellular (highly expressed in antenna and involved in detoxification), neuro/developmental (neural tissues in antennae) and secreted esterases (expressed in different tissues and associated with specific hormonal and pheromonal functions) [32]. The secreted esterase class contains five major subclasses (glutactin, juvenile hormone (JH) esterases, JHEs-like enzymes, β-esterases and semiochemical esterases), and the ODEs are members of the semiochemical esterases, which are potentially involved in the degradation of pheromone compounds and plant volatiles [32]. The secreted esterase are of three different types, antennal enriched (ODEs in antenna), both antennal and PG-enriched (pheromone degradation) and esterases expressed throughout the body (not pheromone specific) [32].

In the present study, we identified 36 transcripts predicted to encode esterases in the E. cautella PG, and the BLASTx results showed that they shared very high amino acid identities with the esterases of S. exigua, S. littoralis, B. mori and D. plexippus (Additional file 9: Table S8). By comparing E. cautella cxes with S. littoralis [67, 78] and S. inferens [40] we identified the E. cautella cxes that are known to be both antennal and PG-enriched, including cxe 4, cxe5, cxe10, cxe11, cxe13 and cxe16 (Fig. 10). All of the esterase transcripts were present at low abundance levels (less than 80 RPKM) in the PG transcriptome (Additional file 9: Table S8); however, cxe13 (EP_contig_15395_AE, EP_contig_17850_AE and EP_contig_ 28382_AE) had the most highly expressed esterase transcript level (RPKM: 58) (Additional file 9: Table S8). Durand et al. [78] reported the ubiquitous expression of cxe13 in S. littoralis with a specific role in pheromone processing. Homologous cxe13s also reported in Antheraea polyphemus and Popilia japonica were found to degrade the pheromone in vitro [68]. Thus, we assumed that cxe13 has a specific role in E. cautella pheromone processing and degradation.

Fig. 10

Maximum likelihood (ML) tree of insect esterases (cxes)

To assign putative functions and correct identifications, an esterase phylogenetic tree was constructed using 36 E. cautella transcripts and other insect (Drosophila melanogaster, Apis mellifera, A. polyphemus, B. mori, S. lottoralis, A. ipsilon and S. inferens) esterases (Fig. 10). The phylogeny showed that cxe10 (EP_Unigen_10), cxe11 (EP_Unigen_5, EP_Unigen_16 and EP_contig_82246), cxe13 (EP_contig_15395, EP_contig_17850 and EP_contig_ 28382), cxe18 (EP_Unigen_7), cxe14 (EP_contig_67903), cxe26 (EP_contig_33478) and cxe19 (EP_Unigen_14) clustered with the corresponding cxes of S. littoralis [78] and S. inferens [40] (Fig. 10). The phylogeny also revealed that EP_contig_38969, EP_contig_57765 and EP_contig_51966 clustered with the JH esterase of D. melanogaster (ACZ94438), integumental esterase of A. polyphemus (AAM14416) and neuroligins of A. mellifera (NP_001139211), respectively (Fig. 10).

Candidate pheromone carrier proteins in the E. cautella PG

In insects, the odorant binding proteins and the chemosensory proteins are involved in olfaction and contact chemosensation [79, 80]. Specific OBPs that are involved in pheromone binding and transport are called PBPs [26, 79]. In moths, OBPs are divided into three main classes based on sequence alignment and their localization in the insect body. PBPs are preferentially expressed in pheromone-sensitive sensilla trichodea from male antennae, GOBPs are mainly found in the female antennae, in particular in the plant odor sensitive sensilla basiconica [81–83] and ABPXs represent the third class [84]. Based on the sensillar distribution, PBPs may be involved in pheromone binding, GOBPs in binding the plant volatiles and ABPXs in binding general odorants. The functions of OBPs are to solubilize the hydrophobic odorant molecules in the aqueous lymph surrounding the dendrites and to protect them from the ‘degrading esterases’ circulating in the lymph [85]. Additionally, the OBPs deliver the odorant stimuli molecules to specific ORs by releasing the odorant upon contact with membrane structures [86].

Although CSPs are expressed all over the insect’s body, they exist mainly in the legs and in contact chemosensory sensilla. CSPs consist of polypeptide chains of about 110 amino acids with a molecular weight of 12–13 kDa. OBPs have six highly conserved cysteines, whereas CSPs have only four cysteines [79]. Several studies have shown that moth sex pheromones are protected against degradation until they are released from the female PG, and it has been proposed that OBPs and CSPs participate in this process [87, 88].

In the E. cautella PG, we identified transcripts of 7 CSPs and 17 OBPs (Additional file 10: Table S9), all containing the typical insect OBP [88, 89] or CSP sequence motifs [87], respectively. One CSP transcript, EP_unigene_4_CSP, appears to be expressed at an extremely high level (cumulative RPKM: 148,880) in the PG and has a relatively high abundance of transcripts in the PG transcriptome (Additional file 10: Table S9). Phylogenetic analysis shows EP_unigene_4_CSP clustered with B. mori CSP5 (BmorCSP5) [89], H. virescens CSP [38] and A. ipsilon CSP8 (AipsCSP8) [39] (Additional file 12: Figure S11). AipsCSP8 shows a high expression level in the PG and has extremely abundant transcripts in the A. ipsilon PG [39]. Previously, RNAi studies suggested a novel role for a CSP5 in the development of the embryonic integument in A. mellifera and were found to be highly expressed in the ovary [90]. All OBPs have very low expression levels in the PG, but EP_contig_6721_OBP and EP_contig_8460_OBP were comparatively highly expressed OBP transcripts (RPKM: 103) (Additional file 10: Table S9). To assign putative identifications, an OBP phylogenetic tree was constructed with E. cautella OBP transcripts and the B. mori OBPs [91] (Fig. 11). The phylogeny identified E. cautella OBPs homologous to B. mori OBP39 (EP_unigene_11_OBP), OBP37 (EP_contig_8460), OBP44 (EP_contig_2298), OBP43 (EP_contig_6721), OBP31 (EP_unigene_2), OBP20 (EP_unigene_3), OBP18 (EP_unigene_4), OBP17 (EP_unigene_1), OBP15 (EP_unigene_8), OBP14 (EP_unigene_9) and OBP1 (EP_unigene_10) (Fig. 11). A six-cysteine signature is the most typical feature of classical OBPs [88, 89] and E. cautella OBPs carry most of the conserved cysteine residues (data not shown). The B. mori OBPs reported above were found to express in multiple tissues including the terminal abdominal segments (ovary, hind gut, fat body, Malpighian tubule and PG) [91]. Further studies are needed to clarify the roles of the CSPs and OBPs in protecting against the degradation of the E. cautella pheromone prior to its release from the female PG.

Fig. 11

Maximum Likelihood (ML) tree of the Odorant Binding Proteins (OBPs)

We found very low expression levels of three ABP transcripts (less than 60 RPKM) and two PBP transcripts (less than 20 RPKM) in the E. cautella PG (Additional file 10: Table S9). Widmayer et al. [77] detected PBP2 in the H. virescens ovipositor tip, and we found a homologous sequence in the E. cautella PG (EP_contig_73451_PBP; hereinafter EcauPBP2) showing an 89.7 % amino acid identity with A. transitella (ACX47890). In H. virescens, the response to the major pheromone component (Z11-16: Al) is mediated by the PBP2 and the pheromone receptor HR13 [77].

The expression of odorant receptor proteins has been shown to be necessary and sufficient for odor detection in insects [25]. Widmayer et al. [77] detected the pheromone receptors HR2, HR6 and HR13 in the H. virescens ovipositor tip, and HR13 along with PBP2 mediated abdominal responses to the emitted pheromones. In the E. cautella PG, we identified 21 putative ORs, 3 candidate IRs and 2 candidate SNMPs (Additional file 11: Table S10). To assign putative identifications, an OR phylogenetic tree was constructed using E. cautella OR transcripts, H. virescens HR2, HR6 and HR13 [77] and H. armigera ORs [92] (Additional file 13: Figure S12). Based on the phylogenetic analysis, EP_unigene_5_OR (hereinafter EcauOR13) is closely related to the pheromone receptor of H. virescens H13 [77] and OR13 of H. armigera (HarmOR13) [92] (Additional file 13: Figure S12). Recently, Liu et al., [92] reported HarmOR13 responds to the H. armigera pheromone compound (Z11-16: Ald) using calcium imaging studies, and they also found significant gene expression levels in the terminal abdominal segments (TASs) of H. armigera. Hence, we assume that in the E. cautella PG, the response to the pheromone component (Z9,E12-14:OAc) is mediated by the pheromone binding protein, EcauPBP2 and the receptor type EcauOR13, and it may have an important role in the transport and release of the pheromone molecule. Further studies on the functional characterization of EcauPBP2 and the receptor type EcauOR13 will be necessary to prove the hypothesis. Phylogenetic comparisons of E. cautella ORs with those of H. armigera [92] identified the E. cautella receptor proteins expressed in the terminal abdominal segment (Additional file 13: Figure S12). All the receptor protein transcripts were present in very low abundance (less than 50 RPKM) in the PG (Additional file 11: Table S10). It is noteworthy that a SNMP, EP_contig_ 1371_SNMP, which has a greater than 70 % amino acid identity with an O. nubilalis SNMP (GenBank: ADQ73889), was expressed highly in the TAS (RPKM: 102) (Additional file 11: Table S10). Further studies on these PG-expressed ORs, SNMPs and IRs involved in E. cautella pheromone binding and transport need to be performed.

Conclusions

The tropical warehouse moth, E. cautella, listed as a major storage pest, is a serious threat to the date and chocolate factories of the Middle East and Europe. Our study provides comprehensive information on the pheromone molecules, biosynthetic pathways and genes expressed in the PG that are related to pheromone biosynthesis, degradation, transport and release. Our study provides information on the E. cautella sex pheromone and precursors in the PG, and shows two possible pheromone biosynthetic pathways. Both pathways initiate from C16:acyl-CoA, and one involves ∆14 and ∆9 desaturation to generate Z9,E12-14:acyl, while the other pathway involves the chain shortening of C16:acyl-CoA to C14:acyl-CoA, followed by ∆12 and ∆9 desaturation to generate Z9,E12-14:acyl-CoA. Finally, reduction and acetylation generate Z9,E12-14:OAc. Using the Illumina sequencing of the PG transcriptome, we identified candidate genes: PBAN receptor isoforms A and C, FATPs, ACCs, FASs, DESs, several β-oxidation enzymes, FARs and FATs. Two transcripts, EP_unigene_9_1286/EP_contig_343 and EP_contig_349/EP_contig_145 are highly expressed, form a cluster with moth desaturases and might be involved in the ∆14 or ∆12 and ∆9 desaturation processes. The highly expressed β-oxidation enzymes, dehydrogenases, oxidases, hydratases, thiolases, enoyl and dienoyl isomerases, which are involved in the chain shortening, have been identified. Three possible candidate FARs have also been identified, EP_contig_72742, EP_contig_45618 and EP_contig_79681, which form a cluster with moth pgFAR, and thus, might be involved in the reduction step of Z9,E12-14:acid to Z9,E12-14: alcohol. Two possible FATs, EP_unigene_7 and 10, which clustered with yeast alcohol acetyltransferases, are good candidates for gene expression studies. We found many promising candidate PBEs, and further research using heterologous gene expression or RNAi could reveal the significance of these genes in E. cautella pheromone biosynthesis. Several candidate esterases have also been identified, including cxe13 (contig 15395; 17850 and 28382), which may be involved in signal inactivation by removing the pheromone molecules. The CSP (EP_unigene_4) is the most highly abundant transcript of the E. cautella PG, and, together with two OBPs (Contig 6721 and 8460) and one EcauPBP2 (EP_contig_73451), it may have an important functional role in protecting sex pheromones from the activities of esterases, as well as in the transport and release of the pheromone molecules. The ORs, EcauOR13 (EP_unigene_5) and EcauPBP2, meditate abdominal responses to the emitted pheromone in E. cautella, and may have important roles in the transport and release of the pheromone molecules. Our study provides strong background information on the enzymes involved in pheromone biosynthesis that will be useful for the in vitro production of E. cautella sex pheromones. The study also provides information on novel genes involved in the transport, release and degradation of pheromone compounds, increases the understanding of the sex pheromone detection system, and it may provide potential targets for disrupting the pheromone-based communication system in E. cautella for control purposes.

Methods

Chemicals

C14:COOMe and C16:COOMe were purchased from Sigma. (Z,E)-9,12- tetradecadienyl acetate was purchased from Pherobank. E11-13:OH, E14-16:COOMe, Z9-14:COOMe, E12-14:COOMe, Z9-16:COOMe, E9-16:COOMe, Z11-16:COOMe, E11-16:COOMe and Z9E12-14:COOMe were purchased from Pest Control of India Private Limited (Mumbai, India). All standard compounds were of > 98.0 % purity and diluted in n-hexane (250 ng/μl).

Insects

The E. cautella individuals (dried date fruit strain) were originally collected from the Al Hasa date factory (Saudi Arabia) (25°38′ N, 49°60′ E), and were established on an artificial diet comprised of dried broken wheat, peptone and sucrose as the main components. Pupae were sexed, and the female pupae were placed in cages at 24 ± 2.0 °C under a L16:D8 photoperiod. The female pupae were collected separately, and the newly emerged adults were maintained in a vial under the same conditions.

Sex PG extraction and fatty-acyl precursor analysis

The terminal abdominal segments (TASs) (segments 8–10) of individual virgin E. cautella one day before adult eclosion, and of 0-, 1-, 2-, and 3-day-old female moths at mid-scotophase, were excised with micro-scissors. Each gland was extracted for 30 min at room temperature (RT) in a glass insert vial containing 50 μL n-hexane (Sigma) and 250 ng/μL of E11-13:OH as an internal standard (IS). The individual PG extracts were stored at −20 °C until GC-MS analysis.

Total lipid and residue extracted from E. cautella PG were subjected to base methanolysis to convert fatty acyl moieties to the corresponding methyl esters [93]. Two glands were homogenized with a glass rod and transferred into a conical glass vial, and the lipid content was extracted in 500 μL methanol: chloroform (1: 2, v: v), vortexed vigorously for a few minutes and incubated at RT for 30 min. Later, the organic phase was transferred into a new glass tube, and the solvent was evaporated under a gentle stream of nitrogen. Then, 1 mL of 2 % H₂SO₄ (in methanol) was added and incubated at 90 °C for 1 h. Later, 1 mL of milliQ water and 1 mL of n-hexane were added and vortexed vigorously for a few seconds. The upper hexane portions were then transferred into a new glass vial and stored at −20 °C prior to GC-MS analysis.

Both pheromone extracts and the methyl ester samples were subjected to GC-MS analysis on a Agilent 7850A GC coupled to a mass detector (Agilent 5975C) and equipped with a medium-polar INNOWax column (100 % polyethylene glycol, 30 × 0.25 mm I.D., film thickness 0.25 mm, Agilent Technologies, USA). The GC-MS was operated in electron impact mode (70 eV), the injector was configured in split-less mode at 220 °C, and helium was used as carrier gas (velocity: 30 cm/s). The oven temperature was set to 80 °C for 1 min, then increased at a rate of 10 °C/min up to 210 °C, followed by a hold at 210 °C for 15 min, and then increased at a rate of 10 °C/min up to 230 °C, followed by a hold at 230 °C for 20 min.

RNA isolation, cDNA synthesis and library construction

The TASs of ~100 virgin 2- to 3-day-old female E. cautella moths at mid-scotophase were excised. The total RNA of E. cautella PGs was prepared using a NORGEN purification kit (NORGEN Biotek Corp., Canada). Purification is based on spin column chromatography using Norgen’s proprietary resin as the separation matrix, and the procedures were performed according to the manufacturer’s instructions. The quantity and quality of the total RNA was validated using Qubit® 2.0 Fluorometer (Invitrogen, Life Technologies), and the RNA integrity was further confirmed using the 2100 Bioanalyzer (Agilent Technologies) with a minimum RNA integrated number value of 6.8.

The paired-end cDNA libraries were prepared using Illumina protocols and sequenced on the Illumina HiSeq platform. Briefly, the cDNA library was constructed using a TruSeq™ RNA Kit (Illumina Inc.), which consists of mRNA purification and fragmentation from total RNA, synthesizing first and second strands of cDNA, performing cDNA end repair and adenelylating the 3′ ends, followed by adapter ligation and cDNA fragment enrichment. These products were purified and enriched using PCR to create the final cDNA library. Finally, the cDNA library quantity was validated using a Qubit® 2.0 Fluorometer (Invitrogen, Life Technologies), while the quality was validated using an Agilent Technologies 2100 Bioanalyzer prior to the HiSeq Illumina sequencing.

Illumina sequencing

HiSeq Illumina sequencing was performed at the core sequencing facility of the King Abdulla University of Science and Technology (KAUST), Jeddah, Saudi Arabia. The insert size of the library was ~306 bp. Image deconvolution and quality value calculations were performed using the Illumina GAPipeline1.3. All sequencing reads were submitted to the SRA of NCBI under the accession number SRX646348.

Sequence pre-processing, assembly and analysis

A quality control step was first performed on raw sequencing reads using the NGS QC Toolkit [94]. Standard RNA adapter sequences and regions of poor quality were clipped using the CLC Genomic Server and its tool ‘Trim Sequences’. The de novo assembly was performed by the CLC Genomics Server using the scaffolding option and the mapping reads back to transcripts option. This Transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession GBXH00000000. The resulting de novo assembled transcripts were locally searched against the non-redundant (nr) protein database using the BLASTx algorithm (e ≤ 0.001) implemented in the standalone version of the blast + tool [95] and stored in the BLAST archive format (ASN.1). Later, the results were parsed into the required format (XML, tabular, pairwise) using the blast_formatter tool. XML BLASTx results were imported into the BLAST2GO annotation tool. The RPKM values were calculated for assembled transcripts based on their mapping data according to the formula published by Mortazavi et al. [96].

Gene identification and functional annotation

Following the assembly, each transcript was identified by local or web-based searches using the BLASTx and BLASTn programs of NCBI [97]. Blast hits with e-values less than 1.0E-⁵ were considered as significant [98], and the genes were putatively assigned to each contig based on the BLASTx hit with the highest score value. The BLAST XML files were uploaded to BLAST2GO and the mapping, gene annotation, INTERPRO and KEGG analyses were performed as with BLAST2GO [99, 100]. Each gene was checked in terms of molecular function, biological process or cellular component.

Transcripts containing errors leading to misassemblies were edited using Geneious v7.1.5 (www.geneious.com/), de novo assemblies of isotigs were performed and the open reading frame (ORF) of each unigene was determined using the ORF finder tool (NCBI). INTERPRO analysis terms were assigned by BLAST2GO [101] through a search of the nr databases. To annotate the pooled assembled transcriptome, we performed a BLAST search against the nr databases of NCBI, UniProtKB and KEGG using an e-value cut-off of 1.0E5.

Comparative analysis of PG transcriptome

The A. ipsilon PG transcriptome data were downloaded from NCBI (SRX189143) and assembled in the CLC Genomics Server. The H. virescens PG ESTs (14,112 with accession numbers: GR958232-GR972305 and GT067784-GT067747 [38], and the B. mori PG ESTs (10,501 with accession number: BP184340-BP182009, AV404455-AV403746 and DC552314-DC544856) were downloaded from the dbEST database of NCBI (http://www.ncbi.nlm.nih.gov/nucest) and saved as FASTA files. The comparative analyses of E. cautella, A. ipsilon [39], H. virescens [38] and B. mori [39] PG transcripts were performed based on the best bidirectional hit results (first 10 blast hits) (reciprocal BLASTn, e-value less than 1.0E − 6).

Identification of candidate genes involved in E. cautella pheromone biosynthesis

The search for PBEs in our NGS dataset was based on the candidate genes involved in the pheromone biosynthesis in B. mori. We focused on the following target genes: PBANs, ACCs, FASs, desaturases, β-oxidation enzymes, FARs and FATs.

Identification of putative genes involved in fatty acid transport and pheromone degradation

A fatty acid transport protein, BmFATP, was identified from the PG of the silkworm B. mori, which produces a Type-I sex pheromone (bombykol) [47]. BmFATP was shown to facilitate the uptake of extracellular fatty acids into PG cells for the synthesis of bombykol. We performed BLASTx and BLASTn searches to identify E. cautella FATP (EcFATP) genes in the E. cautella PG NGS dataset.

There were earlier reports that esterases may play a major role in pheromone degradation [67–69]. Therefore, we performed BLASTx and BLASTn searches to identify candidate esterase genes in the E. cautella PG assembled NGS dataset.

Identification of putative genes involved in pheromone transport

Genes encoding OBPs and CSPs were identified through BLASTx and BLASTn searches, as well as by the “OBP sequence motif” C1-X15-39-C2-X3-C3-X21-44-C4-X7-12-C5-X8-C6 [79, 87–89, 91] and the “CSP sequence motif” C1-X6-8-C2-X16-21-C3-X2-C4 [87, 89]. Candidate ORs, IRs and SNMP genes were identified by BLASTx and BLASTn searches. Sequence alignments were performed using the ClustalX program [102].

Phylogenetic analyses

E. cautella desaturase and FAR nucleotide sequences were used as query (BLASTx) in the GenBank database, and the desaturase and FAR sequences from different insect species and their amino acids were retrieved for tree construction. The similarity analyses of DNA and protein sequences and a multiple-sequence alignment were performed using the ClustalX program [102], followed by manual inspection. For the phylogenetic analyses, phylogenetic reconstructions were performed using the Geneious tree builder v7.1.5 (www.geneious.com/). The neighbor-joining algorithm analysis was computed using amino acid sequences (Geneious tree builder, Pam250, Jukes-Cantor and Global alignment).

A dataset of esterase sequences was created by retrieving amino acid sequences from NCBI using BLASTx searches of E. cautella PG esterases, and maximum likelihood trees were constructed using MEGA v6.0 [103]. Similarly, acetyltransferase, CSP, OBP and OR sequences were retrieved from the NCBI database, and maximum likelihood trees were constructed using MEGA v6.0 [103]. The NCBI accession number for each gene is provided in the tree.