Lipases and carboxylesterases are involved in interspecic pheromone differences between two moth species

Moth sex pheromones are a classical model for studying sexual selection. Females produce a species-specic pheromone blend that attracts males. Revealing the enzymes involved in the interspecic variation in blend composition is key for understanding the evolution of these sexual communication systems. The nature of the enzymes involved in the variation of acetate esters, which are prominent compounds in moth pheromone blends, remains unclear. We identied enzymes involved in acetate metabolism in two closely related species: Heliothis (Chloridea) subexa and H. (C.) virescens, which differ in production of acetate esters. Through comparative transcriptomic analyses and CRISPR/Cas9 knockouts, we showed that two lipases and two esterases induce lower levels of acetate esters in female pheromones. To place our ndings in an evolutionary context, we explored the molecular evolution of related lipases and esterases in Lepidoptera. Together, our results show that lipases and carboxylesterases are unexpectedly involved in tuning Lepidoptera pheromones composition. analyses, to identify new genes underlying a phenotype. Through these analyses, we revealed that lipases and carboxylesterases are involved in moth sex pheromone metabolism, specically in the degradation of acetate esters. As the role of these enzyme families in moth sex pheromone biosynthesis was previously unexpected, this paves the way for new biochemical and molecular analyses to gain insights in the specic functions of these enzymes in moth sex pheromone biosynthesis. Our study also highlights that pheromone degradation needs to be considered as an important mechanism in the evolution of moth sex pheromones.


Introduction
One of the most fascinating questions in evolutionary biology is to understand how new species arise.
Sexual signals, through which individuals recognize suitable mating partners, are critical for prezygotic isolation of species. As such, sexual signals offer a seductive model to uncover evolutionary patterns that lead to speciation. Identifying genes underlying the production and detection of sexual signals opens a door to the analysis of their evolution across species.
As many other species, nocturnal moths use sex pheromones for mate recognition 1 . Female moths produce a species-speci c blend of a small number of pheromone components to which conspeci c males are attracted 2 . As moths are one of the most diverse groups of animals, and the sex pheromone of more than 2000 species have been identi ed 3 , moths are exemplary to understand the formation of new species through the evolution of sexual communication.
The species-speci city of moth sex pheromone blends is determined by the ratio of alcohols, aldehydes and acetate esters (hereafter also referred to as "acetates") with carbon backbones of various lengths and degrees of desaturation 2,4 . Biochemical studies have unraveled some of the biosynthetic pathways leading to these pheromone compounds 4,5 . Before emission, alcohols can be acetylated to acetate esters by acetyltransferases, or oxidized to aldehydes by alcohol oxidases. Acetate esterase or aldehyde reductase activities in and around the gland can convert (back) the acetate ester or aldehyde pheromone compounds into alcohols 6,7 . The terminal steps in (de)acetylation of pheromone components determine the nal species-speci c pheromone composition. Surprisingly, even though the sex pheromone blend of almost all moth species contain acetate esters, the genes underlying acetate production remains elusive despite their pivotal role 8 .
To identify the gene(s) underlying acetate production, we conducted several genetic analyses with two closely related sympatric moth occurring species. These noctuid moths have recently been reclassi ed from the paraphyletic genus Heliothis to the monophyletic genus Chloridea 9 but we retain the older name here for continuity with the extensive literature on pheromones in these species. The most striking pheromone difference between these two species is the presence of (Z)-11-hexadecenyl acetate (Z11- 16:OAc) and two other acetate esters (Z7-16:OAc and Z9-16:OAc) in the H. sub exa female pheromone, which are absent in the H. virescens female blend. This difference plays a critical role in the reproductive isolation of these two species, because the presence of these acetates in the pheromone blend inhibits the attraction of H. virescens males while increasing H. sub exa male attraction [10][11][12] . Since these two species can also be hybridized in captivity, they are an excellent model to identify the genes underlying acetate production.
Previously, we identi ed two quantitative trait loci (QTL) underlying the difference in acetate production [13][14][15] . In this study, we introgressed one of these QTLs from H. virescens corresponding to one end of Chromosome 20 into the genomic background of H. sub exa. We then analyzed the transcriptomes of the sex pheromone glands of the introgressed and pure-strain females, followed by RT-qPCR. With these analyses, we found that certain lipases and esterases encoded within the introgressed QTL region were overexpressed in H. virescens and in the introgressed line compared to H. sub exa. Knock-outs of these genes through CRISPR/Cas9 correspondingly caused an increase in acetate amounts in the introgressed line, implicating these enzymes in hydrolysis of acetate esters in the pheromone gland. This reveals an unexpected critical role of speci c lipases and esterases in the evolution of moth sex pheromones.

