Identification and Expression Profiles of Sex Pheromone Biosynthesis and Transport Related Genes in Spodoptera litura

Although the general pathway of sex pheromone synthesis in moth species has been established, the molecular mechanisms remain poorly understood. The common cutworm Spodoptera litura is an important agricultural pest worldwide and causes huge economic losses annually. The female sex pheromone of S. litura comprises Z9,E11-14:OAc, Z9,E12-14:OAc, Z9-14:OAc, and E11-14:OAc. By sequencing and analyzing the transcriptomic data of the sex pheromone glands, we identified 94 candidate genes related to pheromone biosynthesis (55 genes) or chemoreception (39 genes). Gene expression patterns and phylogenetic analysis revealed that two desaturase genes (SlitDes5 and SlitDes11) and one fatty acyl reductase gene (SlitFAR3) showed pheromone gland (PG) biased or specific expression, and clustered with genes known to be involved in pheromone synthesis in other moth species. Furthermore, 4 chemoreception related genes (SlitOBP6, SlitOBP11, SlitCSP3, and SlitCSP14) also showed higher expression in the PG, and could be additional candidate genes involved in sex pheromone transport. This study provides the first solid background information that should facilitate further elucidation of sex pheromone biosynthesis and transport, and indicates potential targets to disrupt sexual communication in S. litura for a novel pest management strategy.


Introduction
Species-specific sex pheromone-elicited behaviors play a key role in sexual communication and reproduction in most moth species, which therefore serve as a good model to study reproductive isolation in animals, from insects to mammals [1][2][3]. Moth sex pheromones are biosynthesized and released by specialized sex pheromone glands (PGs) that are located along the intersegmental membrane between the 8 th and 9 th abdominal segments of females [4,5]. In most moths, sex pheromone components are composed of C10 -C18 unsaturated acyclic aliphatic compounds with a functional group such as formyl, hydroxyl, or acyloxyl [6,7].
Sex pheromone biosynthesis in moths begins with a palmitic or stearic acid moiety that is synthesized de novo in the PG through modifications of the fatty acid biosynthetic pathway [8]. Through a series of enzymatic reactions such as desaturation, chain-shortening reaction, reduction, acetylation, and oxidation, the palmitic or stearic acids are then converted to the final pheromone components in a step-wise manner [5,9,10]. Therefore, different enzymes are likely to be involved in the different reactions, and to date, the genes encoding 4 different classes of enzymes that are essential for this pathway have been functionally identified-desaturases (Des), fatty acid reductases (FARs), fatty acid transport proteins (FATPs), and acyl-CoaA-binding proteins (ACBPs). Among these, Des proteins are the most intensively studied class of enzymes involved in moth sex pheromone biosynthesis, which can introduce double bonds into pheromone precursors. Previous studies have demonstrated the broad functional diversity of these enzymes such as Δ5 [11,12], Δ6 [13], Δ9 [9,14], Δ11 [15,16], Δ10-12 [10], and Δ14 desaturase [17]. FARs are responsible for reducing fatty acids to alcohols, and have also been functionally identified in a few moth species, including Bombyx mori [18], Ostrinia nubilalis [19] and 4 Heliothine species [20]. FATPs and ACBPs have been functionally identified to play roles in the production of the B. mori sex pheromone bombykol based on in vivo RNA interference methods [21,22]. In addition, some other important enzymes that have not been functionally confirmed are postulated to be involved in the sex pheromone biosynthesis pathway. For example, biochemical studies have suggested that acetyltransferase (ACT) and alcohol oxidase (AO) play a role by converting alcohol into acetate ester [23] and by oxidizing alcohol into the corresponding aldehyde component [13,24], respectively.
The common cutworm Spodoptera litura (Lepidoptera: Noctuidae) is an important agricultural pest worldwide and causes huge economic losses annually. The female sex pheromones of S. litura have been identified as a blend of Z9-E11-14:OAc, Z9-E12-14:OAc, Z9-14:OAc, and E11-14:OAc with a ratio of 100:27:20:27 in China [25]. To date, only 9 chemoreception genes have been identified and functionally characterized, including 5 odorant-binding proteins (OBPs) [26,27] and 4 pheromone receptors (PRs) [28]. However, the genes involved in the pheromone biosynthesis of S. litura have not been explored. In present study, we constructed a genetic database of the genes expressed in the female PGs of S. litura using the Illumina HiSeq (TM) 2500 sequencing platform. In total, we identified 94 genes that are possibly related to pheromone biosynthesis (55 genes) or chemoreception (39 genes). Furthermore, tissue expression evaluation and phylogenetic analyses were performed to postulate the functions of the identified genes. The results indicated that some of these genes might play crucial roles in the biosynthesis and transport of S. litura sex pheromones and could be as candidates for further functional studies.

