Skip to main content

Evolution of the codling moth pheromone via an ancient gene duplication

Abstract

Background

Defining the origin of genetic novelty is central to our understanding of the evolution of novel traits. Diversification among fatty acid desaturase (FAD) genes has played a fundamental role in the introduction of structural variation in fatty acyl derivatives. Because of its central role in generating diversity in insect semiochemicals, the FAD gene family has become a model to study how gene family expansions can contribute to the evolution of lineage-specific innovations. Here we used the codling moth (Cydia pomonella) as a study system to decipher the proximate mechanism underlying the production of the ∆8∆10 signature structure of olethreutine moths. Biosynthesis of the codling moth sex pheromone, (E8,E10)-dodecadienol (codlemone), involves two consecutive desaturation steps, the first of which is unusual in that it generates an E9 unsaturation. The second step is also atypical: it generates a conjugated diene system from the E9 monoene C12 intermediate via 1,4-desaturation.

Results

Here we describe the characterization of the FAD gene acting in codlemone biosynthesis. We identify 27 FAD genes corresponding to the various functional classes identified in insects and Lepidoptera. These genes are distributed across the C. pomonella genome in tandem arrays or isolated genes, indicating that the FAD repertoire consists of both ancient and recent duplications and expansions. Using transcriptomics, we show large divergence in expression domains: some genes appear ubiquitously expressed across tissue and developmental stages; others appear more restricted in their expression pattern. Functional assays using heterologous expression systems reveal that one gene, Cpo_CPRQ, which is prominently and exclusively expressed in the female pheromone gland, encodes an FAD that possesses both E9 and ∆8∆10 desaturation activities. Phylogenetically, Cpo_CPRQ clusters within the Lepidoptera-specific ∆10/∆11 clade of FADs, a classic reservoir of unusual desaturase activities in moths.

Conclusions

Our integrative approach shows that the evolution of the signature pheromone structure of olethreutine moths relied on a gene belonging to an ancient gene expansion. Members of other expanded FAD subfamilies do not appear to play a role in chemical communication. This advises for caution when postulating the consequences of lineage-specific expansions based on genomics alone.

Background

Establishing the origin of genetic novelty and innovation is central to our understanding of the evolution of novel traits. While genes can evolve de novo in the genome, the most common mechanism involves duplications of existing genes and the subsequent evolution of novel properties harbored by the encoded gene products. Expansions provide opportunities for specialization and evolution of novel biological functions within a lineage, including the breadth of expression via modifications of promoter architecture [1]. The size of gene families is influenced by both stochastic processes and selection, and particularly large differences in genetic makeup can be indicative of lineage-specific adaptation and potentially associate with traits contributing to phenotypical differentiation between groups [2, 3]. However, the precise functional consequences of such amplifications remain frequently unclear, even if expansions show readily discernible patterns. Therefore, deciphering the underlying genetic and molecular architecture of new phenotypic characters is necessary for a complete understanding of the role played by the accumulation of genetic variation through gene duplication.

For organisms relying on chemical communication, evolution of the ability to produce and detect a new type of molecule could allow for the expansion of the breadth of available communication channels and provide a medium with no or limited interference from other broadcasters. Since the identification of bombykol by Butenandt and co-workers in 1959 [4], a countless number of studies have contributed to revealing the diversity of fatty acid derivatives that play a pivotal role in the chemical communication of insects. As a consequence of homologies with well-studied metabolic pathways, our understanding of the molecular basis of pheromone biosynthesis from fatty acid intermediates has greatly advanced over the past two decades, highlighting the role of several multigene families [5]. Insect semiochemicals are synthesized in specialized cells in which fatty acyl intermediates are converted in a stepwise fashion by a combination of desaturation, chain-shortening, and chain-elongation reactions followed by modifications of the carbonyl group, to cite a few possible steps ([5, 6] and references therein). Desaturation appears particularly important and contributes to producing the great diversity of structures observed in insect pheromones. This derives from the properties of the enzymes that are central to many uncanonical fatty acid synthesis pathways seen in insects. These enzymes can exhibit diverse substrate preference, introduce desaturation in either or both cis (Z) and trans (E) geometry, and give rise to variation in chain-length double-bond position, number, and configuration.

In insects, the group of proteins responsible for catalyzing desaturation reactions are fatty acyl-CoA desaturases (FADs). These membrane-bound acyl-lipid desaturases are biochemically and structurally homologous to the desaturases ubiquitously found in animals, yeast, fungi, and many bacteria where they play important basic biological functions in lipid metabolism and cell signaling and contribute to membrane fluidity in response to temperature fluctuation [7]. In the past decades, the integration of molecular and phylogenetic approaches has greatly advanced our understanding of the function of FAD genes in the biosynthesis of mono- and poly-unsaturated fatty acids in a range of organisms. Moreover, an ever-increasing number of acyl-CoA desaturase genes have been functionally characterized, demonstrating mechanistically their crucial role in the biosynthesis of pheromone and semiochemicals in Drosophila fruitflies, bees, wasps, beetles, and lepidopteran species. Variation in the number and expression of acyl-CoA desaturase genes have been shown to affect the diversity of pheromone signals between closely related species [8,9,10,11,12]. The family is characterized by multiple episodes of expansion and contraction that occurred during the evolution of insects [13, 14]. Consequently, the FAD gene family has become a model to study how structural and regulatory changes act in concert to produce new phenotypes.

Among the taxa available to study the molecular basis of pheromone production in an evolutionary framework, leafroller moths (Lepidoptera: Tortricidae) provide a model system of choice [15]. Leafrollers represent one of the largest families in the Lepidoptera with over 10,000 described species [16]. Their larvae feed as leaf rollers, leaf webbers, leaf miners, or borers in plant stems, roots, fruits, or seeds. Many tortricid species are important pests and, due to their economic impact on human society, became the target of many pheromone identification studies. In 1982, Roelofs and Brown published a comprehensive review on pheromones and evolutionary relationships among Tortricidae [17]. These authors related the patterns in pheromone diversity in a biosynthetic perspective to different postulated phylogenies of the Tortricidae, and specifically the two major groups within Tortricidae, Tortricinae and Olethreutinae. Based on the pheromone identifications available at the time, it was suggested that species in the Tortricinae use mostly 14-carbon pheromone components (acetates, alcohols, and aldehydes) whereas species in the Olethreutinae subfamily use mostly 12-carbon compounds. Building on recent advances towards a robust molecular phylogeny of Tortricidae [18, 19] and incorporating the information for the pheromones identified in 179 species and sex attractants reported for an additional 357 species, we show that the proposed dichotomy is well-supported by the data currently available (Fig. 1). Furthermore, a majority of species in Tortricinae use pheromone components with double bonds in uneven positions, i.e., ∆9 isomers or ∆11 isomers. By contrast, the Olethreutinae pheromones typically contain components with double bonds in even positions (∆8, ∆10), with doubly unsaturated ∆8∆10:12C fatty acyl chains being one of the signature structures of the subfamily. Roelofs and Brown [17] suggested that the use of a Lepidoptera-specific ∆11-desaturase acting on myristic acid (C14) could account biosynthetically for most of the pheromone compounds found in Tortricinae. This hypothesis was later confirmed with the functional characterization of FADs expressed in the female pheromone gland of several Tortricinae representatives. These include a desaturase that makes only the E11-isomer in the light brown apple moth Epiphyas postvittana (E11-14, E11-16, and E9E11-14) [20], as well as desaturases from the redbanded leafroller moth Agryrotaenia velutinana [21] and the obliquebanded leafroller moth Choristoneura rosaceana [22] which both produce a mixture of Z/E11-14:Acids. In the case of Olethreutinae, ∆11 desaturation followed by chain-shortening could account for the ∆9:12C compounds found in several tribes of Olethreutinae. On the other hand, the biochemical pathways leading to the pheromone components with double bonds in even positions, the ∆8 and ∆10 as well as the doubly unsaturated ∆8∆10:12C compounds in the Olethreutinae, were not obvious.

Fig. 1
figure1

Phylogeny of Tortricidae and their associated female sex pheromone components. (left) The maximum likelihood tree was obtained for predicted nucleotide sequences of Tortricidae species (7591 aligned positions). The species represented comprise all tortricids for which pheromones or attractants have been reported plus some outgroups. Typically, one representative species was chosen per genus and contributed molecular evidence for all species in the genus. Outgroup species are represented by red branches whereas tortricid species from the same tribe are represented by branches of the same color. Branch support values were calculated from 1000 replicates using the Shimodaira-Hasegawa-like approximate ratio test (SH_aLRT) and ultrafast bootstrapping (UFboot). Support values for branches are indicated by colored circles, with color assigned based on thresholds of branch selection for SH-aLRT (80%) and UFBoot (95%) supports, respectively. The major subfamilies and represented tribe names are indicated (Phric: Phricanthini; Schoeno: Schoenotenini). (right) Heatmap representing the presence/absence of unsaturated fatty acid structure in bioactive molecules. Attractants correspond to compounds found to be attractive in either field or laboratory experiments; pheromone components correspond to sex attractants produced naturally by the organism and with a demonstrated biological activity on conspecific males. Double-bond positions are annotated in Δ-nomenclature without referring to the geometry. Molecular and trait data retrieved from GenBank and the Pherobase, respectively

The codling moth Cydia pomonella (Linnaeus) (Tortricidae: Olethreutinae: Grapholitini) is one of the most devastating pests in apple and pear orchards worldwide [23]. With the goal of disrupting its reproduction, it has received substantial attention and been at the center of numerous studies focused on characterizing its communication system via sex pheromones. Its pheromone, (E8,E10)-dodecadien-1-ol, also known under the common name codlemone, was first identified using gas chromatography in combination with electroantennogram (EAG) recordings [24]. This identification was later confirmed by fine chemical analysis [25, 26]. The codling moth thus provides a relevant system to unravel the molecular pathway associated with the production of the ∆8∆10:12C typical of olethreutines. Previous studies support the hypothesis that the biosynthesis involves the desaturation of a ∆9 monoene intermediate. First, Arn et al. [27] reported the presence of the unusual (E)-9-dodecenol (E9-12:OH) at about 10% of the doubly unsaturated alcohol in pheromone gland extracts and effluvia from C. pomonella females. Although this monoene does not carry any behavioral activity, its occurrence suggested that E9-12:Acyl could be an intermediate of codlemone biosynthesis in a process analogous to the biosynthesis of 10,12-dienic systems via ∆11 monoene intermediates observed in Bombyx mori, Manduca sexta, and Spodoptera littoralis [28,29,30]. The selective incorporation of deuterium-labeled E9-12 fatty acid precursors into codlemone and its direct precursor, E8,E10-dodecadienoate (E8E10-12:Acyl), supported the presence of an unusual ∆9 desaturase in C. pomonella and the biosynthesis of a conjugated diene system via 1,4-desaturation and the characteristic elimination of two hydrogen atoms at the allylic position of the double bond in the monoene intermediate precursor [31] (Fig. 2). Similar conclusions were reached from a replicate study using C. splendana and C. nigricana females [32]. To date, the FAD(s) central to the biosynthesis of codlemone and related fatty acyl derivates with a ∆8∆10 system has not been characterized.

Fig. 2
figure2