Dominant effect of H.virescens QTL for reduced acetate level
Through continuous backcrossing one of the two major QTLs for acetate production 13-15 , we isolated this QTL (2.4cM at one end of Chromosome 20) into an otherwise complete H. sub exa genomic background. We refer to this introgressed line as DD23. Homozygous (VV) and heterozygous (VS) DD23 females showed signi cantly lower acetate levels than wild-type H. sub exa (SS) (Tukey post-hoc test, P < 0.001) ( Figure 1). Also, the introgression of Chr20-QTL from H. virescens into H. sub exa has a dominant effect, as VS females had similarly low levels of acetates as VV (Tukey post-hoc test, p = 0.585).

Differential expression of lipases putatively acting on esters in the DD23 transcriptome
To determine which genes were differentially expressed in the sex pheromone glands of the different types of females, we analyzed the transcriptomes of these glands after building a reference transcriptome (see Supplementary Results). We considered differentially expressed contigs when we found an absolute and signi cant log2 fold change in expression. This rst ltering revealed 89 differentially expressed contigs, none of which were annotated by Blast2Go with putative functions in pheromone biosynthesis (see Supplementary Table S1). However, among these differentially expressed contigs, ve were annotated as lipases (38088_c0_seq4, 38088_c0_seq5, 38579_c0_seq3, 39873_c1_seq2 and 39873_c1_seq3) and one as putative esterase (comp39435_c1_seq17). As the GO terms predicted that the genes represented by these contigs are putatively involved in acetate degradation, we considered these contigs as candidate transcripts for explaining the phenotypic difference in acetate production between VV and SS.
When we veri ed the differential expression levels in the RNAseq experiment with quantitative real-time PCR experiments on pheromone gland cDNA, we observed higher expression levels of LipX and Est1 in H. virescens and VV individuals compared to H. sub exa and SS females, and this difference was signi cant in comparing VV with H. sub exa and SS individuals (Games-Howell post-hoc test, p < 0.05, Supplementary Table S2

Functional characterization of our candidate genes
To determine the role of the identi ed candidate genes in pheromone biosynthesis, we performed loss-offunction studies by generating mutant DD23 lines using CRISPR/Cas9. To avoid functional compensation between closely related enzymes, we targeted multiple genes at once in each knock-out line. Since we found two esterases in tandem in the QTL region in a draft version of the H. sub exa genome (data not shown), we knocked out both these esterases as well. In the esterase-only knock-out line We did not nd signi cant differences between the ve genotypes in the total amount of pheromone produced, except between VS and the esterase knock-out line, which is most probably due to one VS outlier (Supplementary Figure S3.A). In comparing differences in the three other components that affect male response (Z9-16:Ald, Z11-16:Ald, and Z11-16:OH), we found that Z9-16:Ald levels were different only between VV and SS (Supplementary Figure S3.D). However, the lipase knock-out line as well as the fourgene knock-out line had similar levels of Z11-16:Ald and Z11-16:OH as SS females, but signi cantly lower levels than VS females, while, the esterase knock-out line had similar levels as VS females and signi cantly higher levels than SS females (Supplementary Figure S3.E and S3.F). These results thus show that these lipases and to a lesser extent the esterases are involved in moth sex pheromone acetate production.

Phylogeny of insect acidic lipases
To investigate the evolution of the acidic lipase gene family, to which LipX and LipZ belong, we built a phylogeny of these genes, by exploring the transcriptomes and genomes of 16 insects (including nine Lepidoptera; Figure 4). With this analysis, we discovered that the insect acidic lipases are split across multiple order-speci c supported clades (bootstrap ≥ 70). This is consistent with what has been observed in earlier work on the phylogeny of insect acidic lipases 16 . In this phylogeny, we observed ten Lepidoptera-speci c supported clades annotated from A-J ( Figure 4). Four of these clades include sequences from the DD23 H. sub exa pheromone gland transcriptome.

Phylogeny of LipX and LipZ orthologues in Lepidoptera and analyses of selection pressures
LipX and LipZ both belong to the C clade, which can be divided into three clades encompassing sequences from multiple Lepidoptera super-families (clades C.1 to C.3). These three clades show recent duplication events, as they contain genus or species-speci c sub-clades ( Figure 5.A). To determine selective pressures in the Lepidoptera lipase sequences belonging to the C clade, we built a separate phylogenetic tree using only Lepidoptera sequences ( Figure 5.A). We detected positive selection at the basal node of the C.1 and C.2 subclades (ω > 1, likelihood ratio tests < 0.05, see Supplementary table S4, Supplementary Figure S1). We also observed multiple other instances of positive selection in the tree, most of them are linked with supported clades encompassing genus, or species-speci c duplications.

Phylogeny of Est1 and Est2 orthologues in Lepidoptera and analyses of selection pressures
To study selective pressures in Est1 and Est2 orthologues, we built a phylogenetic tree using esterase sequences from multiple Lepidoptera and used Tsubota and Shiotsuki's carboxylesterase sequences as outliers (see Tsubota and Shiotsuki, 2010a, 2010b) ( Figure 6.A). We observed multiple events of positive selection (ω > 1, likelihood ratio tests p < 0.05, see Supplementary Table S5, Supplementary Figure S2) across the tree. Interestingly, we found a Noctuidae-speci c clade composed of Est1 and Est2 homologues that are clustered in two separated sub-clades. The Est1 clade encompasses genomic sequences from T. ni (XP_026728424.1) and S. litura (XP_022818360.1), and transcriptomic sequences from H. armigera (ATJ44487.1) and H. assulta (ATJ44554.1) pheromone glands, as well as from of S. littoralis (ACV60233) and Sesamia inferens (AII21982) male antennae. The Est2 clade includes genomic sequences from T. ni (XP_026728423.1), S. litura (XP_022818359) and H. armigera (XP_021196515.1), as well as transcripts from Sesamia nonagrioides (ABH01082) and S. littoralis (ACV60239) male antennae. We found evidence of positive selection at the base of the Est1 and Est2 clades and along the evolution of the Est1 H. armigera and H. assulta homologues.

Discussion
Enzymes degrading acetates inHeliothis pheromone gland By conducting RNAseq on female sex pheromone glands from two closely related species with acetate esters (H. sub exa) or without them (H. virescens) as well as acetate QTL introgression lines (VV and SS), we identi ed genes involved in moth sex pheromone acetate production. Surprisingly, we found that Heliothis females with low acetate levels over-express transcripts of two lipases and an esterase, indicating that acetate esters can be degraded in the species without acetates. This result was unexpected, as the general assumption is that species producing acetate esters possess an acetyltransferase that non-acetate-producing species lack; so, there should be no need to degrade an acetate ester that is never made. Using a gene inactivation approach, we demonstrated that lipases and esterases indeed reduced the amount of acetate esters in the pheromone of Heliothis females. These observations could be explained by a model where acidic lipase and esterase enzymes hydrolyze the acetate pheromone compounds that are produced in the gland of both species.
Acetates are produced but degraded in the non-acetate producing species Even though our results and model seem counter-intuitive, as not producing a compound seems more parsimonious than degrading it, the presence of an acetate ester hydrolysis activity in the pheromone gland of H. virescens and Choristoneura fumiferana (Tortricidae) had been suggested before [21][22][23] . We can think of three possible evolutionary scenarios that would have led to the selection of acetate hydrolysis enzymatic activities in pheromone glands in Heliothis spp. or other moths: i) A trace acetates removing scenario, where acetate hydrolysis might counterbalance the spontaneous acetylation of pheromone compounds thus ensuring that acetates are absent in the pheromone, ii) The anti-basal acetate production scenario where acetate enzymes leading to production of acetate pheromones might be involved in other physiological enzymatic pathways and cannot be under-expressed without undesirable side-effects thus leading to signi cant production of acetate pheromone compounds and making degradation of those compounds mandatory, iii) The ne tuning of pheromone composition scenario where acetate hydrolysis would serve to ne-tune the species-speci c relative amounts of acetate-esters, alcohols and aldehydes in the pheromone blend. Below, we discuss the implications of gene families we highlighted in the evolution of moth sex pheromone biosynthesis.
Lipases and pheromone biosynthesis So far, lipases have been rarely studied in relation to moth pheromone metabolism. The function of lipases has been mainly considered in digestion, following work in vertebrates. For example, the human gastric lipase is the most well-known acidic lipase and predominantly hydrolyzes triglycerides and cholesterol esters 24 . The expression of gastric-type and pancreatic-type lipases is dependent on diet in the lighbrown apple moth 25 . Lipases in general are known to act on water-insoluble lipids 26 which is consistent with the involvement of lipases in pheromone hydrolysis. Furthermore, there is a body of evidence supporting that lipases are involved in moth pheromone biosynthesis: 1) Lipases have been shown to be involved in the biosynthesis of long carbon chain pheromones in other insects 27,28 . 2) In moths, the hormone PBAN activates lipase activity in the pheromone gland in some Apoditrysia (Lepidoptera) species 29,30 . 3) Lipase enzymes have been used in hydrolysis reactions of Lymantriinae moth pheromone precursors 31,32 . 4) Du et al. highlighted the expression of some lipases, including three acidic lipases, in the female pheromone glands of Bombyx mori 33 . Using RNAi, they found that two of the acidic lipases increase the proportion of bombykol in the female pheromone 33 . Thus, even though our results are unexpected, there are precedents for the involvement of acidic lipases in Lepidoptera sex pheromone biosynthesis.