Transcriptome Sequencing and Sequence Assembly
We carried out a next-generation sequencing analysis using a cDNA library constructed from the female PGs of S. litura using the Illumina HiSeq (TM) 2500 platform. The transcriptome sequencing provided approximately 63 million reads (6.3 Gb). After clustering and redundancy filtering, we finally acquired 42,646 unigenes with an N50 length of 2,093 bp (Table 1). We called these 42,646 ones unigenses according to some recently published papers [29,30], although each of them may not necessarily represents a unique gene. Of the 42,646 unigenes, those with a sequence length of more than 500 bp accounted for 40.08% of the total transcriptome assembly (Fig 1).
GO annotation was used to classify the transcripts into functional groups according to GO categories. Of the 42,646 unigenes, 11,692 (27.41%) could be annotated based on sequence homology. In the molecular function category, the genes expressed in the PG were mostly enriched to binding, catalytic activity and transporter activity. Cellular and metabolic processes were the most highly represented in the biological process categories, and cell, cell part, and organelle were most abundantly represented in the cellular component category (Fig 2). In addition, the Des gene SlitDes5 was the most abundant of all unigenes, and two CSP genes (SlitCSP3 and SlitCSP2) also showed very high abundance (Fig 3 and S1 Table).
Of the 94 identified genes (Tables 3 and 4), the sequences of 5 genes were identical to those already deposited in GenBank