Pathway of codlemone biosynthesis. Palmitic acyl (C16) is first formed by the fatty acid synthesis pathway. Lauric acyl (C12) is then produced following two recurring reactions of beta-oxidation. A first desaturation occurs, introducing an E9 double bond between carbons 9 and 10 from the carboxyl terminus. The resulting monoenoic fatty acyl is then the substrate of a second desaturation via the removal of hydrogen atoms at the carbons 8 and 11, causing the formation of a E8E10 conjugated diene system. Finally, an unknown fatty acyl-CoA reductase converts the fatty acyl precursor into the fatty alcohol (E8,E10)-8,10-dodecadienol commonly known as codlemone. The steps characterized in this study are highlighted in red

The availability of a high-quality draft genome for C. pomonella provides an opportunity to comprehensively annotate and analyze relevant genes [33]. Here we annotated a total of 27 FAD genes and performed a phylogenetic analysis, revealing expansions of different ages in the ∆9(KPSE) clade and in the ∆10/∆11(XXXQ/E) clade. Using a transcriptomic approach, we determine the breadth of expression of all FAD genes and identified genes with upregulated expression in the pheromone gland of female C. pomonella, where the biosynthesis of codlemone takes place. We tested the function of these desaturases in heterologous expression systems and identified one FAD gene conferring on the cells the dual desaturase functions playing a key role in the biosynthesis of (E8,E10)-dodecadienol in C. pomonella and a signature structure of many Olethreutinae moth pheromones.

Results

Expansion of FADs in the genome of C. pomonella

We identified candidates potentially involved in fatty acid synthesis by searching in the genome of C. pomonella for genes encoding fatty acid desaturases, which are characterized by a fatty acid desaturase type 1 domain (PFAM domain PF00487). First, we improved the genome annotation by incorporating data from the pheromone gland transcriptome, a tissue which was not part of the panel of tissues used to generate the available annotation of C. pomonella. Searching our improved annotation, we identified 27 genes harboring the signature domain of FADs. While the vast majority of genes identified in our exhaustive search contain open-reading frames of 300 amino acids or longer, a small number of genes (i.e., 2–3) may represent pseudogenes or assembly errors. We adopted the nomenclature proposed by Knipple et al. [34] in which genes are named based on the composition of 4 amino acid residues at a signature motif. Next, we looked at their genomic organization. We found that FAD genes are distributed across 12 of the 27 autosomes, with no FAD genes on either Z or W sex chromosome (chr1 and chr29, respectively) (Fig. 3). Two genes, Cpo_QPVE and Cpo_MATD(2), were found on unplaced scaffolds. As is typical for members of a multigene family, we identify several clusters corresponding to tandemly duplicated genes. Moreover, several genes are present as the single member of the family on a given autosome.

Fig. 3
figure3

Genomic organization of fatty acyl desaturase genes in Cydia pomonella. The positions of the 27 predicted FAD genes were mapped to the genome. Twenty-five genes could be placed on 12 autosomes; 2 genes, Cpo_MATD(2) and Cpo_QPVE, are located on unplaced scaffolds (not drawn). The distribution analysis showed that there are 5 gene clusters which contain 2 or more FAD genes. Names refer to chromosome names in the genome assembly. Chromosomes are drawn to scale and the bands represent gene density, with darker bands depicting gene-poor regions

Phylogenetic analyses place expansions in the ∆9(KPSE) and ∆10/∆11(XXXQ/E) clades

In order to determine which of the genes identified above could be important for codlemone biosynthesis, we assessed the functional classes represented by these genes using phylogenetic analyses. To that end, we aligned the predicted protein sequences of the C. pomonella genes with a panel of FAD sequences from Lepidoptera species including genes for which function has been previously characterized using in vitro heterologous expression. We found representatives of the eight insect acyl-CoA desaturase subfamilies sensu Helmkampf et al. [14] (Fig. 4; Additional file 1: Figure S1).

Fig. 4
figure4

Phylogeny of Lepidoptera FAD genes. The maximum likelihood tree was obtained for the predicted amino acid sequence of 114 FAD genes (805 aligned positions) of 28 species, with branch support values calculated from 1000 replicates using the Shimodaira-Hasegawa-like approximate ratio test (SH_aLRT) and ultrafast bootstrapping (UFboot). Support values for branches are indicated by colored circles, with color assigned based on SH-aLRT and UFBoot supports using 80% and 95% as thresholds of branch selection for SH-aLRT and UFBoot supports, respectively. The major constituent six subfamilies of First Desaturase (A1 to E) and two subfamilies of Front-End (Cyt-b5-r) and Sphingolipid Desaturases (Ifc), respectively, are indicated following the nomenclature proposed by Helmkampf et al. [14]. For First Desaturases, the different shades correspond to the indicated putative biochemical activities and consensus signature motif (if any). Triangles indicate sequences from C. pomonella (see also Additional file 1: Figure S1 for the extended version of this tree). The scale bar represents 0.5 substitutions per amino acid position

All subfamilies form highly supported groups, with the exception of the relationship between Desat A1 and A2, which appear weakly or strongly supported, depending on the set of genes used in our analyses. With a characteristic domain architecture (PF08557 in front of PF00487), Cpo_TYSY encodes a putative Sphingolipid Delta-4 desaturase and is homologous to interfertile crescent in D. melanogaster (Ifc). Two tandemly duplicated genes, Cpo_RIDY and Cpo_DRKD, contain a Cytochrome b5-like heme binding domain (PF00173) in front of the fatty acid desaturase domain and group with putative Front-End desaturases homologous to Cytochrome b5-related in D. melanogaster (Cyt-b5-r). These three genes bear little similarity with the other 24 FAD genes which encode First Desaturases of the subfamilies Desat A1 through E.

Desat C, D, and E are present as single-copy genes in C. pomonella (Cpo_QPVE, Cpo_KSTE, and Cpo_KYPH, respectively). Cpo_QPVE groups with the ∆14 desaturase identified in male and female corn borer pheromone biosynthesis (Lepidoptera: Crambidae) [9, 35].

Desat B, which is particularly expanded in Hymenopterans and in Bombyx mori [14], has only two representatives in C. pomonella. These genes, Cpo_VKVI and Cpo_TPVE, are homologous to the E6 desaturase identified in Antheraea pernyi (Lepidoptera: Saturnidae) [36] and the Z5 desaturase identified in Ctenopseustis obliquana (Lepidoptera: Tortricidae) [37], respectively.

Desat A2 subfamily is characterized by a single-copy gene, Cpo_GATD, which seems to be typical of Lepidoptera. The homolog identified in Choristoneura parallela (Lepidoptera: Tortricidae) is associated with the formation of ∆9 desaturation in saturated acyl moieties of a range of length (C14-C26) [38].

Finally, Desat A1 forms the largest group and experienced a particularly dynamic evolutionary history. With 17 genes, C. pomonella harbors twice as many genes as Bombyx mori (7 in the most recent version of the genome SilkDB 3.0 [39]). While ∆9 C16<C18 (NPVE) is represented by a single-copy gene, we found two significant expansions in the group corresponding to ∆9 C16>C18 (KPSE) and ∆11, ∆10, and bifunctional (XXXQ/E) FADs. The tandemly arrayed cluster found in chr4 corresponds to an extension within the A1 subfamily encoding ∆9 C16>C18 FADs. Nine genes are found in the Lepidoptera-specific A1 subfamily encoding ∆11, ∆10, and bifunctional FADs. These genes are scattered across 6 private chromosomes, none of which containing members of the other functional classes.

Based on their placement in the FAD gene tree, these genes can be sub-divided into four groups (Fig. 4; Additional file 1: Figure S1). Group1 contains Cpo_SPTQ(1), Cpo_SPTQ(2), Cpo_SATK, and Cpo_TSTQ, which are homologous to the ∆10 desaturase identified in the New Zealand tortricid Planotortrix octo (Lepidoptera: Tortricidae) [40]. Group2 is formed with Cpo_RATE, Cpo_MATE, and Cpo_TPSQ, which are homologous to the FAD with ∆6 activity isolated from Ctenopseustis herana (Lepidoptera: Tortricidae) [12]. Group3 is represented by Cpo_RPVE, which groups with desat2 from Ctenopseustis obliquana and allied species (Lepidoptera: Tortricidae) [12] and Lca_KPVQ from Lampronia capitella (Lepidoptera: Prodoxidae) [41]. Finally, Cpo_CPRQ falls near Dpu_LPAE from Dendrolimus punctatus (Lepidoptera: Lasiocampidae) [42]. Based on our phylogenetic reconstructions, several genes are plausible candidates in the biosynthesis of codlemone. First, genes from the ∆9 C16>C18 (KPSE) subfamily have been implicated in pheromone biosynthesis and carry activities similar to those expected in C. pomonella. Specifically, the Dpu_KPSE from Dendrolimus punctatus (Lepidoptera: Lasiocampidae) produces a range of ∆9-monounsaturated products including Z9- and E9-12:Acyl [42]. This exemplifies that enzymes in this subfamily can introduce cis and trans double bonds in lauric acid (C12). Furthermore, when supplemented with Z7- and E7-14:Acyl, Dpu_KPSE can introduce a second desaturation to produce ∆7∆9-14:Acyl, illustrating the formation of conjugated double bonds by this type of FAD. The expansion of the KPSE clade we report in C. pomonella is an interesting coincidence. Second, the abundance of representatives in the ∆11, ∆10, and bifunctional (XXXQ/E) subfamily supports the possibility that one or several could be involved in codlemone biosynthesis. This group contains many genes involved in the biosynthesis of conjugated double bonds, e.g., Bmo_KATQ in Bombyx mori (Lepidoptera: Bombycidae) [43], Mse_APTQ in Manduca sexta (Lepidoptera: Sphingidae) [44], Lca_KPVQ in Lampronia capitella (Lepidoptera: Prodoxidae) [41], and Dpu_LPAE in Dendrolimus punctatus [42]. Cpo_RPVE groups with Lca_KPVQ, which encodes an enzyme involved in the desaturation of 16:Acyl and Z9-14:Acyl to produce the conjugated Z9Z11-14:Acyl pheromone precursor of L. capitella [41]. Another interesting candidate is Cpo_CPRQ, which clusters near Dpu_LPAE, an enzyme capable of catalyzing the production of E9Z11-16:Acyl and E9E11-16:Acyl [42].

Patterns of gene expression identify Cpo_CPRQ and Cpo_SPTQ(1) as candidates

To identify candidates for codlemone production, we compared the expression levels of the twenty-seven genes found in the genome using published RNA sequencing (RNA-Seq) data from various tissues and life stages augmented with two new pheromone gland data sets (this study). Since codlemone is found in the pheromone gland of females, which is located at the tip of the abdomen, we hypothesized that its desaturase(s) would be enriched or even exclusively expressed in this sex and in the population of cells forming this tissue.

We found strong qualitative differences in the transcription profile of the 27 examined FAD genes (Fig. 5). The genes divide into three clusters. The first two correspond to genes that appear ubiquitously expressed. These include the two putative Front-End desaturases Cpo_DRKD and Cpo_RIDY (Cyt-b5-r), the Sphingolipid ∆4 desaturase Cpo_TYSY (Ifc), and five First desaturases, namely Cpo_NPVE (A1; ∆9 C16<C18), Cpo_KPSE (A1; ∆9 C16>C18), Cpo_RPVE (A1; ∆11/∆10), Cpo_KSTE (D), and Cpo_KYPH (E). With the exception of Cpo_RPVE, these genes exhibit a high evolutionary stability characterized by a single-copy gene with rare cases of gains and losses across insects. Their broad expression profiles suggest that these genes fulfill a fundamental metabolic function, making them unlikely candidates for pheromone biosynthesis. By contrast, the other First Desaturase genes found in the third cluster are characterized by a sparse expression profile and greater tissue and/or life-stage specificity: head: Cpo_TPSQ (A1; ∆11/∆10); ovary/female abdomen: Cpo_TPVE (B; ∆5); pheromone gland: Cpo_CPRQ (A1; ∆11/∆10), Cpo_SPTQ(1) (A1; ∆11/∆10); and eggs: Cpo_NPRE (A1), Cpo_VKVI (B; ∆6). These expression patterns are suggestive of more specialized functions. Specifically, their expression profiles and levels indicated that Cpo_CPRQ and Cpo_SPTQ(1) (average FPKM in pheromone gland, 6082.0 and 1451.5, respectively) represent primary candidates for functional testing.