Duplication of LipX and LipZ and acetate production in Noctuidae
Our phylogenetic study revealed two orthologues of LipX and LipZ in some Heliothinae moths, but not in the closely related noctuid Spodoptera species. Genomic data support a scenario where the two lipase genes result from a duplication event that happened in a common ancestor of H. virescens and H. armigera. This duplication may have allowed a faster evolution of LipZ compared to LipX and positive selection speci cally in LipZ.
In A. segetum, S. exigua and H. assulta, only the homologue of LipX is expressed, which matches the expression pro le observed in H. sub exa. Interestingly, all these species display acetates in their sex pheromones 3 . However, in S. Litura no homologue of LipX is present in the pheromone gland transcriptome, while the female sex pheromone blend contains a signi cant amount of acetates 34,35 . Together, these ndings suggest that expression of LipX alone is not su cient to hydrolyze the acetates in female sex pheromone glands. Intriguingly, LipZ is found in two Heliothinae species without acetates in their pheromone blend, H. virescens and H. armigera, while their close relatives with acetates (H. sub exa and H. assulta, respectively) do not appear to express a LipZ homologue in their pheromone gland 36,37 . Thus, LipZ expression correlates with the absence of acetates in Noctuidae pheromone glands and vice versa. This supports that the function of LipZ evolved following the evolutionary "trace acetates removing" scenario or "anti-basal acetate production" scenario. Having information on the expression of HarmLipZ in the female pheromone gland could further support this correlation, unfortunately, Li et al.
didn't report the expression of lipases in their work 37 .