Expression Profile of the Putative Genes Related to Sex Pheromone Biosynthesis and Transport
To investigate the general expression profiles of the candidate genes, reverse transcriptionpolymerase chain reaction (RT-PCR) analyses were conducted for all 94 genes (Figs 4 and 5), and the expression levels of 16 selected genes were further quantified with qPCR ( Fig 6) to validate the RT-PCR results. The overall relative expression profiles of these genes in different tissues obtained were similar with the two methods. Most of the pheromone biosynthesis-related genes (>70%) displayed PG-biased or specific expression, and 4 chemoreception related genes (SlitOBP6, SlitOBP11, SlitCSP3, and SlitCSP14) were also more highly expressed in the PG than in other tissues of S. litura (Table 2).
Desaturase (Des) and Fatty Acyl Reductase (FAR). The results of RT-PCR and qPCR showed that 9 of the 12 Des genes showed PG-biased or specific expression, which is greater than the proportion reported for other moth species (Table 2); the other Des genes, SlitDes7, SlitDes9, and SlitDes12, were detected in both the PG and other parts of the body ( Acyl-CoA Binding Protein (ACBP), Fatty Acid Transport Protein (FATP) and Acetyl-CoA Carboxylase (ACC). Of the 2 SlitACBP and 3 SlitFATP genes, only SlitACBP2 and Slit-FATP3 were highly expressed in the PG, whereas the SlitACC gene encoding an ACC with 94% identity to the ACC of Agrotis segetum (GenBank Accession No.: AID66639.1) did not show PG-biased or specific expression ( Table 2, Figs 4D and 6).
Odorant Binding Protein (OBP) and Chemosensory Protein (CSP). We identified a total of 25 OBP genes from the PG of S. litura in this study, including 1 PBP (OBP25) and 24 other OBPs. Of the 25 SlitOBPs, 4 genes (SlitOBP8, SlitOBP9, SlitOBP13, and SlitOBP14) displayed a very wide range of tissue distribution in all 4 tissues examined, whereas 8 other genes (SlitOBP3, SlitOBP4, SlitOBP5, SlitOBP10, SlitOBP16, SlitOBP19, SlitOBP24, and SlitOBP25) were expressed predominately or specifically in the adult antennae. Furthermore, 2 genes (Sli-tOBP6 and SlitOBP11) were expressed at much higher levels in the PG than in other tissues.
Unlike the SlitOBPs, most of the SlitCSPs were expressed at similar levels among the 4 tissues examined; only SlitCSP3 and SlitCSP14 showed significantly higher expression in the PG than in other tissues (
The OBP phylogenetic tree showed that SlitOBP25 (SlitPBP3) was clustered into the PBP/ GOBP clade. The other SlitOBPs and all SlitCSPs were clustered with at least one lepidopteran orthologous gene (S1 and S2 Figs).

Discussion
In the present study, we sequenced and analyzed the transcriptome of the PG of S. litura. Among the 42,646 unigenes identified, only 27.41% could be annotated to one or more GO Agrotis ipsilon   terms, which is similar to other lepidopteran species [37][38][39], indicating that a large number of S. litura genes are either non-coding or are homologous with genes that do not have any GO term. Importantly, we identified 89 novel genes that are involved in sex pheromone biosynthesis and transport in S. litura. Our results not only provide an important foundation for further elucidation of the molecular mechanisms of sex pheromone metabolism but also provide general insight into insect physiology and development of a novel pest control strategy [40].
S. litura has 4 sex pheromone components: Z9,E11-14:OAc, Z9,E12-14:OAc, Z9-14:OAc, and E11-14:OAc (at a ratio of 100:27:20:27) in China [25]. According to some previous studies related to the biosynthesis pathway of Δ11or Δ9-containing sex pheromones in several moth species, the defined pathway involves a step of Δ11 or Δ9 desaturation that is catalyzed by Δ11 or Δ9 desaturase, respectively [16,[41][42][43][44]. On the other hand, Liu et al. [45] reported that a Δ9 desaturase participates in the production of the Δ11-containing pheromone component by introducing a Δ9-double bond at 14:CoA, followed by carbon chain elongation to Δ11-16: CoA. In sum, these findings indicate that the Δ11 and Δ9 desaturases are likely responsible for introduction of the Δ11 and Δ9-double bonds in all 4 pheromones of S. litura. In this study, we obtained 9 desaturase genes showing PG-biased or specific expression. Further phylogenetic analysis showed that SlitDes5 was clearly assigned to the Δ11 desaturase group, and was most closely related to Z/E11 of S. littoralis (GenBank Accession No. Q6US81). SlitDes9 and SlitDes11 were allocated to the Δ9 (18C>16C) and Δ9 (16C>18C) desaturase groups, respectively, but only SlitDes11 displayed a PG-biased expression pattern; the other genes were allocated to other desaturase groups. Therefore, SlitDes5 and SlitDes11 are very likely involved in the desaturation step from saturated acids (14C) to unsaturated acids, with a double bond introduced at the 11 th and 9 th positions of the carbon chain, respectively. To date, there is no report on Δ12 desaturase genes in moths; therefore, the type of enzyme that S. litura uses to introduce a double bond at the 12 th position needs to be further studied.
In the process of sex pheromone biosynthesis, once the specific unsaturated fatty acid precursors are produced, they will be converted into corresponding alcohols by FAR, which has been demonstrated in different moth species [18][19][20]46]. In this study, we identified 11 FAR genes with PG-biased or specific expression, but only SlitFAR3 was clustered in the moth pgFAR group, suggesting that this gene plays a crucial role in the biosynthesis of precursor alcohols. ACTs are essential in the biosynthesis of acyloxyl components, although no related gene had been functionally characterized to date. In the present study, 15 ACT genes were found to be highly or specifically expressed in the PG by RT-PCR and qPCR, indicating that these genes may participate in the progress of ACT to produce corresponding acetate ester. FATPs have been functionally confirmed to bind to and transport fatty acids across the insect hemolymph into PG cells for pheromone biosynthesis in B. mori [22] and Eilema japonica [47], and ACBPs have been functionally confirmed to serve as carriers or cellular deposits for the acyl-CoAs used in pheromone biosynthesis [22]. Therefore, the PG-biased expression of SlitFATP3 and SlitACBP2 may indicate that they play a similar role in the pheromone biosynthesis of S. litura. Some previous studies identified the presence of chemosensilla on the ovipositor [48,49], which may function in the chemoreception of plant odors, ovipositor-deterring pheromones, and sex pheromones, suggesting that female moths may receive and transport pheromone compounds or their precursors via their ovipositor [35,[50][51][52]. This finding further suggests that there may be a feedback loop in the moth's PG (including the ovipositor) to achieve accurate control of the biosynthesis pathway and release of sex pheromones. To date, many studies have demonstrated the role of OBPs and CSPs in the binding and transportation of hydrophobic molecules, including plant volatiles, sex pheromones and their precursors [50,[53][54][55][56][57]; therefore, these two proteins appear to be essential for the sex pheromone biosynthesis pathway. Similarly to other moths [32,34,35], we identified a total of 39 genes in the S. litura PG, including 25 OBP and 14 CSP genes, but only 4 genes (SlitOBP6, SlitOBP11, SlitCSP3, and SlitCSP14) showed PG-biased expression, indicating that these genes may play important roles in the binding and transport of sex pheromone compounds and plant volatiles.
In conclusion, through sequencing and transcriptome analyses, we obtained an extensive set of putative genes that may be related to the biosynthesis and transport of the sex pheromone of S. litura. As the first step towards understanding the functions of these genes, we conducted a comprehensive and comparative examination of gene expression patterns and conducted a phylogenetic analysis with sequences of other species. We identified a number of genes with PG-biased or specific expression, indicating their involvement in the biosynthesis and transport of the sex pheromone. Further studies are needed to explore the functions of these genes with integrated functional studies.

Insects Rearing
S. litura were reared on an artificial diet comprising wheat germ flour and soybean flour, and were sexed as pupae and kept separately in cages for eclosion. The rearing conditions were 27°C, with a 14-h light:10-h dark photoperiod, and 65 ± 5% relative humidity. Adults were provided with a cotton swab dipped in 10% honey solution that was renewed daily.

Tissue Collection
For transcriptome sequencing, 20-25 PGs (with the ovipositor) were collected from 3-day-old virgin female adults at 6-7 h into the scotophase, since 3-day-old moths show particularly high mating activity [25]. For the tissue expression study, 25-30 female antennae (FA), 25-30 male antennae(MA), 20-25 PGs (with the ovipositor), 10-15 whole insect body without PGs(B), and 10-15 whole insect body without PGs and antennae (Bo) were also collected under the same conditions. All samples were immediately frozen and stored at −70°C until use.

cDNA Library Construction
Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). cDNA library construction and Illumina sequencing of the samples were performed at Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). The mRNA was purified from 3 μg of total RNA using oligo (dT) magnetic beads and fragmented into short sequences in the presence of divalent cations at 94°C for 5 min. Then, the first-strand cDNA was generated using random hexamer-primed reverse transcription, followed by synthesis of the second-strand cDNA using RNaseH and DNA polymerase I. After end repair and ligation of adaptors, the products were amplified by PCR and purified using the QIAquick PCR Purification Kit to create a cDNA library, and the library quality was assessed on the Agilent Bioanalyzer 2100 system.

Clustering and Sequencing
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq (TM) 2500 platform, and paired-end reads were generated.

De novo Assembly of Short Reads and Gene Annotation
Clean short reads were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from the raw reads. Transcriptome de novo assembly was carried out with these short reads in the assembling program Trinity (r20140413p1) [58,59] with min_ kmer_cov set to 2 by default, and all other parameters at the default settings. The resulting sequences of Trinity were deemed to be unigenes. Unigenes larger than 150 bp were first aligned to protein databases by BlASTX, including Nr, Swiss-Prot, KEGG, and COG (evalue < 10 −5 ), to retrieve proteins with the highest sequence similarity to the obtained unigenes along with their protein functional annotations. Then, we used the Blast2GO program [60] to get the GO annotation of the unigenes, and the GO functional classification was obtained using WEGO software [61].

Expression Abundance Analysis of the Unigenes
The expression abundance of these unigenes were calculated based on the reads per kilobase per million mapped reads (RPKM) method [62], using the formula: RPKM (A) = (10,00,000 × C × 1,000)/(N × L), where RPKM (A) is the abundance of gene A, C is the number of reads that uniquely aligned to gene A, N is the total number of reads that uniquely aligned to all genes, and L is the number of bases in gene A. The RPKM method is able to eliminate the influence of different gene lengths and sequencing discrepancy in the calculation of expression abundance.

RNA Isolation and cDNA Synthesis
Total RNA was extracted with the SV 96 Total RNA Isolation System (Promega, Madison, WI, USA) following the manufacturer's instructions, in which DNaseI digestion was included to avoid genomic DNA contamination. RNA quality was checked with a spectrophotometer (NanoDropTM 1000, Thermo Fisher Scientific, USA). The single-stranded cDNA templates were synthesized from 1 μg of total RNA from various tissue samples using the PrimeScript RT Master Mix (TaKaRa, Dalian, China).

Phylogenetic Analyses
The phylogenetic trees were constructed for phylogenetic analyses of SlitOBPs, SlitCSPs, SlitDes and SlitFAR, based on these genes (the signal peptides of sequences had been removed of the putative chemosensory genes) and the sequences of other insects. The OBP dataset contained 25 sequences from S. litura, 19 from Manduca sexta [37,64], 15 from S. littoralis [64], 23 from Sesamia inferens [65], and 43 from B. mori [66]. The CSP dataset contained 14 sequences from S. litura, 14 from M. sexta [37], 9 from S. littoralis [64], 13 from S. inferens [65], and 14 from B. mori [67]. The Des dataset contained 12 sequences from S. litura and 59 from other insects [33]. The FAR dataset contained 13 sequences from S. litura and 56 from other insects [33]. The names and accession numbers of the genes used for phylogenetic tree construction are listed in S2 Table. Amino acid sequences were aligned with ClustalX 1.83 [68], and unrooted trees were constructed by MEGA5.0 [69] using the neighbor-joining method, with Poisson correction of distances. Node support was assessed using a bootstrap procedure base on 1000 replicates.

Reverse Transcription-PCR Analyses
Gene-specific primers across the ORFs of predicted chemosensory genes were designed using Primer Premier 5.0 (PREMIER Biosoft International, CA, USA). The sequences of these primers are listed in S3 Table. PCR experiments including negative controls (no cDNA template) were carried out under the following conditions: 94°C for 4 min, 30-35 cycles at 94°C for 30 sec, 60°C for 30 sec, and 72°C for 40 sec, and final extension for 10 min at 72°C. The reactions were performed in a total volume of 25 μL, containing 12.5 μL of 2×EasyTaq PCR SuperMix (TransGene, Beijing, China), 0.4 μM for each primer, 1 μL of sample cDNA (15 ng/μL), 9.5 μL of sterilized H 2 O. PCR products were analyzed by electrophoresis on 1.5% w/v agarose gel in TAE buffer (40 mM Tris-acetate, 2 mM Na 2 EDTAÁH 2 O). The gene encoding S. litura glyceraldehyde-3-phosphate dehydrogenase (SlitGAPDH) (GenBank Accession No.: HQ012003.2) was used as a reference gene for checking the integrity of the cDNA template. Each PCR reaction was done at least twice.

Quantitative Real Time-PCR Validation
The expression profile of 16 putative genes with PG-enriched or specific expression was evaluated to validate the accuracy of the RT-PCR results using quantitative real time-PCR (qPCR) experiments. The qPCR was performed on an ABI 7300 (Applied Biosystems, Foster City, CA, USA) using a mixture of 10 μL 2× SYBR Green PCR Master Mix, 0.4 μL each primer (10 μM), 2.5 ng of sample cDNA, and 6.8 μL sterilized ultrapure H 2 O. The reaction programs were at 95°C for 30 sec, followed by 40 cycles of 95°C for 5 sec, and 60°C for 31 sec. The results were analyzed using the ABI 7300 analysis software SDS 1.4. The qPCR primers (S3 Table) were designed using Beacon Designer 7.9 (PREMIER Biosoft International, CA, USA). The mRNA levels were measured by qPCR using SYBR Premix ExTaq (TaKaRa, Dalian, Liaoning, China). Subsequently, fluorescence was measured throughout a 55-95°C melting curve in order to detect a single gene-specific peak and to check the absence of primer dimer peaks; single and discrete peaks were detected for all primers tested. Negative controls were non-template reactions (replacing cDNA with H 2 O).
The expression levels of 16 genes were calculated relative to the reference gene SlitGAPDH (GenBank Accession No.: HQ012003.2) and SlitEF (elongation factor-1 alpha) (GenBank Accession No.: DQ192234.1) using the Q-Gene method in the Microsoft Excel-based software Visual Basic [70,71] For each sample, three biological replications were performed with each biological replication measured in three technique replications.

Statistical Analyses
Data (mean ± SE) form various samples were subjected to one-way nested analysis of variance (ANOVA) followed by the least significant difference test (LSD) for mean comparison, and two-sample analysis was performed by the Student t-test using SPSS Statistics 17.0 software (SPSS Inc., Chicago, IL, USA).  Table. The tree was constructed with MEGA5.0, using the neighbour-joining method. Values at the nodes are results of bootstrap with 1000 replicates. (TIF) S1