Fig. 5
figure5

Expression profiling of FAD genes in different tissues of Cydia pomonella. The heatmap represents the absolute expression value (log FPKM) of all FAD genes in the corresponding tissues. Genes are identified by their signature motif. Automatic hierarchical clustering of FAD genes distinguishes three clusters: clusters I and II contain genes with ubiquitous expression across tissues and developmental stages; cluster III contains genes with more divergent expression and higher tissue specificity. Automatic hierarchical clustering of tissues indicates that pheromone gland samples have a distinguishable profile, which is characterized by the overexpression of Cpo_CPRQ and Cpo_SPTQ(1). Estimates of abundance values were obtained by mapping reads against the genome. The NCBI SRA accession numbers of all RNA-Seq data sets used are given in Additional file 6: Table S4

Functional characterization demonstrates the E9 FAD activity of Cpo_CPRQ and Cpo_SPTQ(1)

Previous studies have shown that the baker’s yeast Saccharomyces cerevisiae is a convenient system for studying in vivo the function of FAD genes and it has become the eukaryotic host of choice for the heterologous expression of desaturases from diverse sources including moths. Conveniently, yeast cells readily incorporate exogenous fatty acid methyl esters and convert them to the appropriate coenzyme A thioester substrates required by FADs.

First, we cloned Cpo_CPRQ and Cpo_SPTQ(1) into a copper-inducible yeast expression vector to generate heterologous expression in S. cerevisiae. We carried out assays with precursor supplementation to ensure availably of the medium-chain saturated fatty acids, i.e., lauric (C12) and myristic acids (C14), which are typically less abundant. Cpo_CPRQ conferred on the yeast the ability to produce a small amount of E9-12:Acyl and had no detectable activity outside of the C12 substrate (Fig. 6). Cpo_SPTQ(1) showed ∆9 desaturase activity, producing E9- and Z9-12:Acyl, Z9-14:Acyl, and Z9-16:Acyl. Double-bond positions in these products were confirmed by dimethyl disulfide (DMDS) derivatization (Fig. 6). In addition to retention times matching those of synthetic standards, spectra of the DMDS derivatives of the methyl esters displayed the following diagnostic ions: Z9-16:Me (M+, m/z 362; A+, m/z 145; B+, m/z 217), Z9-14:Me (M+, m/z 334; A+, m/z 117; B+, m/z 217), and E9- and Z9-12:Me (M+, m/z 306; A+, m/z 89; B+, m/z 217). No doubly unsaturated products were detected in these experiments.

Fig. 6
figure6

Functional characterization of desaturase activity of candidate genes in yeast. Total ion chromatograms of fatty acid methyl ester (FAME) products of Cu2+-induced ole1 elo1 S. cerevisiae yeast supplemented with saturated acyl precursors and transformed with a empty expression vector (control), b pYEX-CHT-Cpo_CPRQ, and c pYEX-CHT-Cpo_SPTQ(1). Confirmation of the identity of enzyme products was obtained by comparison of retention time and mass spectra of synthetic standards. Confirmation of double-bond positions by mass spectra of DMDS adducts (d, e)

In a second round of experiments, we assayed 7 other First desaturases that could be amplified from pheromone gland cDNA. As predicted from their respective placement in the FAD phylogeny, Cpo_NPVE, Cpo_KPSE, and Cpo_GATD encode ∆9 desaturases (Additional file 2: Figure S2). Cpo_NPVE is a ∆9 desaturase with preference for myristic (C14) and palmitic (C16) acid, showing only weak activity on lauric acid (C12). Cpo_GATD and Cpo_KPSE showed ∆9 desaturase activities mainly on C16 but also some activity on C14. The other 4 desaturases had no detectable activity in the yeast system.

Finally, we supplemented the growth medium with the monounsaturated methyl esters E9-12:Me and Z9-12:Me. The doubly unsaturated codlemone precursor, E8E10-12:Acyl, was inconsistently detected in some replicates of the assays with Cpo_CPRQ. No doubly unsaturated products were detected with the other constructs.

Functional analysis identifies Cpo_CPRQ as the FAD catalyzing ∆8∆10 desaturation

In order to evaluate the presumptive ∆8∆10 desaturase activity of Cpo_CPRQ, we further investigated the activity of the enzyme in an insect cell line. We reasoned that the inconsistent activity of the enzyme in yeast could be due to the cellular environment of the expression system. When we expressed Cpo_CPRQ using the Sf9 cell system and supplementing the culture medium with lauric acid, we could observe the consistent production of E9-12 as well as E8E10-12, demonstrating that Cpo_CPRQ catalyzes the biosynthesis of E8E10-12:Acyl and its monounsaturated intermediate E9-12:Acyl (Fig. 7). As expected from the previous experiments in yeast, we did not detect activity on longer chain substrates. The retention time and mass spectrum of the E8E10-12:Me peak were identical to those of the synthetic standard. Double-bond positions of the conjugated system were confirmed by derivatization with 4-methyl-1,2,4-triazoline-3,5-dione (MTAD) (Fig. 7). Finally, we carried out assays to further evaluate the stereospecificity of the enzyme. When provided with E9-12:Me, Cpo_CPRQ produced the doubly unsaturated product (Fig. 7). In contrast, we found no evidence for activity on the corresponding diastereomer, Z9-12:Me (Fig. 7).

Fig. 7
figure7

Functional characterization of desaturase activity of Cpo_CPRQ in insect cells. Total ion chromatograms of methyl esters (FAME) samples from Sf9 cells supplemented with lauric methyl ester (C12) infected with a empty virus (control) or b recombinant baculovirus expressing Cpo_CPRQ. Cpo_CPRQ produces large amount of E9-12:Me and E8E10-12:Me. Sf9 insect cells infected with bacmid expressing Cpo_CPRQ in the medium in the presence of the monoenic intermediate c (E)-9-dodecenoic methyl ester (E9-12:Me) or d (Z)-9-dodecenoic methyl ester (Z9-12:Me). The retention time and mass spectrum of the E8E10-12:Me peaks observed after addition of 12:Me and E9-12:Me were identical with those of the synthetic standard. The relatively long retention time (compound eluting later than 14:Me) is in agreement with what is expected from a diene with conjugated double bonds. e Spectrum of the E8E10-12:Me peak in b. The relatively abundant molecular ion m/z 210 is in agreement with the expectation for a diene with conjugated double bonds. f Analyses of MTAD-derivatized samples displayed diagnostic ions at m/z 323 (M+), m/z 308, and m/z 180 (base peak), confirming the identification of the conjugated double-bond system

Altogether, these results show that Cpo_CPRQ exhibit both the ability to catalyze the trans desaturation of lauric acid to form E9-12:Acyl and the 1,4-desaturation of the latter to form the E8E10-12 conjugated pheromone precursor of codlemone.

Discussion

In this study, we show the key role of a desaturase with a dual function in the biosynthesis of codlemone. We confirm the results suggested by the labeling experiments reported by Löfstedt and Bengtsson [31] and demonstrate that Cpo_CPRQ can conduct both the initial desaturation of 12:Acyl into E9-12:Acyl, and then transform E9-12:Acyl into E8E10-12:Acyl, the immediate fatty acyl precursor of E8E10-12:OH, the sex pheromone of C. pomonella (Fig. 2). Phylogenetically, Cpo_CPRQ clusters in the Lepidoptera-specific XXXQ/E clade (Fig. 4), a lineage which encodes pheromone biosynthetic enzymes with diverse properties [34]. Several FADs found in that clade have been demonstrated to be involved in the pheromone biosynthesis of many moth species, including the species of Tortricinae using ∆11 (Epiphyas postvittana, Choristoneura rosaceana, Agryrotaenia velutinana) and ∆10 desaturations (Planotortrix octo) [20,21,22, 40]. Our phylogenetic reconstruction indicates that Cpo_CPRQ does not appear orthologous to any of those genes (Additional file 1: Figure S1), suggesting the recruitment of a different ancestral duplicate early in the evolution of Olethreutinae. In addition to Cpo_CPRQ, Cpo_SPTQ(1) clustered in the same clade and had a similar activity profile on lauric acid, although it produces more Z9- than E9-12:Acyl. With an expression at about 25% the level of Cpo_CPRQ, Cpo_SPTQ(1) is the second most highly expressed FAD gene in the pheromone gland, a tissue in which this FAD appears exclusively expressed. Although Cpo_CPRQ appears sufficient to produce E8E10-12:Acyl, we cannot exclude the possibility that Cpo_SPTQ(1) contributes to the production of the monoene intermediate and is involved in the pheromone biosynthesis of codlemone or E9-12:OH, a pheromone component eliciting antennal response but with no apparent behavioral effect [45]. Interestingly, even though they are phylogenetically clustered within the ∆11/∆10-desaturase subfamily, Cpo_CPRQ and Cpo_SPTQ(1) exhibit primarily ∆9-desaturase activity, which likely represent the ancestral state. Whether this represents a conservation or a reversal to the ancestral function awaits further investigations. This finding provides exciting opportunities to study the structure-function relationships in FADs, in particular the structural determinant of regioselectivity.

Expansion of gene families is usually regarded as playing a key role in contributing to phenotypic diversity. The acyl-CoA desaturase gene family is characterized by a highly dynamic evolutionary history in insects [13, 14, 34]. Bursts of gene duplication have provided many opportunities for evolution to explore protein space and generate proteins with unusual and novel functions. In moths, many examples highlight the role played by the diversification in FAD functions in allowing the expansion of the multi-dimensional chemical space available for communication via sex pheromones. Careful re-annotation of the C. pomonella genome and phylogenetic analyses of identified desaturase genes allowed us to identify 27 genes representative of all desaturase subfamilies previously identified in insects. Interestingly, all functional classes previously identified in moths (i.e., ∆5, ∆6, ∆9, ∆10, ∆11, ∆14) have homologs in C. pomonella. This indicates that FAD genes tend to be preserved in the genomes for long periods of time. The availability of several chromosome-level assemblies in Lepidoptera has revealed that the genome architecture is generally conserved, even among distant species, with few to no structural rearrangements, i.e., inversions and translocations, being observed (but see [46]). This is attested by the high level of synteny existing between C. pomonella and Spodoptera litura (Lepidoptera: Noctuidae), two representatives of lineages which shared their last common ancestor ca. 120 MYA [33, 47]. This stability over long evolutionary time means that we can draw tentative conclusions from the genomic organization of gene families. Specifically, the genomic location of FAD genes in C. pomonella suggests that the diversification of this gene family is ancient. First Desaturase genes from all groups sensu Helmkampf et al. [14] have representatives in species that have shared their last common ancestor over 300 MYA. Moreover, we found several expansions with differing patterns within ∆9 C16>C18 (KPSE) and the Lepidoptera-specific gene subfamily ∆11, ∆10, and bifunctional (XXXQ/E). Genes clustered in tandem arrays represent evolutionary recent duplications, as exemplified by the expansion in ∆9 C16>C18 (KPSE). These duplicates are likely the product of imperfect homologous recombination. By contrast, with 9 representatives spread across 6 chromosomes, the ∆11/∆10 subfamily comprises both recent duplicates (e.g., Cpo_SPTQ(1) and Cpo_SPTQ(2)) as well as unlinked FAD genes (such as Cpo_CPRQ). The latter genes are likely to be the products of more ancient events of duplications and chromosomal rearrangements. The location of the majority of genes within this clade suggests that the diversification within this subfamily essential for the pheromone biosynthesis occurred early in the evolution of the Lepidoptera. As genome data from more species become available, it will be possible to follow the evolutionary trajectories of this gene family which has been a key contributor to pheromone evolution.