LipX homologues and pheromone biosynthesis in other Lepidoptera
While LipX doesn't appear su cient to hydrolyze acetates in the pheromone gland of noctuid species, two acidic lipases of the C.1 clade have been found expressed in Lepidoptera pheromone gland transcriptomes: Px001710 and GBXH01009737.1, which respectively belong to Plutella xylostella (Plutellidae) 38 and Cadra cautella (Pyralidae) 39 . As both these species have acetate esters in their pheromone blend, these expression patterns indicate a potential role of LipX Lepidoptera orthologues in pheromone metabolism. Moreover, in our phylogeny of the C clade of acidic lipases, we observed duplications in multiple super-families of Lepidoptera. These observations suggest an adaptive role of these acidic lipases in ne-tuning moth sex pheromone ratios ( ne tuning of pheromone composition scenario) or to counterbalance unavoidable production of acetate pheromone compounds (anti-basal acetate production evolutionary scenario).

Acidic lipases and pheromone biosynthesis in Lepidoptera
The numerous supported clades of acidic lipases encompassing lipases from one order at a time underlines a diversi cation of the acidic lipase family during insects' evolution. The high number of Lepidoptera acidic lipases clades in the phylogeny supports that the diversity of these enzymes was already high in the basal species of Lepidoptera, suggesting a diversity of functions early in their evolution, which makes it plausible that of them might be involved in pheromone biosynthesis.