Expansion of gene families is often followed by a relaxation of selection, which translates into higher rates of fixation of mutations, including non-synonymous nucleotide substitutions [48, 49]. Looking at the impact of selection on desaturase genes in insects, previous studies found little evidence for positive selection but rather purifying selection acting on all desaturase genes in insects, even in expanded subfamilies [14, 34]. This could be indicative that large constraints along the coding sequences to not change sites affecting the proper enzymatic function outweigh the potential for diversifying selection. Desaturases encoded by genes in the ∆11/∆10 (XXXQ/E) lineages display a great variety of enzymatic functions, which are paralleled by very high levels of sequence divergence. The lack of evidence for positive selection could result from a lack of power in the statistical analysis currently available under when analyzing genes selected for very diverse properties [34]. Modification of the regulation of duplicated genes is also expected, with some paralogs specializing to a limited set of tissues or developmental stages [49]. Data from a panel representing a relatively large variety of tissue types and developmental stages reveal a large expression divergence between FAD genes. Interestingly, we see that genes present in the recent expansion (i.e., ∆9 C16>C18 (KPSE)), of which members form tandem arrays in the genome, have low to no expression in the large panel of samples we analyzed. The evidence at hand suggests that these may have limited physiological relevance. This advocates for caution when interpreting the importance and functional consequences of expansions. Some expansions can contribute to lineage-specific adaptation and an increased demand for variability in associated traits. However, it is almost certain that many expansions will not reflect a response to changes and be responsible for the evolution of novel traits but rather have their evolutionary origin in the stochastic nature of the gene duplication process. Functional testing is paramount to advance our mechanistic understanding of the proximate consequences of variation in the size of multigene families.

Previously characterized FADs that catalyze the formation of conjugated fatty acids among moth sex pheromones precursors fall into two groups: (1) double bonds can be inserted via sequential 1,2-desaturation to produce a diene, either by the action of two different desaturases operating consecutively as postulated in Dendrolimus punctatus (Lasiocampidae) [42], or by the same desaturase operating twice with an intermediate chain-shortening step as suggested for Epiphyas postvittana (Tortricidae) [20], L. capitella [41], and Spodoptera litura (Noctuidae) [50]; (2) bifunctional desaturases can introduce the first double bond and then rearrange the monoene to produce a diene with conjugated double bonds via 1,4-desaturation, as seen in Bombyx mori (Bombycidae) [43], S. littoralis (Noctuidae) [51], and Manduca sexta (Sphingidae) [44]. The C. pomonella desaturase falls into the latter category: Cpo_CPRQ introduces first an E9 double bond in the 12:Acyl saturated substrate and then transforms the E9 unsaturation into the E8E10 conjugated double bonds. The independent evolution of similar features in enzymes used by distantly related species offers the opportunity to further study mechanistically the determinants of the switch between 1,2 and 1,4-dehydrogenation.

Recent molecular phylogenetic studies [18, 19] provide strong support for the monophyly of Tortricinae and Olethreutinae, the two most speciose subfamilies within Tortricidae, and the paraphyly of Chlidanotinae. Our tribal-level tree for Tortricidae combined with data on moth pheromones and attractants corroborate the pattern of pheromone components in Tortricinae and Olethreutinae reviewed and discussed by Roelofs and Brown [17]. Tortricinae species use mostly 14-carbon pheromone components with unsaturation introduced via ∆11 desaturation (Fig. 8). Molecular and biochemical studies have confirmed the paramount role of ∆11 FADs in the pheromone biosynthesis of not only tortricines, but the vast majority of ditrysian Lepidoptera (see for instance [13, 41, 52,53,54]). No pheromones are reported for the Chlidanotinae. However, the sex attractants reported for three species include Z9-14:OAc and Z11-16:OAc, which supports the idea that ∆11 desaturation in combination with chain-shortening may be typical of this subfamily. By contrast, the Olethreutinae subfamily uses mainly 12-carbon compounds. The ∆9-12 chain structures used by many species could be accounted for the action of a ∆11 enzyme producing ∆11-14:Acyl as an intermediate following by chain-shortening to ∆9-12:Acyl (Fig. 8). For example, the host races of the larch budworm Zeiraphera diniana (Tortricidae: Olethreutinae: Eucosmini) developing on larch (Larix decidua) or Cembra pine (Pinus cembra) used respectively E11-14:OAc and E9-12:OAc as pheromone (Guerin et al. 1984), and it is parsimonious to hypothesize that polymorphism in the activity of the chain-shortening enzyme is responsible for this difference between host races that show little to moderate genomic differentiation [55]. Our data suggest however that an alternative pathway may exist. Indeed, two of the enzymes characterized in the present study (i.e., Cpo_SPTQ(1) and Cpo_CPRQ) can catalyze the production of ∆9-12:Acyl directly from lauric acid (C12), bypassing the need for prior ∆11 desaturation and chain-shortening (Fig. 8).

Fig. 8
figure8

Postulated biosynthetic steps to account for the typical pheromone components found in the Tortricidae. Fatty acid biosynthetic routes for sex pheromone components of Tortricinae (a) and Olethreutinae (b) are illustrated. Starting from palmitic acyl (C16), the pheromone precursors are produced via a combination of chain-shortening (by limited beta-oxidation) and desaturation steps. Genes encoding the fatty acyl desaturases with the postulated activities have been characterized in this study or previously, with the exception of the enzyme(s) involved in the synthesis of ∆8 monoenes in Olethreutinae

In addition to ∆9-12 chain structures, 12-carbon components with ∆8 and ∆8∆10 systems are frequently used by olethreutines. Specifically, these occur in the three major tribes: Grapholitini, Eucosmini, and Olethreutini (Fig. 1). The identification of the bifunctional desaturase catalyzing the biosynthesis of ∆8∆10-12 fatty acid derivatives sheds new light on the pathway associated with the production of pheromones of a large number of tortricids. Monounsaturated compounds with double bonds in ∆8 position frequently co-occur with conjugated ∆8∆10, including in trace amount in C. pomonella [45]. However, since we did not detect any ∆8 activity towards the lauric acid from Cpo_CPRQ or Cpo_SPTQ(1), we anticipate that a different FAD plays a role in the pheromone biosynthesis of ∆8 monoenes. This is supported by the pheromone bouquets identified for several species of olethreutines. In the podborer Matsumuraeses falcana (Tortricidae: Olethreutinae: Grapholitini), the pheromone consists of a mixture of E8-12:OAc and E8E10-12:OAc, plus minor amounts of (E7,Z9)- and (E7,E9)-dodecadienyl isomers, and the monoenes E10-12:OAc and E10-14:OAc [56]. The major pheromone component in Hedya nubiferana (Tortricidae: Olethreutinae: Olethreutini) is E8E10-12:OAc, but in addition, the pheromone contains relatively large amounts of Z8- and E8-12:OAc [57]. The same holds for Epiblema foenella (Tortricidae: Olethreutinae: Eucosmini) in which the major pheromone component was identified as Z8Z10-12:OAc: in addition to the other three geometric isomers of the diene, the females produce also E8- and Z8-12:OAc but no E9-12:OAc [32]. Biosynthesis of these bouquets suggests the involvement of a FAD with ∆8 or ∆10 activity (Fig. 8). A ∆10 FAD is the key enzyme involved in the production of monounsaturated compounds with double bonds in ∆8 position. Specifically, in Planotortrix sp. and Ctenopseustis sp. (Tortricidae: Tortricinae: Archipini), Z8-14:OAc is biosynthesized via ∆10 desaturation of palmitic acid (C16) followed by chain-shortening, whereas the production of Z10-14:OAc involves the same FAD acting on myristic acid (C14) [12, 40, 58] (Fig. 8a). Our phylogenetic analysis of FAD genes shows that the C. pomonella genome contains several genes clustering with the tortricine ∆10 FADs, indicating that this biochemical activity may in theory still be carried out in olethreutines. Additional work is necessary to elucidate the routes towards the biosynthesis of monounsaturated pheromone components with double bond at even positions.

On a practical point of view, pheromones have gradually been added to the pest control toolkit. In the specific case of Cydia pomonella, in 2010, the annual production of synthetic codlemone—the codling moth pheromone—was estimated at 25 metric tons, which were used for the treatment by mating disruption of 200,000 ha of orchards worldwide [59]. The biological production of moth pheromones in plants or microorganisms such as yeast has been suggested as an environmentally friendly alternative to conventional chemical synthesis [60]. Beyond the stage of conceptualization, the approach is proven, and large-scale production is on the horizon [61,62,63,64]. Identification and characterization of the key enzymes underlying production of pheromone is a crucial step to assemble the building blocks central to the engineering of plant and cell factories. In this context, the identification and characterization of a desaturase involved in pheromone biosynthesis in C. pomonella is of particular interest, as it represents an important step towards the biological production of pheromone for the integrated pest management of several economically important olethreutines.

Conclusions

The fatty acid desaturase gene family is central to the pheromone-based communication of many insects. Our integrative approach shows that the evolution of the signature pheromone structure of olethreutine moths involved the recruitment of a gene resulting from an ancient gene duplication within a Lepidoptera-specific desaturase lineage. Members of other expanded FAD subfamilies do not appear to play a role in chemical communication. This advises for caution when postulating the consequences of lineage-specific expansions based on genomics alone.

Methods

Semiochemical data set

To compile a dataset of the semiochemicals used as sex attractants by tortricid moths, we manually retrieved data from the Pherobase, a freely accessible database of pheromones and semiochemicals which compiles information from the literature and peer-reviewed publications [65]. Altogether, we gathered data for 536 taxa belonging to the family Tortricidae. Specifically, we were interested in collecting information about the chemical structure of the bioactive fatty acyl derivatives, i.e., double-bond position and number and length of the aliphatic chain. Following the database convention, we classified bioactive compounds as pheromone or attractant. The term pheromone denotes a component of the female sex attractant, identified from females, and with demonstrated biological activity; attractant characterizes those bioactive compounds found by male responses alone, such as by systematic field screening or by significant non-target capture in traps baited with pheromone lures. Although these compounds are not necessarily used in the natural communication system, they have in most cases a high probability of being pheromone components. Details are in Additional file 3: Table S1.

Molecular data set

To infer the evolutionary relationship among the species represented in our trait dataset, we used sequence data from 137 terminal taxa in our phylogenetic analysis, including 131 species of Tortricidae as ingroup. The outgroup contained six taxa, each a species representative of Cossidae, Galacticidae, Heliocosmidae, Lacturidae, Limacodidae, and Sesiidae, respectively. Our data set is based largely on Regier et al. [18] and Fagua et al. [19], augmented with a few additional Genbank accessions. Following the aforementioned authors, we used regions of five nuclear genes and one mitochondrial gene. The nuclear genes were carbamoyl-phosphate synthetase II (CAD), dopa decarboxylase (DDC), enolase (Eno), period (PER), and wingless (WG); the mitochondrial gene was the fragment of cytochrome oxidase I (COI) used for barcoding. Details are in Additional file 4: Table S2. For each genus, we used one representative species as the type species. Our main goal was to place the major transition in the pheromone communication system of Tortricidae and thus to infer the phylogenetic relationships with good resolution at the subfamily and tribe level. We retrieved a total of 340 Genbank accessions (see Additional file 4: Table S2 for accession numbers).

Phylogenetic tree

We concatenated the molecular data corresponding to one taxon and obtained an alignment of the entire molecular data set using MAFFT [66] as implemented in Geneious (Biomatters). The final alignment comprised 7591 sites. To define the best scheme of partitions using the MAFFT alignment, we used PartitionFinder v2.1.1 [67]. We tested models of nucleotide substitution with partitions representing codon position and genes using greedy search [68], AICC model selection, and setting the phylogeny program to phyml [69] (default value for all other settings). A 16-partition scheme was selected and used to perform the parametric phylogenetic analysis in IQ-TREE 1.6.11 [70, 71]. The maximum likelihood analysis was performed using default settings and by calculating Shimodaira-Hasegawa-like approximate ratio test (SH-aLRT) support and ultrafast bootstrap support (UFBoot [72];) after 1000 replicates each. We used the R package ggtree [73] for visualizing and annotating the phylogenetic tree with associated semiochemical data.

Insects, tissue collection, and RNA-Seq

Pupae of C. pomonella were purchased from Entomos AG (Switzerland) and were sexed upon arrival to our laboratory. Adult males and females emerged in separate jars in a climate chamber at 23 ± 1 °C under a 16 h:8 h light to dark cycle. On the third day after emergence, pheromone glands of female moths were dissected, immersed in TRIzol Reagent (Thermo Scientific), and stored at − 80 °C until RNA extraction. Two pools of 30 glands were collected 3 h before and after the onset of the photophase, respectively. Total RNA was extracted according to the manufacturer’s instructions. RNA concentration and purity were initially assessed on a NanoDrop2000 (Thermo Fisher Scientific). Library preparation and paired-end Illumina sequencing (2 × 100 bp) were performed by BGI (Hong Kong, PRC). We obtained ~ 66 million QC-passing reads per sample.

Genome annotation improvements