Esterases and pheromone biosynthesis in Lepidoptera
In contrast to lipases, carboxylesterases have been extensively studied for their potential role in degradation of acetate esters in the pheromone gland and elsewhere, mainly in male antennae [40][41][42][43][44][45] . Only recently, carboxylesterases have been found to possibly also play a role in pheromone biosynthesis 46 .
Thus, carboxylesterases may not only be involved in the degradation of pheromone components in the context of olfactory reception, but also in the production of Lepidoptera sex pheromone blends.
Our phylogenetic analyses revealed the presence of orthologues of Est1 and Est2 in the Noctuidae clade, except for H. assulta, S. inferens and S. nonagrioides for which we found only one sequence, namely either Est1 or Est2. As the sequences from these three species were obtained from transcriptomic data, it is possible that only one copy is expressed, or one duplicated gene may have been lost. Based on these ndings, the most parsimonious explanation is that a duplication event happened at the basal node of the Noctuidae clade, leading to positive selective pressure in both Est1 and Est2 clades. Since we could not nd expression patterns of these esterases in female pheromone glands that are consistent with acetate esters in the female sex pheromone 37 , the role of each esterases in moth sex pheromone biosynthesis remains elusive. As with LipX homologues in Lepidoptera, this favors that the function of Est1 and Est2 evolved following the " ne tuning of pheromone composition" evolutionary scenarios and "anti-basal acetate production" evolutionary scenarios.

More enzymes involved in acetate synthesis in Heliothis virescens
We and previous workers 23 have never found any trace of acetate esters in H. virescens female pheromone glands. The hypothetical acetyltransferase has eluded all previous attempts at discovery 8 , including our own. Our previous nding of two QTLs underlying acetate production 47 indicates that additional enzymes are involved in acetate production, perhaps in the second QTL. Unfortunately, even after 10 years of effort, we have never been able to introgress this second QTL region from H. virescens into H. sub exa and we are currently exploring different approaches to investigate the presence of potential hydrolyzing genes in this QTL and possible interaction effects.

Conclusion
Our study underlines the effectiveness of introgressing a major QTL into the genomic background of a parental species to conduct transcriptome and loss-of-function analyses, to identify new genes underlying a phenotype. Through these analyses, we revealed that lipases and carboxylesterases are involved in moth sex pheromone metabolism, speci cally in the degradation of acetate esters. As the role of these enzyme families in moth sex pheromone biosynthesis was previously unexpected, this paves the way for new biochemical and molecular analyses to gain insights in the speci c functions of these enzymes in moth sex pheromone biosynthesis. Our study also highlights that pheromone degradation needs to be considered as an important mechanism in the evolution of moth sex pheromones.

Insects rearing
Moths were lab-reared as described in detail in 48 . The so-called DD23 strain was generated by hybridizing H. virescens and H. sub exa, followed by continuous backcrossing to H. sub exa, as described in detail in [13][14][15] . This resulted in the introgression of Chr20-QTL from H. virescens into a complete H. sub exa genomic background, which led to signi cantly lower acetate levels in the pheromone blend 13 . In this study, we used animals homozygous (VV), heterozygous (VS) and wild-type (SS) for this QTL. As this QTL is dominant, i.e. one V-copy is enough to reduce the acetate levels 13 (Figure 1), we used VS females in the loss-of-function analyses. Distinction between the three genotypes was made through PCR ampli cation of marker sequence where an EcoRI restriction site is present in V-allele and absent in the Sallele.