To ensure that all FAD genes were correctly annotated in the C. pomonella genome, we used our RNA-Seq data of pheromone gland to improve the annotation (cpom.ogs.v1.chr.gff3) which we obtained from InsectBase (http://www.insect-genome.com/cydia/). We performed low-quality base trimming and adaptor removal using cutadapt version 1.13 [74] and aligned the trimmed read pairs against the genome using HISAT2 version 2.2.0 [75]. The existing annotation was used to create a list of known splice sites using a python script distributed with HISAT2. We used StringTie version 2.1.3b [76] with the conservative transcript assembly setting to improve the annotation and reconstruct a non-redundant set of transcripts observed in any of the RNA-Seq samples. We applied Trinotate version 3.2.1 [77] to generate a functional annotation of the transcriptome data. We identified candidate FAD genes by searching for genes harboring the PFAM domain corresponding to a fatty acid desaturase type 1 domain (PF00487) in their predicted protein sequences. We used the R package KaryoploteR [78] to visualize the location of the FAD genes.

Desaturase phylogenetic analysis

To investigate the evolutionary relationship of the C. pomonella FADs, we carried out a phylogenetic analysis with other moth FAD proteins. Our sampling includes all the functional classes previously identified in moths (i.e., ∆5, ∆6, ∆9, ∆10, ∆11, ∆14) and representatives of Tortricidae as well as other moth families (see Additional file 5: Table S3 for list of taxa). Protein sequences for these enzymes were downloaded from GenBank. Protein sequences for C. pomonella were predicted from the genome using gffread [79] with our improved annotation.

We aligned amino acid sequences using MAFFT version 7.427 [66, 80]. We used the L-INS-i procedure and the BLOSUM30 matrix as the scoring matrix. The final multiple sequence alignment contained 114 sequences with 805 amino acid sites. The phylogenetic tree was constructed in IQ-TREE 1.6.11 [70]. An automatic model search was performed using ModelFInder [81] with the search restricted to models including the WAG, LG, and JTT substitution models. The maximum likelihood analysis was performed using default settings and by calculating Shimodaira-Hasegawa-like approximate ratio test (SH-aLRT) support and ultrafast bootstrap support (UFBoot [72];) after 1000 replicates each. We used the R package ggtree [73] for visualizing and annotating the phylogenetic tree, with midpoint rooting performed using the midpoint.root function implemented in the phytools package [82].

Transcript abundance estimation

In addition to our pheromone gland data sets, we retrieved 27 already published Illumina paired-end RNA-Seq data sets for C. pomonella available on the Sequence Read Archive (SRA) (see Additional file 6: Table S4 for accession numbers). These data correspond to various tissues or life stages. Following QC, low-quality base trimming, and adaptor removal with cutadapt, read pairs were aligned against the genome using HISAT2. The obtained alignment BAM files were used to estimate transcript abundance using StringTie together with our improved annotation. The abundance tables from StringTie were imported into R using the tximport package [83], which was used to compute gene-level abundance estimates reported as FPKM. We used the R package pheatmap [84] to visualize the expression level of FAD genes.

Cloning of desaturases

First-strand cDNAs were synthesized from 1 μg of total RNA using the SuperScript IV cDNA synthesis kit (Thermo Scientific). The cDNA products were diluted to 50 μL. PCR amplification (50 μL) was performed using 2.5 μL of cDNA as template with attB-flanked gene-specific primers (Additional file 7: Table S5) in a Veriti Thermo Cycler (Thermo Fisher Scientific) using Phusion polymerase master mix (Thermo Fisher Scientific). Cycling parameters were as follows: an initial denaturing step at 94 °C for 5 min, 35 cycles at 96 °C for 15 s, 55 °C for 30 s, 72 °C for 1 min followed by a final extension step at 72 °C for 10 min. Specific PCR products were cloned into the pDONR221 vector (Thermo Fisher Scientific) in the presence of BP clonase (Invitrogen). Plasmid DNAs were isolated according to standard protocols and recombinant plasmids were subjected to sequencing using universal M13 primers and the BigDye terminator cycle sequencing kit v1.1 (Thermo Fisher Scientific). Sequencing products were EDTA/ethanol-precipitated, dissolved in formamide, and loaded for analysis on an in-house capillary 3130xl Genetic analyzer (Applied Biosystems).

Pheromone compounds and corresponding fatty acid precursors

All synthetic pheromone compounds used as references for identification came from our laboratory collection unless otherwise indicated. Pheromone compounds are referred to as short forms, e.g., (E8,E10)-8,10-dodecadienol is E8E10-12:OH, the corresponding acetate would be E8E10-12:OAc, corresponding acyl moiety (dodecadienoate) is E8E10-12:Acyl (although it may occur as CoA-derivative), and the methyl ester thereof is referred to as E8E10-12:Me. The dodecanoic methyl ester (methyl laurate, 12:Me) and tetradecanoic methyl ester (methyl myristate, 14:Me) were purchased from Larodan Fine Chemicals (Malmö, Sweden). The (E)-9-dodecenol (E9-12:OH) was purchased from Pherobank (Wageningen, The Netherlands) and then converted into the corresponding methyl ester (E9-12:Me). E8E10-12:OAc was purchased from Bedoukian (Danbury, CT, USA) and converted to the corresponding alcohol by hydrolysis using a 0.5-M solution of KOH in methanol. Fatty alcohols were oxidized to the corresponding acid with pyridinium dichromate in dimethylformamide as described by Bjostad and Roelofs [85]. The methyl esters were prepared as described under fatty acid analysis (see below).

Heterologous expression of desaturases in yeast

Plasmids containing the genes of interest were selected as entry clones and subcloned into the copper-inducible pYEX-CHT vector [86]. Sequences of recombinant constructs were verified by Sanger sequencing. The pYEX-CHT recombinant expression vectors harboring C. pomonella FADs were introduced into the double deficient elo1∆/ole1∆ strain (MATa elo1::HIS3 ole1::LEU2 ade2 his3 leu2 ura3) of the yeast Saccharomyces cerevisiae, which is defective in both desaturase and elongase functions [87] using the S.c. EasyComp Transformation kit (Thermo Fisher Scientific). For selection of uracil and leucine prototrophs, the transformed yeast was allowed to grow on an SC plate containing 0.7% YNB (w/o AA; with ammonium sulfate) and a complete drop-out medium lacking uracil and leucine (Formedium, UK), 2% glucose, 1% tergitol (type Nonidet NP-40, Sigma), 0.01% adenine (Sigma), and containing 0.5 mM oleic acid (Sigma) as extra fatty acid source. After 4 days at 30 °C, individual colonies were selected and used to inoculate 10 mL selective medium, which was grown at 30 °C at 300 rpm in a shaking incubator for 48 h. Yeast cultures were diluted to an OD600 of 0.4 in 10 mL of fresh selective medium containing 2 mM CuSO4 for induction, with or without addition of a biosynthetic precursor. Yeast cells contain sufficient quantities of naturally occurring palmitic acid; hence, supplementation with 16:Me was typically not necessary. Lauric acid and myristic acid were added to the medium as methyl esters (12:Me and 14:Me). In subsequent assays, the monounsaturated methyl esters E9-12:Me and Z9-12:Me were also supplemented. Each FAME precursor was prepared at a concentration of 100 mM in 96% ethanol and added to reach a final concentration of 0.5 mM in the culture medium. Yeasts were cultured in 30 °C with Cu2+ induction. After 48 h, yeast cells were harvested by centrifugation at 3000 rpm and the supernatant was discarded. The pellets were stored at − 80 °C until fatty acid analysis. Yeast expression experiments were conducted with three independent replicates.

Heterologous expression of desaturases in insect cells

The expression construct for Cpo_CPRQ in the baculovirus expression vector system (BEVS) donor vector pDEST8_CPRQ was made by LR reaction. Recombinant bacmids were made according to instructions for the Bac-to-Bac® Baculovirus expression system given by the manufacturer (Invitrogen) using DH10EMBacY (Geneva Biotech). Baculovirus generation was done using Spodoptera frugiperda Sf9 cells (Thermo Fisher Scientific), Ex-Cell 420 serum-free medium (Sigma), and baculoFECTIN II (OET). The virus was then amplified once to generate a P2 virus stock using Sf9 cells and Ex-Cell 420 medium. Viral titer in the P2 stock was determined using the BaculoQUANT all-in-one qPCR kit (OET) and found to be 3 × 108 pfu/mL for Cpo_CPRQ.

Insect cell lines Sf9 were diluted to 2 × 106 cells/mL. Heterologous expression was performed in 20-mL cultures in Ex-Cell 420 medium and the cells infected at an MOI of 1. The cultures were incubated in 125-mL Erlenmeyer flasks (100 rpm, 27 °C), with fatty acid methyl-ester substrates supplemented at a final concentration of 0.25 mM after 1 day. After 3 days, 7.5-mL samples were taken from the culture and centrifuged for 15 min at 4500g at 4 °C. The pellets were stored at − 80 °C until fatty acid analysis. Sf9 expression experiments were conducted in three replicates.

Fatty acid analysis

Total lipids from cell pellets were extracted using 3.75 mL of methanol/chloroform (2:1 v/v) in clear glass tubes. One milliliter of 0.15 M acetic acid and 1.25 mL of water were added to wash the chloroform phase. Tubes were vortexed vigorously and centrifuged at 2000 rpm for 2 min. The bottom chloroform phase (~ 1 mL), which contains the total lipids, was transferred to a new glass tube. Fatty acid methyl esters (FAMEs) were produced from this total lipid extract. The solvent was evaporated to dryness under gentle nitrogen flow. One milliliter of sulfuric acid 2% in methanol was added to the tube, vortexed vigorously, and incubated at 90 °C for 1 h. After incubation, 1 mL of water was added and mixed well, and 1 mL of hexane was used to extract the FAMEs.

The methyl ester samples were subjected to GC-MS analysis on a Hewlett Packard 6890 GC (Agilent) coupled to a mass selective detector HP 5973 (Agilent). The GC was equipped with an HP-INNOWax column (30 m × 0.25 mm × 0.25 μm; Agilent), and helium was used as carrier gas (average velocity 33 cm/s). The MS was operated in electron impact mode (70 eV), and the injector was configured in splitless mode at 220 °C. 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. The methyl esters were identified based on mass spectra and retention times in comparison with those of our collection of synthetic standards.

Double-bond positions of monounsaturated compounds were confirmed by dimethyl disulfide (DMDS) derivatization [88], followed by GC-MS analysis. FAMEs (50 μL) were transferred to a new tube, 50 μL DMDS were added, and the mixture was incubated at 40 °C overnight in the presence of 5 μL of iodine (5% in diethyl ether) as catalyst. Hexane (200 μL) was added to the sample, and the reaction was neutralized by addition of 50–100 μL sodium thiosulfate (Na2S2O3; 5% in water). The organic phase was recovered and concentrated under a gentle nitrogen stream to a suitable volume [41]. DMDS adducts were analyzed on an Agilent 6890 GC system equipped with a HP-5MS capillary column (30 m × 0.25 mm × 0.25 μm; Agilent) coupled with an HP 5973 mass selective detector. The oven temperature was set at 80 °C for 1 min, raised to 140 °C at a rate of 20 °C/min, then to 250 °C at a rate of 4 °C/min and held for 20 min [36].

Double-bond positions of di-unsaturated compounds with conjugated double bonds were confirmed by derivatization with 4-methyl-1,2,4-triazoline-3,5-dione (MTAD) [89]. MTAD adducts were prepared using 15 μL of the methyl ester extracts transferred into a glass vial. The solvent was evaporated under gentle nitrogen flow, and 10 μL of dichloromethane (CH2Cl2) was added to dissolve the FAMEs. The resulting solution was treated with 10 μL of MTAD (Sigma; 2 mg/mL in CH2Cl2). The reaction was ran for 10 min at room temperature, and 2 μL was subjected to GC-MS analysis on an Agilent 6890 GC system equipped with a HP-5MS capillary column (30 m × 0.25 mm × 0.25 μm; Agilent) coupled with an HP 5973 mass selective detector. The injector temperature was set to 300 °C and the oven was set to 100 °C and then increased by 15 °C/min up to 300 °C followed by a hold at 300 °C for 20 min.

Availability of data and materials

All data generated or analyzed during this study are included in this published article and its supplementary information files.

Abbreviations

FAD:

Fatty acyl desaturase

FAME:

Fatty acid methyl ester

DMDS:

Dimethyl disulfide

MTAD:

4-Methyl-1,2,4-triazoline-3,5-dione

References

  1. 1.

    Lespinet O, Wolf YI, Koonin EV, Aravind L. The role of lineage-specific gene family expansion in the evolution of eukaryotes. Genome Res. 2002;12(7):1048–59. https://doi.org/10.1101/gr.174302.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  2. 2.

    Hahn MW, Bie TD, Stajich JE, Nguyen C, Cristianini N. Estimating the tempo and mode of gene family evolution from comparative genomic data. Genome Res. 2005;15(8):1153–60. https://doi.org/10.1101/gr.3567505.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  3. 3.

    Hahn MW, Han MV, Han S-G. Gene family evolution across 12 Drosophila genomes. PLoS Genet. 2007;3(11):e197. https://doi.org/10.1371/journal.pgen.0030197.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  4. 4.

    Butenandt A, Beckmann R, Stamm D, Hecker E. Uber den Sexual-Lockstoff des Seidenspinners Bombyx mori - Reindarstellung und Konstitution. Z Natforsch B. 1959;14(4):283–4.

    Google Scholar 

  5. 5.

    Löfstedt C, Wahlberg N, Millar JM. Evolutionary patterns of pheromone diversity in Lepidoptera. In: Allison JD, Cardé RT, editors. Pheromone communication in moths: evolution, behavior and application. Berkeley: University of California Press; 2016. p. 43–82.

    Google Scholar 

  6. 6.

    Tillman JA, Seybold SJ, Jurenka RA, Blomquist GJ. Insect pheromones-an overview of biosynthesis and endocrine regulation. Insect Biochem Mol Biol. 1999;29(6):481–514. https://doi.org/10.1016/S0965-1748(99)00016-8.

    CAS  Article  PubMed  Google Scholar 

  7. 7.

    Sperling P, Ternes P, Zank TK, Heinz E. The evolution of desaturases. Prostaglandins Leukot Essent Fatty Acids. 2003;68(2):73–95. https://doi.org/10.1016/S0952-3278(02)00258-2.

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Takahashi A, Tsaur S-C, Coyne JA, Wu C-I. The nucleotide changes governing cuticular hydrocarbon variation and their evolution in Drosophila melanogaster. Proc Natl Acad Sci U S A. 2001;98(7):3920–5. https://doi.org/10.1073/pnas.061465098.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  9. 9.

    Roelofs WL, Liu W, Hao G, Jiao H, Rooney AP, Linn CE. Evolution of moth sex pheromones via ancestral genes. Proc Natl Acad Sci U S A. 2002;99(21):13621–6. https://doi.org/10.1073/pnas.152445399.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Fang S, Ting C-T, Lee C-R, Chu K-H, Wang C-C, Tsaur S-C. Molecular evolution and functional diversification of fatty acid desaturases after recurrent gene duplication in Drosophila. Mol Biol Evol. 2009;26(7):1447–56. https://doi.org/10.1093/molbev/msp057.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  11. 11.

    Shirangi TR, Dufour HD, Williams TM, Carroll SB. Rapid evolution of sex pheromone-producing enzyme expression in Drosophila. PLoS Biol. 2009;7(8):e1000168. https://doi.org/10.1371/journal.pbio.1000168.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  12. 12.

    Albre J, Liénard MA, Sirey TM, Schmidt S, Tooman LK, Carraher C, Greenwood DR, Löfstedt C, Newcomb RD. Sex pheromone evolution is associated with differential regulation of the same desaturase gene in two genera of leafroller moths. PLoS Genet. 2012;8(1):e1002489. https://doi.org/10.1371/journal.pgen.1002489.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  13. 13.

    Roelofs WL, Rooney AP. Molecular genetics and evolution of pheromone biosynthesis in Lepidoptera. Proc Natl Acad Sci U S A. 2003;100(16):9179–84. https://doi.org/10.1073/pnas.1233767100a.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  14. 14.

    Helmkampf M, Cash E, Gadau J. Evolution of the insect desaturase gene family with an emphasis on social hymenoptera. Mol Biol Evol. 2015;32(2):456–71. https://doi.org/10.1093/molbev/msu315.

    Article  PubMed  Google Scholar 

  15. 15.

    Roe AD, Weller SJ, Baixeras J, Brown J, Cummings MP, Davis D, Kawahara AY, Parr C, Regier JC, Rubinoff D. Evolutionary framework for Lepidoptera model systems. Genetics and Molecular Biology of Lepidoptera. Boca Raton: CRC Press; 2009. p. 1–24.

  16. 16.

    Gilligan TM, Baixeras J, Brown JW. T@RTS: online world catalogue of the Tortricidae (Ver. 4.0); 2018.

    Google Scholar 

  17. 17.

    Roelofs WL, Brown RL. Pheromones and evolutionary relationships of Tortricidae. Annu Rev Ecol Syst. 1982;13(1):395–422. https://doi.org/10.1146/annurev.es.13.110182.002143.

    CAS  Article  Google Scholar 

  18. 18.

    Regier JC, Brown JW, Mitter C, Baixeras J, Cho S, Cummings MP, Zwick A. A molecular phylogeny for the leaf-roller moths (Lepidoptera: Tortricidae) and its implications for classification and life history evolution. PLoS One. 2012;7(4):e35574. https://doi.org/10.1371/journal.pone.0035574.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  19. 19.

    Fagua G, Condamine FL, Horak M, Zwick A, Sperling FAH. Diversification shifts in leafroller moths linked to continental colonization and the rise of angiosperms. Cladistics. 2017;33(5):449–66. https://doi.org/10.1111/cla.12185.

    Article  Google Scholar 

  20. 20.

    Liu W, Jiao H, Murray NC, O'Connor M, Roelofs WL. Gene characterized for membrane desaturase that produces (E)-11 isomers of mono- and diunsaturated fatty acids. Proc Natl Acad Sci U S A. 2002;99(2):620–4. https://doi.org/10.1073/pnas.221601498.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Liu W, Jiao H, O’Connor M, Roelofs WL. Moth desaturase characterized that produces both Z and E isomers of Δ11-tetradecenoic acids. Insect Biochem Mol Biol. 2002;32(11):1489–95. https://doi.org/10.1016/S0965-1748(02)00069-3.

  22. 22.

    Hao G, O’Connor M, Liu W, Roelofs WL. Characterization of Z/E11- and Z9-desaturases from the obliquebanded leafroller moth, Choristoneura rosaceana. J Insect Sci. 2002;2(1):26.

    Article  Google Scholar 

  23. 23.

    Witzgall P, Stelinski L, Gut L, Thomson D. Codling moth management and chemical ecology. Annu Rev Entomol. 2008;53(1):503–22. https://doi.org/10.1146/annurev.ento.53.103106.093323.

    CAS  Article  PubMed  Google Scholar 

  24. 24.

    Roelofs W, Comeau A, Hill A, Milicevic G. Sex attractant of the codling moth: characterization with electroantennogram technique. Science. 1971;174(4006):297–9. https://doi.org/10.1126/science.174.4006.297.

    CAS  Article  PubMed  Google Scholar 

  25. 25.

    Beroza M, Bierl BA, Moffitt HR. Sex pheromones: (E,E)-8,10-dodecadien-1-ol in the codling moth. Science. 1974;183(4120):89–90. https://doi.org/10.1126/science.183.4120.89.

    CAS  Article  PubMed  Google Scholar 

  26. 26.

    McDonough LM, Moffitt HR. Sex pheromone of the codling moth. Science. 1974;183(4128):978.

    CAS  Article  Google Scholar 

  27. 27.

    Arn H, Guerin PM, Buser HR, Rauscher S, Mani E. Sex pheromone blend of the codling moth, Cydia pomonella: evidence for a behavioral role of dodecan-1-ol. Experientia. 1985;41(11):1482–4. https://doi.org/10.1007/BF01950048.

    CAS  Article  Google Scholar 

  28. 28.

    Yamaoka R, Taniguchi Y, Hayashiya K. Bombykol biosynthesis from deuterium-labeled (Z)-11-hexadecenoic acid. Experientia. 1984;40(1):80–1. https://doi.org/10.1007/BF01959112.

    CAS  Article  Google Scholar 

  29. 29.

    Tumlinson JH, Teal PEA, Fang N. The integral role of triacyl glycerols in the biosynthesis of the aldehydic sex pheromones of Manduca sexta (L.). Biorg Med Chem. 1996;4(3):451–60. https://doi.org/10.1016/0968-0896(96)00025-9.

    CAS  Article  Google Scholar 

  30. 30.

    Navarro I, Mas E, Fabriàs G, Camps F. Identification and biosynthesis of (E,E)-10,12-tetradecadienyl acetate in Spodoptera littoralis female sex pheromone gland. Biorg Med Chem. 1997;5(7):1267–74. https://doi.org/10.1016/S0968-0896(97)00072-2.

    CAS  Article  Google Scholar 

  31. 31.

    Löfstedt C, Bengtsson M. Sex pheromone biosynthesis of (E,E)-8,10-dodecadienol in codling moth Cydia pomonella involves E9 desaturation. J Chem Ecol. 1988;14(3):903–15. https://doi.org/10.1007/BF01018782.

    Article  PubMed  Google Scholar 

  32. 32.

    Witzgall P, Chambon J-P, Bengtsson M, Unelius CR, Appelgren M, Makranczy G, Muraleedharan N, Reed DW, Hellrigl K, Buser H-R, Hallberg E, Bergström G, Tóth M, Löfstedt C, Löfqvist J. Sex pheromones and attractants in the Eucosmini and Grapholitini (Lepidoptera, Tortricidae). Chemoecology. 1996;7(1):13–23. https://doi.org/10.1007/BF01240633.

  33. 33.

    Wan F, Yin C, Tang R, Chen M, Wu Q, Huang C, Qian W, Rota-Stabelli O, Yang N, Wang S, Wang G, Zhang G, Guo J, Gu L(A), Chen L, Xing L, Xi Y, Liu F, Lin K, Guo M, Liu W, He K, Tian R, Jacquin-Joly E, Franck P, Siegwart M, Ometto L, Anfora G, Blaxter M, Meslin C, Nguyen P, Dalíková M, Marec F, Olivares J, Maugin S, Shen J, Liu J, Guo J, Luo J, Liu B, Fan W, Feng L, Zhao X, Peng X, Wang K, Liu L, Zhan H, Liu W, Shi G, Jiang C, Jin J, Xian X, Lu S, Ye M, Li M, Yang M, Xiong R, Walters JR, Li F. A chromosome-level genome assembly of Cydia pomonella provides insights into chemical ecology and insecticide resistance. Nat Commun. 2019;10(1):4237. https://doi.org/10.1038/s41467-019-12175-9.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  34. 34.

    Knipple DC, Rosenfield C-L, Nielsen R, You KM, Jeong SE. Evolution of the integral membrane desaturase gene family in moths and flies. Genetics. 2002;162(4):1737–52.

    CAS  PubMed  PubMed Central  Google Scholar 

  35. 35.

    Lassance J-M, Löfstedt C. Concerted evolution of male and female display traits in the European corn borer, Ostrinia nubilalis. BMC Biol. 2009;7(1):10.

    Article  Google Scholar 

  36. 36.

    Wang H-L, Liénard MA, Zhao C-H, Wang C-Z, Löfstedt C. Neofunctionalization in an ancestral insect desaturase lineage led to rare Δ6 pheromone signals in the Chinese tussah silkworm. Insect Biochem Mol Biol. 2010;40(10):742–51. https://doi.org/10.1016/j.ibmb.2010.07.009.

  37. 37.

    Hagström ÅK, Albre J, Tooman LK, Thirmawithana AH, Corcoran J, Löfstedt C, Newcomb RD. A novel fatty acyl desaturase from the pheromone glands of Ctenopseustis obliquana and C. herana with specific Z5-desaturase activity on myristic acid. J Chem Ecol. 2014;40(1):63–70. https://doi.org/10.1007/s10886-013-0373-1.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  38. 38.

    Liu W, Rooney AP, Xue B, Roelofs WL. Desaturases from the spotted fireworm moth (Choristoneura parallela) shed light on the evolutionary origins of novel moth sex pheromone desaturases. Gene. 2004;342(2):303–11. https://doi.org/10.1016/j.gene.2004.08.017.

    CAS  Article  PubMed  Google Scholar 

  39. 39.

    Lu F, Wei Z, Luo Y, Guo H, Zhang G, Xia Q, Wang Y. SilkDB 3.0: visualizing and exploring multiple levels of data for silkworm. Nucleic Acids Res. 2020;48(D1):D749–D55. https://doi.org/10.1093/nar/gkz919.

    CAS  Article  PubMed  Google Scholar 

  40. 40.

    Hao G, Liu W, O’Connor M, Roelofs WL. Acyl-CoA Z9- and Z10-desaturase genes from a New Zealand leafroller moth species, Planotortrix octo. Insect Biochem Mol Biol. 2002;32(9):961–6. https://doi.org/10.1016/S0965-1748(01)00176-X.

    CAS  Article  PubMed  Google Scholar 

  41. 41.

    Liénard MA, Strandh M, Hedenström E, Johansson T, Löfstedt C. Key biosynthetic gene subfamily recruited for pheromone production prior to the extensive radiation of Lepidoptera. BMC Evol Biol. 2008;8(1):270. https://doi.org/10.1186/1471-2148-8-270.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  42. 42.

    Liénard MA, Lassance J-M, Wang H-L, Zhao C-H, Piskur J, Johansson T, Löfstedt C. Elucidation of the sex-pheromone biosynthesis producing 5,7-dodecadienes in Dendrolimus punctatus (Lepidoptera: Lasiocampidae) reveals ∆11- and ∆9-desaturases with unusual catalytic properties. Insect Biochem Mol Biol. 2010;40(6):440–52. https://doi.org/10.1016/j.ibmb.2010.04.003.

  43. 43.

    Moto K, Suzuki MG, Hull JJ, Kurata R, Takahashi S, Yamamoto M, Okano K, Imai K, Ando T, Matsumoto S. Involvement of a bifunctional fatty-acyl desaturase in the biosynthesis of the silkmoth, Bombyx mori, sex pheromone. Proc Natl Acad Sci U S A. 2004;101(23):8631–6. https://doi.org/10.1073/pnas.0402056101.

  44. 44.

    Matoušková P, Pichová I, Svatoš A. Functional characterization of a desaturase from the tobacco hornworm moth (Manduca sexta) with bifunctional Z11- and 10,12-desaturase activity. Insect Biochem Mol Biol. 2007;37(6):601–10. https://doi.org/10.1016/j.ibmb.2007.03.004.

    CAS  Article  PubMed  Google Scholar 

  45. 45.

    Witzgall P, Bengtsson M, Rauscher S, Liblikas I, Bäckman A-C, Coracini M, Anderson P, Löfqvist J. Identification of further sex pheromone synergists in the codling moth, Cydia pomonella. Entomol Exp Appl. 2001;101(2):131–41. https://doi.org/10.1046/j.1570-7458.2001.00898.x.

  46. 46.

    Hill J, Rastas P, Hornett EA, Neethiraj R, Clark N, Morehouse N, Celorio-Mancera MdlP, Cols JC, Dircksen H, Meslin C, Keehnen N, Pruisscher P, Sikkink K, Vives M, Vogel H, Wiklund C, Woronik A, Boggs CL, Nylin S, Wheat CW. Unprecedented reorganization of holocentric chromosomes provides insights into the enigma of lepidopteran chromosome evolution. Sci Adv. 2019;5(6):eaau3648.

  47. 47.

    Kawahara AY, Plotkin D, Espeland M, Meusemann K, Toussaint EFA, Donath A, Gimnich F, Frandsen PB, Zwick A, dos Reis M, Barber JR, Peters RS, Liu S, Zhou X, Mayer C, Podsiadlowski L, Storer C, Yack JE, Misof B, Breinholt JW. Phylogenomics reveals the evolutionary timing and pattern of butterflies and moths. Proc Natl Acad Sci U S A. 2019;116(45):22657–63. https://doi.org/10.1073/pnas.1907847116.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  48. 48.

    Bielawski JP, Yang Z. Maximum likelihood methods for detecting adaptive evolution after gene duplication. In: Meyer A, Van de Peer Y, editors. Genome evolution: gene and genome duplications and the origin of novel gene functions. Dordrecht: Springer Netherlands; 2003. p. 201–12. https://doi.org/10.1007/978-94-010-0263-9_20.

    Chapter  Google Scholar 

  49. 49.

    Taylor JS, Raes J. Duplication and divergence: the evolution of new genes and old ideas. Annu Rev Genet. 2004;38(1):615–43. https://doi.org/10.1146/annurev.genet.38.072902.092831.

    CAS  Article  PubMed  Google Scholar 

  50. 50.

    Xia Y-H, Zhang Y-N, Ding B-J, Wang H-L, Löfstedt C. Multi-functional desaturases in two Spodoptera moths with ∆11 and ∆12 desaturation activities. J Chem Ecol. 2019;45(4):378–87. https://doi.org/10.1007/s10886-019-01067-3.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  51. 51.

    Serra M, Piña B, Bujons J, Camps F, Fabriàs G. Biosynthesis of 10,12-dienoic fatty acids by a bifunctional Δ11desaturase in Spodoptera littoralis. Insect Biochem Mol Biol. 2006;36(8):634–41. https://doi.org/10.1016/j.ibmb.2006.05.005.

  52. 52.

    Jurenka RA. Biochemistry of female moth sex pheromones. In: Blomquist GJ, Vogt R, editors. Insect pheromone biochemistry and molecular biology. The biosynthesis and detection of pheromones and plant volatiles; 2003. p. 53–80.

    Chapter  Google Scholar 

  53. 53.

    Liénard MA, Wang H-L, Lassance J-M, Löfstedt C. Sex pheromone biosynthetic pathways are conserved between moths and the butterfly Bicyclus anynana. Nat Commun. 2014;5(1):3957. https://doi.org/10.1038/ncomms4957.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Liénard MA, Löfsted C. Functional flexibility as a prelude to signal diversity? Role of a fatty acyl reductase in moth pheromone evolution. Commun Integr Biol. 2010;3(6):586–8. https://doi.org/10.4161/cib.3.6.13177.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Emelianov I, Marec F, Mallet J. Genomic evidence for divergence with gene flow in host races of the larch budmoth. Proc R Soc Lond B Biol Sci. 2004;271(1534):97–105. https://doi.org/10.1098/rspb.2003.2574.

    Article  Google Scholar 

  56. 56.

    Wakamura S. Identification of sex-pheromone components of the podborer, Matsumuraeses falcana (Walshingham) (Lepidoptera: Tortricidae). Appl Entomol Zool. 1985;20(2):189–98. https://doi.org/10.1303/aez.20.189.

    CAS  Article  Google Scholar 

  57. 57.

    Frerot B, Priesner E, Gallois M. A sex attractant for the green budworm moth, Hedya nubiferana. Z Naturforsch C J Biosci. 1979;34(12):1248–52. https://doi.org/10.1515/znc-1979-1229.

    Article  Google Scholar 

  58. 58.

    Foster SP, Roelofs WL. Sex pheromone differences in populations of the brownheaded leafroller, Ctenopseustis obliquana. J Chem Ecol. 1987;13(3):623–9. https://doi.org/10.1007/BF01880104.

    CAS  Article  PubMed  Google Scholar 

  59. 59.

    Witzgall P, Kirsch P, Cork A. Sex pheromones and their impact on pest management. J Chem Ecol. 2010;36(1):80–100. https://doi.org/10.1007/s10886-009-9737-y.

    CAS  Article  PubMed  Google Scholar 

  60. 60.

    Löfstedt C, Xia Y-H. Biological production of insect pheromones in cell and plant factories. In: Blomquist GJ, Vogt RG, editors. Insect pheromone biochemistry and molecular biology. 2nd ed. London: Academic; 2021. p. 89–121. https://doi.org/10.1016/B978-0-12-819628-1.00003-1.

    Chapter  Google Scholar 

  61. 61.

    Nešněrová P, Šebek P, Macek T, Svatoš A. First semi-synthetic preparation of sex pheromones. Green Chem. 2004;6(7):305–7. https://doi.org/10.1039/B406814A.

    Article  Google Scholar 

  62. 62.

    Hagström ÅK, Wang H-L, Liénard MA, Lassance J-M, Johansson T, Löfstedt C. A moth pheromone brewery: production of (Z)-11-hexadecenol by heterologous co-expression of two biosynthetic genes from a noctuid moth in a yeast cell factory. Microb Cell Factories. 2013;12(1):125. https://doi.org/10.1186/1475-2859-12-125.

    CAS  Article  Google Scholar 

  63. 63.

    Ding B-J, Hofvander P, Wang H-L, Durrett TP, Stymne S, Löfstedt C. A plant factory for moth pheromone production. Nat Commun. 2014;5(1):3353. https://doi.org/10.1038/ncomms4353.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  64. 64.

    Holkenbrink C, Ding B-J, Wang H-L, Dam MI, Petkevicius K, Kildegaard KR, Wenning L, Sinkwitz C, Lorántfy B, Koutsoumpeli E, França L, Pires M, Bernardi C, Urrutia W, Mafra-Neto A, Ferreira BS, Raptopoulos D, Konstantopoulou M, Löfstedt C, Borodina I. Production of moth sex pheromones for pest control by yeast fermentation. Metab Eng. 2020;62:312–21. https://doi.org/10.1016/j.ymben.2020.10.001.

    CAS  Article  PubMed  Google Scholar 

  65. 65.

    El-Sayed AM. The Pherobase: database of pheromones and semiochemicals; 2020. https://www.pherobase.com/.

  66. 66.

    Katoh K, Misawa K, Kuma K-I, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30(14):3059–66. https://doi.org/10.1093/nar/gkf436.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  67. 67.

    Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B. PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol Biol Evol. 2017;34(3):772–3. https://doi.org/10.1093/molbev/msw260.

    CAS  Article  Google Scholar 

  68. 68.

    Lanfear R, Calcott B, Ho SYW, Guindon S. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol Biol Evol. 2012;29(6):1695–701. https://doi.org/10.1093/molbev/mss020.

    CAS  Article  PubMed  Google Scholar 

  69. 69.

    Guindon S, Dufayard J-F, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010;59(3):307–21. https://doi.org/10.1093/sysbio/syq010.

    CAS  Article  PubMed  Google Scholar 

  70. 70.

    Nguyen L-T, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32(1):268–74. https://doi.org/10.1093/molbev/msu300.

    CAS  Article  Google Scholar 

  71. 71.

    Chernomor O, von Haeseler A, Minh BQ. Terrace aware data structure for phylogenomic inference from supermatrices. Syst Biol. 2016;65(6):997–1008. https://doi.org/10.1093/sysbio/syw037.

    Article  PubMed  PubMed Central  Google Scholar 

  72. 72.

    Hoang DT, Chernomor O, von Haeseler A, Minh BQ, Vinh LS. UFBoot2: improving the ultrafast bootstrap approximation. Mol Biol Evol. 2018;35(2):518–22. https://doi.org/10.1093/molbev/msx281.

    CAS  Article  PubMed  Google Scholar 

  73. 73.

    Yu G, Smith DK, Zhu H, Guan Y, Lam TT-Y. ggtree: an r package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol. 2017;8(1):28–36. https://doi.org/10.1111/2041-210X.12628.

    Article  Google Scholar 

  74. 74.

    Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011;17(1):10–2. https://doi.org/10.14806/ej.17.1.200.

    Article  Google Scholar 

  75. 75.

    Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. 2019;37(8):907–15. https://doi.org/10.1038/s41587-019-0201-4.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  76. 76.

    Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5. https://doi.org/10.1038/nbt.3122.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  77. 77.

    Bryant DM, Johnson K, DiTommaso T, Tickle T, Couger MB, Payzin-Dogru D, Lee TJ, Leigh ND, Kuo TH, Davis FG, Bateman J, Bryant S, Guzikowski AR, Tsai SL, Coyne S, Ye WW, Freeman RM Jr, Peshkin L, Tabin CJ, Regev A, Haas BJ, Whited JL. A tissue-mapped axolotl de novo transcriptome enables identification of limb regeneration factors. Cell Rep. 2017;18(3):762–76. https://doi.org/10.1016/j.celrep.2016.12.063.

  78. 78.

    Gel B, Serra E. karyoploteR: an R/Bioconductor package to plot customizable genomes displaying arbitrary data. Bioinformatics. 2017;33(19):3088–90. https://doi.org/10.1093/bioinformatics/btx346.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  79. 79.

    Pertea G, Pertea M. GFF utilities: GffRead and GffCompare. F1000Res. 2020;9:304.

    Article  Google Scholar 

  80. 80.

    Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30(4):772–80. https://doi.org/10.1093/molbev/mst010.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  81. 81.

    Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat Methods. 2017;14(6):587–9. https://doi.org/10.1038/nmeth.4285.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  82. 82.

    Revell LJ. phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol Evol. 2012;3(2):217–23. https://doi.org/10.1111/j.2041-210X.2011.00169.x.

    Article  Google Scholar 

  83. 83.

    Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Res. 2016;4:1521.

    Article  Google Scholar 

  84. 84.

    Kolde R. pheatmap: pretty heatmaps. 1.0.12 ed. 2019.

  85. 85.

    Bjostad LB, Roelofs WL. Sex pheromone biosynthetic precursors in Bombyx mori. Insect Biochem Mol Biol. 1984;14:275–8. https://doi.org/10.1016/0020-1790(84)90060-X.

  86. 86.

    Patel O, Fernley R, Macreadie I. Saccharomyces cerevisiae expression vectors with thrombin-cleavable N- and C-terminal 6x(His) tags. Biotechnol Lett. 2003;25(4):331–4. https://doi.org/10.1023/A:1022384828795.

    CAS  Article  PubMed  Google Scholar 

  87. 87.

    Schneiter R, Tatzer V, Gogg G, Leitner E, Kohlwein SD. Elo1p-dependent carboxy-terminal elongation of C14:1Δ9 to C16:1Δ11 fatty acids in Saccharomyces cerevisiae. J Bacteriol. 2000;182(13):3655–60. https://doi.org/10.1128/JB.182.13.3655-3660.2000.

  88. 88.

    Buser HR, Arn H, Guerin P, Rauscher S. Determination of double bond position in mono-unsaturated acetates by mass spectrometry of dimethyl disulfide adducts. Anal Chem. 1983;55(6):818–22. https://doi.org/10.1021/ac00257a003.

    CAS  Article  Google Scholar 

  89. 89.

    Marques FA, Millar JG, McElfresh JS. Efficient method to locate double bond positions in conjugated trienes. J Chromatogr A. 2004;1048(1):59–65. https://doi.org/10.1016/S0021-9673(04)01147-1.

    CAS  Article  PubMed  Google Scholar 

Download references

Acknowledgements

We thank Wendell Roelofs for his comments on an earlier version of this manuscript and the Lund Protein Production Platform (LP3) for technical support and advice on the expression in insect cell lines. The computations in this paper were run on the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University.

Funding

This project has received funding from the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 760798, Olefine) (CL), the Swedish Foundation for Strategic Research (grant no. RBP 14–0037, Oil Crops for the Future) (CL), the Swedish Research Council (CL), Formas (CL), and the Royal Physiographic Society in Lund (JML). Open Access funding provided by Lund University.

Author information

Affiliations

Authors

Contributions

JML, BJD, and CL designed the research; BJD collected tissue and processed samples for RNA-Seq; JML performed analyses related to phylogenetics and bioinformatics; BJD performed functional assays and chemical analyses; JML prepared the figures; JML wrote the paper with input from BJD and CL. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Christer Löfstedt.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors declare the following competing interests: CL is a founder of SemioPlant AB; BJD and CL have filed patents related to the use of Cpo_CPRQ for pheromone production. JML declares that he has no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

Phylogeny of lepidoptera FAD genes. Extended version of the maximum likelihood tree displayed in Fig. 4. The tree was obtained for predicted amino acid sequence of 114 FAD genes (805 aligned positions) of 28 species, with branch support values calculated from 1000 replicates using the Shimodaira-Hasegawa-like approximate ratio test (SH_aLRT) and ultrafast bootstrapping (UFboot). Support values for branches are indicated by colored circles, with color assigned based SH-aLRT and UFBoot supports with 80% and 95% as thresholds of branch selection for SH-aLRT and UFBoot supports, respectively. The major constituent six subfamilies of First Desaturase (A1 to E) and two subfamilies of Front-End (Cyt-b5-r) and Sphingolipid Desaturases (Ifc), respectively, are indicated following the nomenclature proposed by Helmkampf et al. (2015). For First Desaturases, the different shades correspond to the indicated putative biochemical activities and consensus signature motif (if any). Triangles indicate sequences from C. pomonella. The scale bar represents 0.5 substitutions per amino acid position. Species are indicated by three- or four-letter prefixes (see Additional file 5: Table S3 for details). Biochemical activities (or signature motif) are indicated after the abbreviated species name, followed by accession number in parenthesis.