Pheromone extraction
For all pheromone analyses in this study, we used 2 to 4 days old virgin female moths. Females, whose glands were used in the RNAseq and the real-time quantitative PCR, were injected with pheromone biosynthesis activating neuropeptide (PBAN) 1-2 hours before gland extraction to stimulate pheromone production, following 11,15,49 . To minimize possible side-effects of PBAN injections in our subsequent functional analysis experiments, pheromone extractions were performed only on calling females, 3-4 hours into scotophase. Pheromone gland extractions were conducted as described in detail in 48 50 ) were calculated relative to our 200 ng pentadecane internal standard and divided by the total amount to get the relative amount. To break data interdependency and to correct for individual variation in the amount of pheromone produced that play a role in male response (see 51 ), the total amount of the three acetate-esters, Z9-16:Ald, Z11-16:Ald and Z11-16:OH were divided by the amount of 14:Ald. We chose 14:Ald as divisor, because of its not affecting male response behavior 50 and was present in comparable amounts in all extracts, except SS and the esterase knock-out line (Supplementary gure S1.b and Supplementary table S3). Log transformation was performed when this corrected residuals non-normality. Samples with less than a total pheromone amount of 40 ng were excluded, because the ratios of the minor compounds could not be accurately quanti ed in the chromatograms in these samples.