Additional file 2: Figure S2.

Functional characterization of desaturase activity of First-desaturases. Total ion chromatograms of fatty acid methyl ester (FAME) products of Cu2+-induced ole1 elo1 S. cerevisae yeast supplemented with saturated acyl precursors and transformed with (A & C) empty expression vector (control), (B) pYEX-CHT-Cpo_NPVE, (D) pYEX-CHT-Cpo_KPSE, and (E) pYEX-CHT-Cpo_GATD.

Additional file 3: Table S1.

List of terminal taxa and corresponding bioactivity data used in this study.

Additional file 4: Table S2.

List of terminal taxa and corresponding GenBank accession numbers used in this study.

Additional file 5: Table S3.

List of species represented in the FAD phylogeny.

Additional file 6: Table S4.

List of SRA accessions corresponding to RNA-Seq data sets for Cydia pomonella.

Additional file 7: Table S5.

Primer sequences used for amplification and cloning of desaturase genes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Lassance, JM., Ding, BJ. & Löfstedt, C. Evolution of the codling moth pheromone via an ancient gene duplication. BMC Biol 19, 83 (2021). https://doi.org/10.1186/s12915-021-01001-8

Download citation

Keywords

  • Fatty acyl desaturase
  • Gene family evolution
  • Bifunctional
  • Conjugated double bond
  • Tortricidae