DD23 strain pheromone gland transcriptome assembly
For our transcriptome analysis, we extracted RNA from the pheromone extracted pheromone glands of 10 females, ve of VV genotype and ve of SS, using Ambion™ TRIzol™ Reagent (Thermo Fisher Scienti c, Waltham, Massachusetts, USA) and Direct-zol™ RNA MiniPrep Kit (Zymo Research, Irvine, California, USA). The RNA was sequenced using Illumina HiSeq2500 by BaseClear (Leiden, The Netherlands). Raw reads of used to assemble the pheromone gland transcriptome were deposited on NCBI Sequence Read Archive under the bioproject number PRJNA493752. Resulting reads were trimmed of adapters and for a Phred quality score of 28 or more using Trimmomatic 52 . We then removed ribosomal reads using Ribopicker 53 . We also removed reads without a mate pair. The remaining reads were normalized and de novo assembled using Trinity 54 . We then aligned the Trinity output with the R1 reads used for assembly with Bowtie2 55,56 and removed all contigs matching less than ve reads. The resulting contigs will be further referred to as the DD23 pheromone gland un ltered transcriptome. To improve the quality of the transcriptome, we sequentially used three tools as described in 57 and summarized here. We rst ltered out the contigs without an ORF of 100 amino or more using Transdecoder (http://transdecoder.github.io/). We then used RSEM 58 to remove all contigs with an TPM expression value lesser than one. Finally, we clustered and removed redundant contigs using CD-HIT EST (90% similarity threshold and word size of eight) 59,60 . To assess this reference transcriptome quality, we used BUSCO 61 and Transrate 62 .

Differential expression analysis and differentially expressed transcripts annotation
Differences in expression levels were explored by aligning the reads used for the assembly onto the reference transcriptome, using RSEM 58 , after which we conducted differential expression analysis, using Deseq2 62 and edgeR 63 tools. We selected contigs differentially expressed between SS and VV genotypes by ltering contigs with: a ±2 log2 fold change in expression value (Deseq2) and an associated adjusted P-value (from DEseq2 and edgeR) of 0.05 or lower. We annotated these differentially expressed transcripts by blasting them against the NR NCBI database, using blastx, and attributed GO terms using Blast2GO.

Sequencing of candidate transcripts and quantitative realtime PCR
To con rm our in sillico differential expression results, we rst had to con rm the sequences of the candidate transcripts, which we did by extracting RNA from H. virescens and H. sub exa female pheromone glands, as described above. RNA was retro-transcribed using Verso cDNA Synthesis Kit (Thermo Fisher Scienti c), the candidate gene transcripts were ampli ed using DreamTaq™ DNA Polymerase (Thermo Fisher Scienti c) and primers were designed to amplify an approximately 1000 nucleotides region on cDNA from H. virescens and H. sub exa pheromone glands (see Supplementary  Table S6). Amplicons were sequenced using Sanger sequencing by Marcogen (Amsterdam, the Netherlands). Sequences resulting from genomic DNA Sanger sequencing are available on Genbank, accession numbers OK556469 to OK556476. From these sequences, we designed primers for real-time PCR to get approximately 100 bp amplicons from the alleles of both species for each target gene.
To assess the in vivo expression of the candidate genes, we extracted RNA as described above from 5-6 pheromone-extracted pheromone glands of H. sub exa, H. virescens, SS and VV females. We used 20 µl reaction of 5x HOT FIREPol® EvaGreen® qPCR Mix Plus (ROX) (Solis BioDyne, Tartu, Estonia), 2 ng of cDNA and 1 µM of primer couples for the target genes and used RPS18 as a reference gene. For primer used, see Supplementary Table S6. All biological replicates were also technically replicated. qPCR reactions and measurements were made using an Applied Biosystem 7500 Real-Time PCR System (Thermo Fisher Scienti c).
Knock-out of the candidate genes using CRISPR/Cas9 To functionally characterize the candidate genes discovered in the previous steps, CRISPR/Cas9 knockout experiments were performed as follows. The IDT system for CRISPR/Cas9 was used (Integrated DNA Technologies). Gene-speci c crRNAs are annealed with tracrRNAs and incubated with Cas9 enzyme to form ribonucleotide particles for injection into eggs. Gene-speci c guide RNAs were designed for the two esterases (VVest1-T1 and VVest1-T2 for est1, VVest2-T1 and VVest2-T2) and two lipases (VVLipX-T1 and VVLipX-T2 for LipX, VVLipZ-T1 and VVLipZ-T2 for LipZ, all sequences are in Supplementary Table S6). First, VV eggs were injected with Est1-Est2 guide RNA to knock-out the two esterases. In a second step, VV and VV-Est1-Est2-ko eggs were injected with LipX-LipZ guide RNAs to obtain two lipase knock-outs and four gene knock-outs, respectively. All eggs were 30 min-1 hour old and micro-injected using a Femtojet with a solution previously loaded in home-made glass needles. The solution was prepared by rst dissolving each four scRNAs (Supplementary Table S6) in 1 nmol of tracrRNA to get a nal concentration on 50 µM, after which 200 pmol of this scRNA+tracrRNA mix was combined with 100 pmol of IDT Alt-R Cas9.
The gauze with the injected embryos was kept in Petri-dishes at lab temperature (20°C) and checked daily. After three days, newly hatched neonates were separated in individual cups with wheat germ/soy our-based diet (BioServ Inc., Newark, DE, USA). To identify mutations in the genes, we performed DNA extractions on newly emerged adults by soaking one foreleg in 30 µl of 10% Chelex and 2.5 µl of Proteinase K during three hours at 56 degrees C, then 8 min at 98 degrees C and by nally freezing the mixture before usage (Adapted protocol from Biorad). Distinction between the genotypes was made through PCR ampli cation with gene-speci c primers: VVest1-scr-Fd and VVest1-scr-R for est1, VVest2scr-F and VVest2-scr-R for est2, VVLipX-scr-Fd and VVLipX-scr-R for LipX, and VVLipZ-scr-F and VVLipZ- CDS nucleotide sequences of the proteins used for building the tree were obtained from NCBI NR database, genome databases or the sequencing data of this study. These sequences where used to generate a nucleotide alignment based on the alignment used for the phylogenetic tree of lipases or esterases, using PAL2NAL 68 . This nucleotide alignment and the phylogenetic tree were used as input for selective pressure analysis. We used codeML to infer episodic change in natural selection acting on protein sequence 69 . We used branch site models with default parameters and following the procedure suggested by the user manual of codeML, for each tested clade, we calculated a posterior probability using the procedure described in 70 .

Declarations
Competing interests   Maximum-likelihood phylogeny of insect acidic lipases.